RNA-Seq analysis of differentially expressed genes of Staphylococcus epidermidis isolated from postoperative endophthalmitis and the healthy conjunctiva

Staphylococcus epidermidis (S. epidermidis) is one of the primary pathogens in postoperative endophthalmitis, which is a devastating complication of cataract surgery and often results in irreversible visual loss and even blindness. Meanwhile, it is the most frequently isolated commensal bacterium in the healthy conjunctiva. In this study, we investigated the differentially expressed genes (DEGs) of S. epidermidis isolated from the patients with postoperative endophthalmitis and the healthy conjunctiva to predict their functions and pathways by Illumina high-throughput RNA sequencing. Using genome-wide transcriptional analysis, 281 genes (142 upregulated and 139 downregulated genes) were found to be differentially expressed (fold change ≥ 2, p ≤ 0.05) in the strains from endophthalmitis. Ten randomly selected DEGs were further validated by quantitative reverse transcription polymerase chain reaction (qRT-PCR). GO enrichment analysis suggested that more DEGs were associated with the thioredoxin system and iron ion metabolism. KEGG pathway analysis revealed that more DEGs were associated with the pathways of the two-component system and pyruvate metabolism. Moreover, the gene SE1634 code for staphylococcal toxin was significantly upregulated in S. epidermidis strains of the endophthalmitis, which might be directly responsible for the pathogenesis of endophthalmitis. In conclusion, this research is helpful for further investigations on genes or pathways related with the pathogenesis and therapeutic targets of S. epidermidis endophthalmitis.

www.nature.com/scientificreports/ All clinical samples were inoculated onto blood agar plates (Auto Biotechnology, Zhengzhou, China) and incubated at 37 °C with 5% carbon dioxide for 24 h (Thermo, Massachusetts, USA). S. epidermidis isolates were considered when the incubation exhibited single colony morphology, Gram stain was positive, Catalase test was positive, and Oxidase test was negative. Then the ten pure strains of S. epidermidis were identified by an automatic microbiological identification and susceptibility analysis system (Beckman Coulter WalkAway-96 plus, California USA), which uses colorimetry and fluorescence to identify both gram-negative and gram-positive organisms through a variety of biochemical reactions.
The study was approved by the ethics committee of Shandong Eye Institute. Based on the tenets of the Declaration of Helsinki, written informed consents were obtained from every enrolled participant.

RNA extraction.
For the ten strains of S. epidermidis, the phenol water (Sinopharm Chemical Reagent Co., Shanghai, China) was added and shaken vigorously to pyrolyze each sample. All centrifugal tubes were incubated with an oscillating metal bath at 65 °C for a maximum speed shock for 30 to 60 min. After 5 min of cooling on ice, centrifugation was performed at 12,000×g, 4 °C for 10 min, and then the supernatant water was transferred into a new centrifugal tube. The purification of RNA was performed using the TRK-1002 (Lianchuan Bio, Hangzhou, China) and following the manufacturer's recommendation. The quality of each extracted RNA was checked using the Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA, USA) with an RNA integrity number (RIN) > 7.0 for cDNA library preparation, and the required RNA samples were stored in a refrigerator at − 80 °C before they were sent to Lianchuan Bio for RNA-Seq. cDNA library construction and sequencing. Five-μg RNA of each sample was used to build a cDNA library. Ribosomal RNA (rRNA) was depleted using the Ribo-Zero Magnetic Kit (Epicentre, Madison, WI, USA) before the library was generated using the Illumina Truseq RNA Sample Preparation Kit (Illumina, San Diego, CA, USA) according to the manufacturer's protocol. On the basis of the cleaved RNA fragments, cDNA was synthesized using Uracil-N-Glycosylase (UNG), and a cDNA cluster was generated according to a cBot User Guide to PCR amplification after ligation of adaptors. Illumina sequencing was carried out on an Illumina HiSeq™ 2000 platform (Lianchuan Bio).

