Proteomic and Transcriptomic Analyses of Swine Pathogen Erysipelothrix rhusiopathiae Reveal Virulence Repertoire

E. rhusiopathiae is the causative agent of erysipelas in animals and erysipeloid in humans, but its pathogenicity is poorly understood. To identify virulence factors associated with E. rhusiopathiae and screen engineered vaccine candidates, we used proteomics and transcriptomics to compare the highly virulent strain HX130709 with an isogenic avirulent derivative, HX130709a. 1,299 proteins and 1,673 transcribed genes were identified and 1,292 of the proteins could be associated with genes. In a comparison between HX130907 and HX130709a, 168 proteins and 475 genes exhibited differences in regulation level. Among these, levels for 61 proteins and transcripts were positively or negatively correlated. Gene Ontology (GO) analysis suggests that many of the down-regulated proteins in the attenuated strain have catalytic or binding functions. Potential protein-protein interactions suggest that some of the down-regulated proteins may regulate PTS, GMP synthase and ribosomal proteins. Morphological results showed that HX130709 and HX130709a have similar colony and capsule morphology. Growth curves and pyruvate measurements suggest that TCA cycle and saccharide phosphorylation levels were decreased and gluconeogenesis was increased in HX130709a. Our study confirms that SpaA and neuraminidase, but not hyaluronidase and capsule, are associated with virulence in E. rhusiopathiae. We conclude that the virulence of E. rhusiopathiae may be associated with slow reactions of the TCA cycle and down-regulation of selected proteins.


Introduction
Erysipelothrix rhusiopathiae (E. rhusiopathiae) is a small gram-positive, slender, and straight rod-shaped bacterium that can cause erysipelas in swine and other animals, including sheep, fish, reptiles and birds. E. rhusiopathiae is also an important pathogen in humans and is the causative agent of erysipeloid, a skin disease [1]. The bacterium can be isolated from sick and healthy animals (pork, seafood, chicken) and from environment in which they live [2,3]. E. rhusiopathiae belong to the genus Erysipelothrix along with E. tonsillarum and two other tubes, flash frozen in liquid nitrogen, and submitted to BGI tech for proteomic and transcriptomic analyses. Three biological replicates were prepared independently for each sample.

Quantitative transcriptomics (RNA-seq)
RNA isolation and mRNA purification. Total RNA was isolated with TRIzol reagent using the standard protocol, and dissolved in 200 μL RNase-free water. The concentration of total RNA was determined using a NanoDrop spectrophotometer (Thermo Scientific, USA), and the RNA integrity value (RIN) was determined using the RNA 6000 Pico LabChip and an Agilent 2100 Bioanalyzer (Agilent, USA). Total RNA was incubated with 10 U DNase I at 37°C for 1 h, and then nuclease-free water was added to bring the sample volume to 250 μL. Messenger RNA was further purified by digesting ribosomal RNA and tRNA with Terminator 5΄phosphate-dependent exonuclease. The resulting RNA samples were quantified using a DU800 spectrophotometer (Beckman Coulter, USA) and mixed with fragmentation buffer to generate short mRNA fragments.
cDNA synthesis and Illumina sequencing. cDNA was synthesized using the mRNA fragments as templates. The short cDNA fragments were purified and resuspended in EB buffer for end repair and single nucleotide A (adenine) addition, then ligated to adapters. After agarose gel electrophoresis, appropriately sized products were selected as templates for PCR amplification. During the QC steps, the Agilent 2100 Bioanalyzer and ABI StepOnePlus Real-Time PCR System were used to monitor sample library quantity and quality. The library was sequenced using the Illumina HiSeq™ 2000 high-throughput sequencing system.
Bioinformatics analysis. Sequencing reads were mapped to the reference genome Fujisawa (NC_015601) using BLASTN with a threshold e value of 0.00001 and the "-F F" parameter [17], which allowed alignments with up to two mismatches. Reads that mapped to rRNA genes or that failed to align using these parameters were excluded from further analysis. The read totals were expressed as RPKM (Reads Per Kilo bases per Million reads) [18] for each gene, and then differently regulated genes were identified using the DEGseq package and the MARS (MA-plot-based method with Random Sampling model) method [19]. We used FDR 0.001 and an absolute value of log 2 Ratio 1 as thresholds to judge the significance of differences in gene expression.

Quantitative proteomics
Protein preparation. Harvested bacteria were washed three times with ice-cold phosphate-buffered saline (137 mM NaCl, 2.7 mM KCl, 10.1 mM Na 2 HPO 4 , 1.8mM KH 2 PO 4 , pH 7.4). The supernatant was discarded after the final centrifugation at 12,000 × g for 30 min. The pellets were extracted with lysis buffer (7 M Urea, 2 M Thiourea, 4% CHAPS, 40 mM Tris-HCl, pH 8.5) with a final concentration of 1mM PMSF and 2mM EDTA. After 5 min of vortexing, DTT was added to the samples to a final concentration of 10 mM. The samples were sonicated at 200 W for 15 min and then centrifuged at 4°C, 30,000 × g for 15 min. The samples were mixed well with 5× volume of ice-cold acetone containing 10% (v/v) TCA and incubated at -20°C overnight. After centrifugation at 4°C, 30,000 × g, the supernatant was discarded. The precipitate was washed with ice-cold acetone three times. The pellet was vacuum-dried and dissolved in lysis buffer (7M urea, 2 M thiourea, 4% NP40, 20mM Tris-HCl, pH 8.0-8.5). The samples were sonicated at 200 W for 15 min and centrifuged at 4°C, 30,000 × g for 15 min, and then the supernatant was transferred to a new tube. To reduce disulfide bonds in proteins, 10 mM DTT (final concentration) was added and samples were incubated at 56°C for 1 h. Subsequently, 55 mM IAM (final concentration) was added to block the cysteines, and samples were incubated for 1 h in the dark. The supernatant was mixed well with 5× volume of ice-cold acetone for 2 h at -20°C to precipitate proteins. After centrifugation at 4°C, 30,000 × g, the supernatant was discarded and the pellet was vacuum-dried for 5 min. The samples were then dissolved in 500 μL 0.5 M TEAB and sonicated at 200 W for 15 min. Finally, samples were pelleted at 4°C, 30,000 × g for 15 min. The supernatant was transferred to a new tube and protein content quantified. Protein preparations were stored at -80°C for later analysis.
iTRAQ (Isobaric tag for relative and absolute quantitation) labeling and SCX fractionation. 100 μg of total protein was withdrawn from each sample and then digested at 37°C for 16 hours with Trypsin Gold, at a protein: trypsin ratio of 30: 1. After digestion, peptides were vacuum-dried and re-dissolved in 0.5M TEAB, then processed with 8-plex iTRAQ reagent, according to the manufacturer's protocol. Briefly, one unit of iTRAQ reagent was thawed and mixed with 24 μL isopropanol. The peptides labeled with the isobaric tags were pooled then dried by vacuum centrifugation.
SCX chromatography was performed using a LC-20AB HPLC Pump system (Shimadzu, Kyoto, Japan). The peptide mixtures labeled with isobaric tags were reconstituted in 4 mL buffer A (25 mM NaH2PO4 in 25% ACN, pH 2.7) and loaded onto a 4.6×250 mm Ultremex SCX column containing 5-μm particles (Phenomenex). The peptides were eluted with a gradient of buffer A for 10 min, at a flow rate of 1 mL/min, then eluted with 5-60% buffer B (25mM NaH 2 PO4, 1 M KCl in 25% ACN, pH 2.7) for 27 min, and then with 60-100% buffer B for 1 min. Prior to the next sample injection, the system was maintained in 100% buffer B for 1 min, then equilibrated with buffer A for 10 min. Elution was monitored and fractions were pooled every 1 min. The eluted peptides were grouped into 20 fractions and desalted with a Strata X C18 column (Phenomenex) then vacuum-dried.
LC-ESI-MS/MS analysis based on Q EXACTIVE. Each fraction was re-dissolved in buffer A containing 2% ACN and 0.1% FA and centrifuged at 20000 × g for 10 min. The final concentration of peptide was about 0.5 μg/μL. 10 μL supernatant was loaded by autosampler onto a 2 cm C18 trap column for analysis in a LC-20AD nanoHPLC (Shimadzu, Kyoto, Japan). The peptides were eluted onto an analytical C18 column (inner diameter 75 μm) packed inhouse. The samples were loaded at 8 μL/min for 5 min, then the column was run with a gradient from 2% to 35% in buffer B (98% ACN, 0.1% FA) at 300 nL/min for 35 min. A linear gradient to 80% was run through the column for 2 min, 80% buffer B was maintained for 4 min then returned to 5% for 1 min.
The peptides subjected to nanoelectrospray ionization were analyzed by tandem mass spectrometry (MS/MS) in a Q EXACTIVE instrument (Thermo Fisher Scientific, San Jose, CA) coupled online with the HPLC. Orbitrap was used to detect intact peptides at a resolution of 70,000. MS/MS was used to select peptides in high-energy collision dissociation (HCD) operating mode with a normalized collision energy setting of 27.0 and Orbitrap was set at a resolution of 17,500. The 15 most abundant precursor ions above a threshold ion count of 20,000 in the MS scan with a subsequent Dynamic Exclusion duration of 15 s were analyzed with a datadependent procedure that alternated between one MS scan followed by 15 MS/MS scans. The electrospray voltage was set as 1.6 kV. Automatic gain control (AGC) was applied to optimize the spectra generated by the Orbitrap. The AGC target for full MS was 3e6 and 1e5 for MS2. For MS scans, the m/z scan range was between 350 and 2,000 Da. For MS2 scans, the m/z scan range was between 100 and 1,800.
Data analysis. Raw data files obtained from the Orbitrap were converted into MGF format using Proteome Discoverer 1.2 (PD 1.2, Thermo), [5600 msconverter]. Proteins were identified using the Mascot search engine (Matrix Science, London, UK; version 2.3.02) against a database (downloaded from NCBI http://www.ncbi.nlm.nih.gov/protein/?term=txid1648[Organism:exp]) containing 3,409 sequences. For protein identification, a mass tolerance of 20 Da (ppm) was used for intact peptide masses and 0.05 Da for fragmented ions, and one missing trypsin cleavage was permitted in the trypsin digests. Gln->pyro-Glu (N-term Q), Oxidation (M), and Deamidated (NQ) were specified as potential modifications, and Carbamidomethyl (C), iTRAQ8plex (N-term), and iTRAQ8plex (K) were set as invariable modifications. Peptide charge states were set to +2 and +3. We included an automatic decoy database search as part of the analysis. In this option, whenever a protein sequence from the target database is tested, a random decoy sequence of the same length is also tested. To reduce false positives, proteins required support from at least one unique peptide, and only peptides that exceeded "identity" at 99% confidence were accepted.
For protein quantitation, at least two unique peptides were required. The median ratio by Mascot was used to weight and normalize quantitative protein ratios. Proteins with p-values < 0.05 and fold changes >1.5 were considered to be significantly different between samples.
Functional annotation. The functional annotation of proteins was conducted using the Blast2GO application and the non-redundant (NR) protein database at NCBI. The KEGG database (http://www.genome.jp/kegg/) and the COG protein database (http://www.ncbi.nlm. nih.gov/COG/) were used to classify and group the identified proteins.

Growth curves
A single colony of HX130709 (F0) and HX130709a (F70) was used to inoculate overnight cultures in BHI medium. 300 μL overnight culture was inoculated into 30 mL BHI medium and incubated at 200 rpm. 500 μL samples were collected in triplicate after 1 h and centrifuged at 12,000 rpm for 3 min. The supernatant was discarded, pelleted bacteria were re-suspended in an equal volume of PBS, and optical density (OD600 nm) was measured.

Morphological observations
HX130709 and the attenuated derivative strain HX130709a were streaked onto BHI agar containing 0.1% Tween 80, incubated overnight, and colony morphology observed. Bacteria cultured to log phase were pelleted and fixed with 2.5% glutaraldehyde and then sectioned. Thin sections were stained with uranyl acetate and lead citrate and examined with a Hitachi H-7500 transmission electron microscope.

Lactate dehydrogenase activity
Bacteria were cultured to mid-logarithmic phase in BHI medium at 37°C. Triplicate samples were collected for each of 3 time points for the virulent strain HX130709 (5h, 6h, 7h) and the attenuated strain HX13079a (8h, 9h, 10h). Each 10 mL sample was centrifuged at 6,000 rpm for 10 min, then washed in PBS 3 times. Pelleted bacteria were resuspended in 500 μL PBS and lysed by sonication. The lysate was centrifuged at 12,000 rpm for 1 min and supernatant collected. 50 μL supernatant was added to wells in 96-well plates, lactate dehydrogenase (LDH) substrate was added, and the plates were incubated at 37°C for 30 min, following the protocol accompanying the CytoTox 96 Non-Radioactive Cytotoxicity Assay kit. Color development was halted by the addition of stop solution and optical density at 490 nm (OD490) was recorded using an ELISA plate reader (BioTek, USA). In parallel, sampled bacteria were plated at a final dilution of 10 −7 on BHI agar plates to determine bacterial concentration (colony forming units, CFU). Results were expressed as LDH/cfu.

Results and Discussion
Summary of RNA-seq and iTRAQ data Transcriptome data were obtained by high-throughput RNA-sequencing as described in materials and methods. Three RNA samples were analyzed from HX130709 and its attenuated derivative HX130709a (six samples total). After excluding low-scoring sequence reads, the average read length was 90 bp and the total reads obtained were 3,675,074, 4,692,139 and 3,839,099 for samples HX130709-1, -2, -3 respectively, and 3,251,469, 4,755,689 and 5,151,752 for samples HX130709a-1, -2, -3, respectively. Correlation analysis among the replicates of HX130709, and among the replicates of HX130709a, generated Pearson's correlation coefficients over 0.95 except in two cases. Because the coefficients for HX130709a-1 vs. HX130709a-2 (0.206) and HX130709a-1 vs. HX130709a-3 (0.2008) were unacceptably low, the HX130709a-1 data set was discarded. In spite of this loss, the mean sequencing coverage for both strains was 193-259 fold. Reads were mapped to the E. rhusiopathiae genome sequence (Fujisawa) using Soap (v2.01). The results indicated that the transcriptional percentages for ORFs encoded by the Fujisawa strain were 96.37% and 84.58% for HX130709 and HX130709a, respectively. Surprisingly, no SNPs were found in HX130709 or HX130709a. Gene differential regulation was analyzed using the DEGseq software package. Out of 1,673 genes examined, 197 genes exhibited increased transcription in HX130709a and 278 genes exhibited decreased transcription, relative to the levels in HX130709 (S1 Table).
Three biological replicates for each strain were included in the iTRAQ experiment comparing virulent HX130709 and avirulent HX130709a. After trypsinization and labeling with distinct isobaric tags, protein fragments in the six samples were separated and identified by LC-ESI-MS/ MS. Because the error distribution analysis among the replicates of HX130709, and among the replicates of HX130709a, showed that the mean errors for HX130709a-1 with replicas -2 and -3 were 1.37 and 1.43, respectively, the data obtained from HX130709a-1 was again set aside. In the five remaining samples, a total of 312,350 mass spectra were generated. After excluding lowscoring spectra, 79,681 unique spectra that matched to specific peptides were obtained. Mascot analysis identified 12,154 peptides, 11,655 unique peptides, and 1,299 proteins in total, in the five samples (S1 File). Covariance distribution analysis for these samples yielded a mean covariance of 0.063 (S1 Fig), indicating good reproducibility between biological replicates in each group. A protein with 1.5-fold difference and p-value 0.05 was regarded as being differentially regulated. Relative to protein levels in the virulent strain HX130709, the abundance of 101 proteins increased, and 67 proteins decreased, in HX130709a (Table 1).
Functional annotation of the 1299 proteins was based on the three principal classifications developed by the gene ontology (http://www.geneontol.ogy.org). With respect to biological processes, most proteins were associated with metabolic processes (32.18%) and cellular processes (26.68%). When classified according to cellular component, 30.76% of proteins were described as located in the cell, 30.76% in cell parts, 13.14% in membranes, and 9.45% in membrane parts. Classification by molecular function showed that most proteins had catalytic (46.97%) and binding (36.54%) activities (Fig 1). Overall, the gene ontology classifications for the differentially expressed proteins are very similar to those obtained for all 1,299 proteins identified by Mascot. However, when the differentially expressed proteins are classified by molecular function, catalytic activity (50.75% vs. 44.07%) and binding (34.33% vs. 29.66%) are overrepresented.
The COG database (http://www.ncbi.nlm.nih.gov/COG) was used to classify and group the identified proteins (Fig 2). The ten most common categories were as follows:

Correlation between mRNA and protein expression profiles
A coupled transcriptomics-proteomics project provides a unique opportunity to investigate whether protein regulation is correlated with gene transcription. We investigated the correlation between mRNA and protein profiles and found that the levels of 1219 proteins could be correlated, either negatively or positively, with mRNA levels. Among the 168 differentially regulated proteins, only 61 could be correlated with mRNA level variations. These proteins could be clustered into four groups based on the pattern of changes in mRNA and protein levels: Group I, the mRNA and protein levels are positively correlated (45 proteins); Group II, the mRNA level remains almost unchanged while the protein level is decreased (43 proteins); Group III, the mRNA level remains almost unchanged but the protein level is increased (64 proteins); Group IV, the mRNA level is decreased but the protein level is increased (16 proteins) ( Table 1). Group I includes two subgroups: both the mRNA and protein levels are increased synchronously (21 proteins); both the mRNA and protein levels are decreased synchronously (24 proteins) ( Table 1).

Protein-protein interaction analysis
Protein-protein interactions play an important role in bacterial pathogenicity and metabolism. We therefore examined the 67 down-regulated proteins for potential protein interactions using the STRING database [20]. Associations were predicted to exist among these proteins, and especially between them and components of the phosphotransferase system, GMP synthases and ribosomal proteins (Fig 3). These interactions may function in reducing nucleotide and protein synthesis and saccharide phosphorylation.

Phosphotransferase system (PTS)
The phosphoenolpyruvate sugar phosphotransferase system (PTS) is widespread among microorganisms. The PTS couples carbon source uptake and substrate phosphorylation, generating intracellular sugar-phosphate [21]. Pathway analysis suggests that the phosphorylation levels of glucose, maltose, D-glucosamine, N-acetylgalatosamine and galactosamine were down-regulated (1.6-5 fold) in HX130709a, resulting in reduced glucose 6-phosphate and glyceraldehydes-3-phophate synthesis, consistent with the up-regulation of gluconeogenesis (Fig 4).

TCA Cycle
A pathway analysis focused on the TCA cycle showed that the pathway components responsible for the conversion of pyruvate and citrate to oxaloacetate were up-regulated in HX130709a relative to HX130709. In contrast, components involved in the metabolism of succinate to fumarate were down-regulated. We hypothesize that the lower levels of medium-supplied carbohydrates were used by the attenuated E. rhusiopathiae, but the cells compensated by up-regulating production of oxaloacetate and shunting excess oxaloacetate into the gluconeogenesis pathway. However, pathway components involved in the conversion of pyruvate to acetyl-CoA were not changed. Combining this fact with the down-regulation of components involved in the metabolism of succinate to fumarate, we suggest that the TCA cycle as a whole is down-regulated in HX130709a (Fig 4). Growth curves also showed that HX130709 grows more rapidly than HX130709a (Fig 5). Compared with HX130709, pyruvate levels in HX130709a may be slightly higher, but the difference is not significant (Fig 6). The colony morphology of HX130709a was convex and irregular, similar to that of HX130709. Electron microscopy  showed capsule material as an electron-dense layer outside the outer membrane and no significant morphological difference between the virulent and attenuated strains (Fig 7), which is inconsistent with the result reported by Yoshihiro et al. [13,22]. The role of the capsule in virulence has been clearly demonstrated using isogenic mutants with defined mutations. However, differences in growth curves and pyruvate metabolism between parent strains and their avirulent derivatives have not been reported, nor have differences in mRNA and protein levels. Our results demonstrate that capsule is not closely associated with virulence in E. rhusiopathiae.

Virulence factors
In a previous study, Kwok and Li sequenced the genome of a virulent strain of E. rhusiopathiae (SY1027) and used BLASTP and VFDB to identify 37 potential virulence factors [23]. Using a similar strategy, we examined the proteins and genes identified by proteomic and transcriptomic analyses in our study and found 13 virulence factor candidates. To our surprise, compared with HX130709, no potential virulence factors were differentially expressed in HX130709a. However, it is possible that these factors are not associated with pathogenicity in E. rhusiopathiae. Virulence factors such as neuraminidase, hyaluronidase, RspA/RspB, SpaA, and hemolysin have been suggested to play a role in E. rhusiopathiae pathogenicity. Among the proteins and genes identified in our transcriptomic and proteomic analyses, transcriptional activity was detected for neuraminidase, hyaluronidase and spaA, which were decreased 2.7, 2.35 and 15.8 fold in HX130709a vs. HX130709, respectively. Unexpectedly, neuraminidase was not detected as a protein product, nor were transcripts detected that corresponded to the rspA/rspB genes. Since pelleted bacteria were used in the proteomic and transcriptomic analyses, it is possible that neuraminidase was secreted into the culture medium and that the mRNAs for rspA/rspB were unstable. Protein levels for hyaluronidase and RspA were stable, but protein levels for RspB and SpaA were increased and decreased 3.364 and 4 fold, respectively. Hemolysin was not detected in mRNA or protein forms, but a hemolysin-like protein was identified and its protein level was down-regulated.
Neuraminidase is an enzyme responsible for cleavage of sialic acids from sialo-glycoconjugates such as glycoproteins, glycolipids, and oligo-and polysaccharides. The down-regulation of neuraminidase has been associated with decreased virulence in E. rhusiopathiae [10]. We found that the mRNA level of neuraminidase was decreased 2.7 fold in attenuated strain HX130907a, which is consistent with this result. Shimoji et al. conducted transposon mutagenesis with Tn916 to construct mutants defective in hyaluronidase production and reported that hyaluronidase was not associated with virulence in E. rhusiopathiae [24]. Although mRNA levels for hyaluronidase decreased 2.35 fold in our study, protein expression were essentially identical in HX130709 and HX130907a, providing additional support that hyaluronidase is not associated with virulence in E. rhusiopathiae. RspA expression was also identical in the two strains, but RspB was increased in the attenuated strain HX130907a, suggesting that RspA and RspB are not associated with virulence in E. rhusiopathiae. Finally, Borrathybay et al. deleted the spaA gene from wild-type E. rhusiopathiae strain C43065 and found that the virulence of the ΔspaA mutant decreased more than 76 fold [25], but they did not compare growth curves for the wild-type and mutant strains. We also found that spaA transcription and SpaA protein levels decreased in HX130907a, supporting the hypothesis that SpaA is associated with virulence in E. rhusiopathiae.
In summary, the transcriptomes and proteomes of an attenuated E. rhusiopathiae and its parent strain were compared. 475 genes and 168 proteins were found to be up-or down-regulated and the levels for 61 proteins could be correlated with gene transcription levels. The growth of the attenuated strain is slower than its parent strain, but pyruvate metabolism appears to be unaffected. Our data are consistent with other studies showing that SpaA and neuraminidase, but not hyaluronidase and capsule, are associated with virulence in E. rhusiopathiae. We conclude that the down-regulation of the TCA cycle and the down-regulation of several proteins are associated with virulence in this organism.