RNA-Seq data analysis.
After sequencing, image data was transformed into raw reads and stored in FASTQ format for per sample. The resulting clean reads were obtained by removing the low-quality, adapter, poly-N containing, and shorter-than-70 bp reads, and mapped to the reference genome of S. epidermidis strain ATCC12228 (https ://www.ncbi.nlm.nih.gov/genom e/155?Genomeassembly_id = 299299) using software Rockhopper 26 . Gene expression was normalized by calculating Reads per Kilobase per Million Mapped Reads (RPKM) as described by Mortazavi and colleagues 27 , which could eliminate the effects of gene length and the differences of sequencing quantity when calculating the gene expression. DEGs were determined using the edgeR package of R software 28 , and the false discovery rate ≤ 0.05 and the absolute value of log 2 fold change ≥ 1 (|log 2 FC|≥ 1) were used as the threshold to determine the statistically significant differences in gene expression 29 .
To obtain as much information as possible, GO 30 and KEGG 31-33 pathway analyses were performed to identify significant functions of the DEGs. The Fisher's exact test was used to determine the enrichment in categories.

Quantitative reverse transcription polymerase chain reaction (qRT-PCR).
To verify the accuracy of the RNA-Seq data, qRT-PCR analysis of ten randomly selected DEGs (including five upregulated and five downregulated genes) was performed using the same extracted total RNA as the RNA-Seq analysis in each of the ten samples. Briefly, one microgram of RNA was used for reverse transcription with the TURE script 1st Stand cDNA SYNTHESIS Kit (Aidlab, Beijing, China) according to the manufacturer's instructions, and genespecific primer pairs were designed with Beacon Designer 7 based on transcriptome-assembled data. The primer sequence pairs are shown in Table 2. PCR amplifications were carried out using the qTOWER Real-Time PCR Thermal Cycler (Analytik JenaAG, Jena, Germany). The reaction was proceeded under the following conditions: 3 min at 95 °C, 10 s at 95 °C, and 30 s at 58 °C before the plate reads were taken, 10 s at 95 °C, 39 cycles, and final Melt curve analysis (60-95 °C, + 1 °C/cycle, holding time 4 s). The relative quantification of gene expression was computed using the 2 −ΔΔCt (C T reference gene − C T target gene ) method with 16S rRNA as the reference gene 34 . All experiments were performed in triplicates to ensure accuracy. Statistical analysis. Statistical analysis was performed using Prism 6 (GraphPad Software, La Jolla, CA, USA). Differences between groups were assessed using the Student's t-test or Mann-Whitney U test. Validated www.nature.com/scientificreports/ data from qRT-PCR are expressed as the mean ± standard deviation (SD). A p value < 0.05 was considered statistically significant.

Results
Reads generation. Ten cDNA libraries prepared from the 10 strains of S. epidermidis were sequenced on an Illumina HiSeqTM 2000 platform. A total of 234,631,384 raw reads were generated, and after the sequencing adapters and low-complexity reads were removed, 213,509,502 clean reads were generated with an average of 21,059,504 reads in the P-Se cDNA libraries and 21,642,396 reads in the NP-Se cDNA libraries. The clean reads which were perfectly mapped to the reference genomes without mismatch were all greater than 70%. A detailed summary of the sequencing results is shown in Table 3.
Differentially expressed gene analysis. A total of 281 DEGs were identified between postoperative endophthalmitis isolates and healthy conjunctival isolates, of which 142 were upregulated and 139 were downregulated. A scatter plot was used to assess the expression variation of the genes between P-Se and NP-Se ( Fig. 2A). Moreover, a volcano plot was constructed to visualize the DEGs between two groups (Fig. 2B). The annotation of the 281 DEGs in the NCBI-nr database shows that the gene codes for Staphylococcal toxin (SE1634) and phenolsoluble modulins β1 (SE0846) were significantly upregulated in P-Se, which might be directly responsible for pathogenesis of S. epidermidis endophthalmitis. A complete list of all DEGs is shown in Supplementary Table S1.

Validation by qRT-PCR.
To validate the results of RNA-Seq, ten DEGs (including five upregulated and five downregulated genes) were randomly selected for qRT-PCR analysis on the basis of the fold-change and p values. The expression levels of SE1763, SE0727, SE1634, SE0124, and SE0125 were upregulated, while those of SE2185, SE2389, SE2184, SE0241, and SE2183 were downregulated. Specifically, the genes SE1763, SE0727, SE1634, and SE0125 were significantly upregulated in the group of P-Se compared to NP-Se ( Fig. 3A-D), while SE0124 (Fig. 3E) had no significant difference. Similarly, the genes SE2389, SE2184, and SE0241 were significantly downregulated in the group of P-Se compared with their levels in the NP-Se group (Fig. 3F-H), while SE2185 (Fig. 3I) and SE2183 (Fig. 3J) were of no significant difference. In addition, we compared the log2 fold change of ten selected DEGs between RNA-Seq and qRT-PCR (Fig. 4). The qRT-PCR data were consistent with the RNA sequencing data, indicating the reliability of the RNA-Seq results.  -TCC TAC GGG AGG CAG CAG T-3′  5′-GGA CTA CCA GGG TAT CTA ATC CTG TT-3′   SE1763  5′-TCA ATG GAG AGT AGC AGA TA-3′  5′-TCC GAA AGT AGA ACC AAT AC-3′   SE0727  5′-AGA ATT TTA TGA CGG TTT AAG AGG -3′ 5′-ATC GCT TAT GAT TGA TAA CTG TCT -3′   SE1634  5′-GGC AGC AGA TAT CAT TTC TAC-3′  5′-ATC CAT TTT ACT AAA TCA CCG-3′   SE0124  5′-AAC GCC CGA ATT ATT TAA TGTC-3′  5′-ACT TTG ATG TGC TGT ATA ACCA-3′   SE0125  5′-GTG AGC GAA AAA GAG TTA TTG-3′  5′-ACG TAT ACT TTG GTT CAC CGT-3′   SE2185  5′-CAA TGG GCA ACA ACA ACA A-3′  5′-TTC GCT TCA TCA TAC TCA GTC-3′   SE2389  5′-TCG GTC TAA TCA CAC   www.nature.com/scientificreports/ GO functional analysis of DEGs. GO annotation is an international standard system of gene function classification that has three main categories to describe biological process (BP), cellular component (CC), and molecular function (MF). To gain an insight into the biological roles of the most significantly up-or downregulated genes, GO enrichment analysis was performed using the Fisher's exact test with p value ≤ 0.05 as the threshold. In brief, all DEGs were enriched to 37 GO terms and classified into categories of BP with 29 GO terms, CC with 1 GO term, and MF with 7 GO terms. As shown in Fig. 5 and Supplementary Table S2, the top three enriched terms in BP were "homeostatic process (GO: 0042592)", "cellular homeostasis (GO: 0019725)", and "cell redox homeostasis (GO: 0045454)", with the number of DEGs being 9, 8, and 5, respectively. In the CC category, "cell (GO: 0005623)" was the only but highly enriched term, with 6 DEGs attached. In the category of MF, genes were significantly enriched with the terms of "recombinase activity (GO: 0000150)", "hydrolase activity, acting on glycosyl bonds (GO: 0016798)", and "alcohol dehydrogenase (NAD) activity (GO: 0004022)", and the number of DEGs being 5, 5, and 4, respectively. Moreover, by analyzing these involved genes we found several DEGs were  www.nature.com/scientificreports/ closely associated with the thioredoxin system, which plays a crucial role in defense against oxidative stress and provides an opportunity to kill bacteria 35 . Several terms were also strongly associated with the BP of iron ions. The terms of "ferric iron binding (GO: 0008199)", "cellular iron ion homeostasis (GO: 0006879)", and "iron ion homeostasis (GO: 0055072)" were significantly enriched. Iron acquisition strategy is required for pathogenic bacterial colonization and subsequent pathogenesis 36 . More information about the involved genes of GO enrichment is shown in Supplementary  Table S2.

KEGG pathway analysis of DEGs. KEGG pathway was used for a biological pathway-based analysis,
in which all genes were mapped to the reference canonical pathways of KEGG. The pathways of "degradation of aromatic compounds (Ko01220)", "naphthalene degradation (Ko00626)", "tyrosine metabolism (Ko00350)", "chloroalkane and chloroalkene degradation (Ko00625)", and "butanoate metabolism (Ko00650)" were significantly enriched (p ≤ 0.05) (Fig. 6B). The pathways of "two-component system (Ko02020)" and "pyruvate metabolism (Ko00620)" were remarkably enriched with more DEGs but without strong significance (p > 0.05) (Fig. 6A). These two pathways should not be ignored because of their importance in exploring the pathogenesis and antibiotic resistance of S. epidermidis 37,38 . More information about the results of KEGG enrichment is presented in Supplementary Table S3.

Discussion
Based on the analysis of 16S rDNA, many different types of microorganisms have been found on the ocular surface, mainly including harmless commensal organisms and opportunistic pathogens 39,40 . S. epidermidis, formerly considered commensal, has been recognized as an important opportunistic pathogen in nosocomial infection although it is commonly isolated from the skin and mucous membranes of healthy individuals 41,42 . The differences between commensal and nosocomial S. epidermidis isolates were reported at the genetic level. A genome sequencing study suggested that commensal skin S. epidermidis isolated from different parts of the body of healthy individuals showed a high amount of genetic variation 43 . As for ocular infection isolates, significant genetic variations were detected between 42 isolates of keratitis and endophthalmitis and 14 healthy conjunctival isolates by using fluorescence-amplified fragment length polymorphism 20 . Differences were also found in genotype and phenotype between healthy conjunctival isolates and ocular infectious isolates 19 . In this study, we used RNA-Seq for transcriptome analysis to explore the DEGs of S. epidermidis isolated from the healthy conjunctiva and postoperative endophthalmitis. In total, 142 significantly upregulated and 139 significantly downregulated genes were identified in the strains from endophthalmitis. After GO and KEGG functional analyses of these DEGs, we detected several critical genes that may serve key roles in the pathogenesis and therapeutic targets of endophthalmitis caused by S. epidermidis.
S. epidermidis is generally believed to produce no aggressive toxins. However, aggressive members of PSMs have been identified in S. epidermidis 16,17 . PSMs, a family of proinflammatory peptides produced by most staphylococci, have multiple functions in the production of proinflammatory cytokines, immune evasion, biofilm development, cytolytic capacity, and even the killing of competing microbes 44 . Consistently, our results showed that the gene code for phenol-soluble modulins β1 (SE0846) was significantly upregulated in the strains from www.nature.com/scientificreports/ postoperative endophthalmitis (log 2 FC≈1.49). Meanwhile, another gene SE1634, which codes for Staphylococcal toxin, was also found highly upregulated in the group of P-Se (log 2 FC≈1.76). Staphylococcal toxins possess broadly cytolytic properties and strong proinflammation, which leads to tissue degradation and cell death 15,45 . The expression level of SE1634 was validated by qRT-PCR. Thus, it is inferred that the two genes may be directly responsible for the postoperative endophthalmitis caused by S. epidermidis. At present, qRT-PCR is still considered the method of choice for validation of gene expression data obtained on high-throughput profiling platforms. In this study, ten DEGs were randomly selected to validate the results of RNA-Seq. However, there are minor differences between both experimental conditions. The RNA-Seq processing workflows, the annotation of the reference transcriptome, reads mapping, primer design, and reagents have impacts on the results 46 . Therefore, careful validation is warranted when expression profiles are evaluated with this specific gene set.
GO functional annotation of the transcripts in the current study revealed that both up-and down-regulated genes were related to the metabolic process. The results of enrichment analysis indicated that the terms of homeostatic process, cellular homeostasis, and cell redox homeostasis were significantly enriched in the BP. By analyzing these involved genes, we found that most of the DEGs presented an upregulated trend. In detail, SE2097 and SE0594 codes for thioredoxin and SE0838 code for thiol reductase thioredoxin participated in the thioredoxin system and were significantly upregulated in the group of P-Se (log 2 FC ≈ 1.59, 1.36, and 1.90, respectively). The www.nature.com/scientificreports/ thioredoxin system, composed of thioredoxin, thioredoxin reductase (TrxR), and nicotinamide adenine dinucleotide phosphate (NADPH), is widely distributed in natural organisms and serves in defense against oxidative stress 35,47 . Moreover, the obvious differences in structure and reaction mechanisms of TrxR between bacteria and mammals have made this system a novel antibiotic target 35 , for example, in Bacillus anthracis as a new drug target to several important human pathogens 48 and in S. aureus as a feasible target for antibacterial drug design 49 . Thus, the thioredoxin system may be speculated as an attractive antibiotic target for treatment of endophthalmitis induced by S. epidermidis. Furthermore, several genes involved in the iron ion metabolism need to be mentioned. Iron is an essential micronutrient for the growth and proliferation of all organisms, including pathogenic bacteria, and plays a pivotal role in colonization and subsequent pathogenesis 36,50 . S. aureus has been confirmed to evolve sophisticated strategies to obtain iron during infection, suggesting the iron acquisition system might be a viable target for therapeutic interventions 36 . In this study, the DEGs SE1764, SE1578, and SE1783, which were involved in the iron ion metabolism, were significantly upregulated in the strains of S. epidermidis isolated from postoperative endophthalmitis. Moreover, the results of KEGG pathway enrichment analysis indicated that pathways of the two-component system and pyruvate metabolism enriched more DEGs. The two-component system contains a membraneassociated histidine kinase and a response regulator to regulate bacterial adaptation, survival, virulence, and biofilm formation 51 . For opportunistic bacterial pathogens, the existence of essential two-component signal transduction systems is a core element of virulence and antibiotic resistance, suggesting that these systems may be targets for antimicrobial interventions 37 . In the present study, seven DEGs were enriched into this pathway, with SE1637 code for histidine kinase and SE1635 code for accessory protein regulator protein B both significantly upregulated in the strains from postoperative endophthalmitis. In addition, data from this study indicated that six DEGs were enriched into the pyruvate metabolism pathway. Pyruvate, a critical metabolite in a variety of anabolic and catabolic pathways including the oxidative metabolism, gluconeogenesis, and tricarboxylic acid cycle, is essential for basic activities of organisms 52 . Harper et al. demonstrated that pyruvate induced robust leucocidin production and enhanced the pathogenicity of S. aureus, so it may serve as a novel regulatory signal to coordinate S. aureus virulence through intricate regulatory networks 38 . It was also reported that pyruvate was a physiological ligand participating in regulating the network of sporulation in Bacillus subtilis 53 . However, the role of pyruvate metabolism in S. epidermidis is still not clear.
Our study has several limitations. First, we only collected five S. epidermidis isolates from post-cataract endophthalmitis because this type of endophthalmitis has an extremely low incidence, ranging from 0.033 to 0.36% 1-3 , and the positive rate of bacterial culture is also not high. Second, we did not analyze the bacterial phenotype of these isolates in details, which is important for a comprehensive understanding of microorganisms 16,54 . Third, we just explored the differences at the transcriptome level of S. epidermidis isolated from different pathological statuses, but did not carry out genome analysis. Fourth, lots of implications of genes are thought to contribute to S. epidermidis endophthalmitis, and the corresponding verification needs to be implemented in the future investigations.
In conclusion, the present study identified 281 DEGs of S. epidermidis isolates from postoperative endophthalmitis and the healthy conjunctiva using RNA-Seq. With GO and KEGG pathway analyses, several important DEGs and their potential biology pathways were revealed. The findings not only help to understand the pathogenesis of S. epidermidis, but also provide a basic resource for identification of therapeutic targets in postoperative endophthalmitis caused by S. epidermidis. www.nature.com/scientificreports/