Deep Sequencing Analysis of the Eha-Regulated Transcriptome of Edwardsiella tarda Following Acidification

Eha is a virulence regulator associated with the replication of Edwardsiella tarda within RAW264.7 macrophages. Eha is required for the bacteria to resist the acid and oxidative stress in macrophages. We herein demonstrate that Eha regulates the resistance of this bacterium against acidification in macrophages and explain the underlying molecular mechanism. Firstly, to find an acid or oxidative condition to induce the strongest Eha activities, we constructed pMP220Pehalac Z plasmid and inspected the lac Z expression regulated by Eha by using a β-galactosidase assay. At exposure of pH6.3 medium 2 h, whole transcriptomic profiling of the wild type and mutant were performed by RNA-sequencing. We identified 147 differentially -expressed genes (DEGs) (|log2 ratio| ≥ 1), 113 and 34 of which were significantly upand down-regulated, respectively, in the mutant compared with the wild type. These findings were validated by qRT-PCR. A CO functional analysis revealed that these genes were divided into 25 categories, including the bacterial processing, localization, metabolism, combination, catalysis, transportation and cellular composition. Based on the KEGG database, these genes were distributed in 55 pathways, such as the two-component system, ABC transporters, and microbial metabolism. At last, the intracellular survival rates and intraphagosomal pH of wild type ET 13 and its eha mutant in bafilomycin-treated and untreated macrophages were measured. The experiment showed that Eha was involved in protecting the bacteria from the effects of acidification within macrophages. The survival rate of the wild was also higher than that of the mutant under acid stress both in vivo and in vitro (P<0.05). Overall, Eha was found to regulate 147 genes that affect bacterial metabolism and virulence, allowing the bacteria to adapt to an acidic environment. These results could be helpful for further investigations of the mechanisms by which Edwardsiella tarda survives in macrophages. M et ab olo mics: Openccess ISSN: 2153-0769 Metabolomics: Open Access Citation: Gao D, Liu N, Li Y, Zhang Y, Liu G, et al. (2017) Deep Sequencing Analysis of the Eha-Regulated Transcriptome of Edwardsiella tarda Following Acidification. Metabolomics (Los Angel) 7: 196. doi:10.4172/2153-0769.1000196


Introduction
Edwardsiella tarda (E. tarda) is a facultative intracellular pathogen that causes Edwardsiellosis in freshwater and marine fish worldwide and a Salmonella-like gastroenteritis in humans [1,2]. Macrophages play critical roles in the defense against invading bacteria. Zhang et al. showed that E. tarda can utilize macrophages as a niche to initiate its virulence and spread systemic infections [3]. However, there are various stressful conditions such as nutrition deprivation, low pH and high reactive oxygen species (ROS) present in the intracellular niche of the phagocytes [4]. There are the strategies that virulent bacterial strains use to evade phagosomes and escape into the cytosol to enable their survival in the cells [5]. E. tarda could live and multiply in macrophages, and escape from the cells. It is important for the bacterium to lead to extra intestinal diseases and systemic infections [6]. E. tarda is capable to detoxify ROS by generating catalase (Kat) and superoxide dismutase (Sod) to survive within macrophages [7,8]. Okuda et al. suggested that the bacterial type III secretion system is able to interfere with the formation of acid stress in phagosomes to facilitate its replication in macrophages [9]. However, little other information is available about how E. tarda is affected by an acidic environment.
Some findings have been reported about the changes that occur in other intracellular pathogens in response to acid stress in macrophages. For example, Rathman et al. suggested that an acidic environment in Salmonella-containing phagosomes is necessary for the bacterial survival and replication within the macrophage [10]. Supporting this finding, it has been shown that some virulence genes in Salmonella typhimurium are activated within acidified phagosomes [11]. It has also been reported that the VirB secretion apparatus in Brucella is a kind of type IV secretion system. The acidification arising from phagosomes is a key intracellular signal that induces virB expression [12].

Bacterial strains, plasmids, growth conditions and cell culture
The ET13 wild type strain of E. tarda was provided by Dr. Janda [18] (California Department of Health Services Microbial Diseases Laboratory). The eha mutant strain and the complemented strain, ehaComp were constructed previously [14]. The pMP220 plasmid [19] was kindly provided by Prof. Mao (Department of Biochemistry, Southeast University). The eha gene was deleted from ET13 using the pHM5 suicide plasmid with homologous recombination to get the eha mutant strain. The complemented strain, ehaComp (the eha mutant containing plasmid pACYC184-eha), was constructed previously [14]. The pMP220-lac Z plasmid was kindly provided by Prof. Mao (Department of Biochemistry, Southeast University). Bacteria were grown in Luria-Bertani (LB) broth (Sunshine Biotechnology, Nanjing, China).

Bacterial survival and replication in bafilomycin-treated macrophages
The assay was performed according to the method reported by Gao [16]. Briefly, RAW264.7 macrophages were seeded at approximately 5.0 X 10 5 cells per well in 24-well tissue culture plates and incubated overnight. The cells in experimental groups were pretreated with 1 μM bafilomycin (BAF) A1 for 30 min, and the cells in the control groups were cultured without the inhibitor. Then the cells were infected for 1 h with the mid-log wild-type or mutant bacteria at a multiplicity of infection (MOI) of 10:1. The cells were incubated in DMEM medium containing gentamicin (100 μg/ml) for 1 h to kill the exocellular bacteria. After designated time intervals (2, 4, 6 h) post-infection, the macrophages were lysed with 1% Triton X-100 for 10 min, then the solutes were diluted, spread onto LB agar plates, and the plates were incubated at 37°C overnight. The numbers of the colony-forming units (cfu) per milliliter of viable intracellular bacteria on the plates were counted.

Measurement of the intraphagosomal pH in E. tardacontaining macrophages
The phagosome acidification assay was performed according to the manufacturer's instructions (SunShineBio, Nanjing, China) and was modified as described previously [20]. First, RAW264.7 cells were incubated with the mid-log wild-type or mutant bacteria at a multiplicity of infection (MOI) of 10:1. These were then labeled with FITC (SunShineBio, Nanjin, China) (pH-sensitive) and iFluorTM647 (AAT Bioquest, USA) (pH-insensitive) fluorescent dyes for 2 h. The cells were washed twice and permeabilized in ice-cold Dulbecco's phosphate-buffered saline at a fixed pH (ranging from 3 to 8) in buffers comprising (DPBS)/ ethylenediaminetetraacetic acid containing 0.05% Triton X-100 in order to calculate the pH of E. tarda-containing vacuoles. The samples were immediately analyzed by flow cytometry to determine the ratio of the mean fluorescence intensity (MFI) of the emission between FITC and iFluorTM Fluor 647. A standard curve was made on the basis of the MFI for different pH in the phagosomes. Second, the cells in experimental groups were pretreated with 1 μM bafilomycin (BAF) A1 for 30 min, while the cells in the control groups were cultured without the inhibitor. After the cells were infected by the bacteria and labeled with the fluorescent dyes, they were washed twice in DPBS buffer without Triton X-100. The values were compared with the standard curve generated from the cells that had phagocytosed E. tarda for different time periods (2 h, 4 h, 6 h) before the flow cytometry analysis.

Recombinant plasmid construction and bacterial electrotransformation
The forward primer: CAGATCTAATGGGTGAACCACGCAAAT and reverse primer: CTCTAGAATCCATAACCCCGATACACC were used to amplify a 358bp fragment from ET13 modified using BglІІ and XbaІ. To construct the pMPE (pMP220-PehalacZ) plasmid, the 358bp eha promoter was inserted into pMP220-LacZ at BglII and XbaI sites. Inserted fragments from the pMPE plasmid were sent to the Sunshine Biotechnology Company for sequencing (Sunshine Biotechnology, Nanjing, China).

Bacterial growth, RNA isolation and mRNA enrichment
After the wild type ET13 and its eha mutant were cultured in LB (pH=7.2) to an OD600 of 0.8 at 37°C, these bacteria were treated with LB (pH=6.3) for 2 h. Total RNA was isolated from the bacteria using a RNA extraction kit (TaKaRa company, Dalian, China). The purity, integrity and RIN (RNA Integrity Number) of the RNA were assessed using an Agilent Bioanalzyer. The bacterial mRNA was enriched by the removal of 16S and 23S rRNA from the total RNA by poly-Toligoattached magnetic beads using a Microb Express Bacterial mRNA Enrichment Kit (Ambion, Shanghai, China).

cDNA library construction and sequencing using the Illumina genome analyzer
Complementary DNA (cDNA) libraries were built as described previously [23]. Briefly, the mRNA in each sample was fragmented into short sequences with divalent cations by heating in a proprietary buffer. Using the short fragments as templates, the first-strand cDNA was synthesized with random hexamer primers and reverse transcriptase (Invitrogen, Shanghai, China). The second-strand DNA was synthesized using RNase H and DNA polymerase I (Invitrogen, Shanghai, China). The amplified fragments were purified with a Qia Quick PCR Purification kit (Qiagen, Hilden, Germany), and poly (A) was added using Illumina PCR Primer Cocktail in a 35-cycle PCR. Two samples were sent to Shenzhen Hua Da Gene Company and sequenced using the Illumina sequencing 2000 platform. All sequences were Metabolomics, an open access journal ISSN: 2153-0769 examined for possible sequencing errors. Clean reads were obtained by removing the raw reads containing adapter, poly-N, and any lowquality results. The RNA-seq data have been deposited in the National Center for Biotechnology Information-Sequence Read Achieve under accession number SRX1898774.

RNA-seq alignment and identification of the eha-dependent genes altered by acid stress
The clean reads were mapped onto the genome sequences of E. tarda ATCC15947 (http://www.ncbi.nlm.nih.gov/nuccore/) to reflect the distribution and coverage of the reads in the reference genome using the Short Oligonucleotide Analysis Package (SOAP). The gene expression of the two samples (wild type ET13 and its eha mutant) was quantified using the Reads Per Kilobase of coding sequence per Million reads (RPKM) algorithm. The genes were considered to be differentially expressed if the difference in RPKM values between the two samples was ≥ 2.0-fold and the p-value was <0.005. The p-value was adjusted by the false discovery rate (FDR).
The Gene Ontology (GO) is an international standardized gene function classification system that is widely used to classify orthologous proteins. GO annotations for all possible eha-dependent genes were obtained from the GO database (http://www.geneontology.org/). The Kyoto Encyclopedia of Genes and Genomes (KEGG) is the major public pathway-related database that is used to identify significantly enriched metabolic pathways or signal transduction pathways (http:// www.genome.jp/kegg/).

qRT-PCR
The wild type ET13 and its eha mutant were cultured in the condition just like preparing RNA-seq. Total RNA was isolated from the bacteria using a RNA extraction kit (TaKaRa company, Dalian, China). The primers shown in Table 1 were used for quantitative realtime PCR (qRT-PCR). Using a SYBR ExScript qRT-PCR kit (Takara, Dalian, China), qRT-PCR was performed in an ABI 7300 real-time detection system (Applied Biosystems Company, USA). Bacterial 16S rRNA was used as a control. The 2 -△△CT method [24] was used to compare the mRNA levels in wild type and mutant bacteria.

Statistical analysis
Probability (P) values were determined by Student's t-test, and analyses were performed using the GraphPad Prism (version 6) software. The P values were considered to be significant when they were <0.05. The data for qRT-PCR were analyzed statistically by using twotailed t-tests. Pearson's correlation was used to evaluate the relationship between the results of qRT-PCR and RNA-seq, and was determined using SPSS. The difference was considered to be strong and significant when the correlation coefficient was >0.75; a moderate correlation was considered to be present when the coefficient was<0.75 and>0.45.

Identification of the conditions leading to the highest Eha activities in EHA MPE transformants
To determine under what conditions the EHAMPE expresses the strongest Eha activities, we examined the β-galactosidase activity (regulated by Eha) of these EHAMPE transformants. Our results showed that the highest β-galactosidase activity was observed after 2 h exposure, when the transformants treated with LB (pH=6.0) for different time ( Figure 1A). Our results also found that the highest β-galactosidase activity was present at pH=6.3 after 2 h exposure ( Figure 1B), when the transformants treated with different pH of LB. At the same time, our results showed that the highest β-galactosidase activity was observed after for 2 h exposure, when the transformants treated with 0.5% H 2 O 2 of LB for different time ( Figure 1C). Our results also found that the highest β-galactosidase activity was present to 0.1% H 2 O 2 of LB after 2 h exposure, when the transformants treated with different concentration of H 2 O 2 of LB ( Figure 1D). Finally, we decided to select the condition was LB (pH=6.3) for 2 h ( Figure 1E).

Evaluation of the transcriptomes of the wild type and eta mutant bacteria
The RNA concentrations isolated from wild type ET13 and its eha mutant were 1.905 and 1.255 μg/μl, respectively. Their purities were OD260/OD280 ≥ 1.9 and OD260/OD230 ≥ 2.0, RIN:10.0. The two kinds of RNA were considered to be of sufficient quality to construct a cDNA library (Table 1S).
Illumina sequencing of the wild type and mutant bacteria cultured under acid stress generated raw reads per library, respectively. After filtration to remove the low-quality reads, a total of 12,233,396 (wild type) and 12,122,800 (mutant) clean reads were separately obtained from the raw reads of the two samples. More than 94% of the clean reads of 11501727 (94.02%) for wild type and 11443179 (94.39%) for mutant completely matched the reference genome of E. tarda ATCC15947. Among the clean reads, there were 7233093 (59.13%; wild type) and 6313604 (51.48%; mutant) perfect matches, and 4268634 (34.89%; wild type) and (34.91%; mutant) had a ≤ 5 bp mismatch (Table 2S).

Gene Name
Primer Sequence (5'→3') Product (bp) It is important to annotate sequences, because it can reveal the molecular mechanisms underlying changes in gene expression. Using the SOAP software to assess the randomness and coverage of the clean reads mapped to the reference genome, our results showed that different sizes of clean reads were well-distributed in the reference genome ( Figure 1S) and the degree of coverage for more than 90% of the clean reads exceeded 92% of the reference genome ( Figure 2S).

Functional classification of the eha-dependent genes
We quantified the transcriptional activity of the genes using the RPKM algorithm. By comparing the expression levels of the eha mutant with the wild type bacteria, we found that 147 genes were significantly differentially transcribed between the two samples under acid stress (|log2 ratio|≥ 1, FDR≤0.001, and q value < 0.005). Among these genes, 113 were up-regulated and 34 were down-regulated in the eha mutant, compared with the wild type.
Functional assignments were defined by Gene GO terms, which provided a broad functional classification of the genes and their products. A GO category analysis revealed that 126 among the 147 genes were classified into 25 functional categories. These genes were mainly involved in binding (59 genes), catalytic activity (83 genes), cell components (37 genes), cell membranes (33 genes), cellular processes (68 genes), localization (30 genes), metabolism (78 genes), and transport (27 genes) (Figure 2). These were considered to be genes that were transcribed in an eha-dependent manner in response to acid stress.

Pathway analysis of the eha-dependent genes
To further investigate the biological functions and interactions of these genes, a pathway-based analysis was conducted using the KEGG pathway database, which contains all known networks of molecular interactions in different species. The KEGG analysis showed that 130 of the 147 total DEGs were associated with 55 pathways, including amino acid metabolism (10 pathways); carbohydrate, lipid and energy metabolism (16 pathways); nucleotide metabolism (four pathways); sulfur metabolism (three pathways); nitrogen metabolism; iron transport and so on. The pathway analysis may help to understand the interactions among these genes. The metabolic pathways of the DEGs were involved in two-component systems (20 genes), ABC transporters (12 genes), microbial metabolism in diverse environments (21 genes), biosynthesis of secondary metabolites (17 genes) and flagellar assembly (3 genes) ( Table 2).

Validation of RNA-seq analysis by qRT-PCR
It is common practice to validate the DEGs found by the RNAseq using another method, such as qRT-PCR. We therefore compared the qRT-PCR expression levels of 15 genes randomly selected from the wild type and mutant bacteria. The mean values of the log2 ratio of the eha mutant to the wild type determined by qRT-PCR were in good agreement with the log2 ratios of the eha mutant to the wild type from the RNA-seq analysis for 14 of the 15 genes (Figure 3), with a significant correlation noted for these genes (R2=0.799>0.75). There was only one exception (ETATCC -RS11570), therefore, the data from the qRT-PCR analysis validated the findings from the RNA-seq analysis.

Eha expression provides E. tarda with a survival advantage
Our previous studies have shown that ET13 can survive within macrophages for extended periods of time [16]. It has been indicated that the acidification of phagosomes containing Yersinia pseudotuberculosis in macrophages can be blocked by BAF A1, a vacuolar proton-ATPase pump inhibitor [25]. To better characterize the persistence of ET13 against acidification in macrophages and to follow the fate of the internalized bacteria after BAF A1 was used to block acidification, we first accurately monitored the survival rate of ET13 in the cells beginning the second hour post-infection by cell lysis and plate counting. As controls, other cells were cultured without the inhibitor. The pretreatment had Figure 1S: The randomness assessment of the clean reads of the eha mutant and its wild type ET13 mapped to the reference genome. The results showed that there were different sizes for the clean reads of the eha mutant (A) and the wild type ET13 (B) which were distributed in the reference genome of E. tarda ATCC15947.  There were 147 DEGs in the eha mutant compared with the wild type bacteria (|log2 ratio|≥1). The functional categories of these genes shown were up-regulated (red color) and down-regulated (green color) by Eha in ET13 bacteria subjected to acid stress.

Pathway
Gene ID Gene description RNA-seq* no observable deleterious effects at the early time points (2 h). However, after four or six hours, the numbers of wild type and mutant bacteria recovered from BAF-treated macrophages were higher than those for the untreated cells (P<0.05) ( Figure 4A). To more precisely determine the pH of the lumen of E. tarda-containing vacuoles, we employed a ratiometric assay in which Alexa Fluor 647/FITC-labelled E. tarda was used as a pH-sensitive probe. The pH in vacuoles containing the wild type or mutant recovered from BAF-treated macrophages was higher than that for the untreated cells at 2 h, 4 h and 6 h post-infection (P<0.05) ( Figure 4B). These results showed that acidification is important for inhibiting the survival and replication of E. tarda in macrophages.
Second, we wanted to better explain how Eha regulates the resistance of E. tarda against acidification in macrophages and to compare the fates of the internalized wild type and mutant bacteria. The CFU of the wild type bacteria in untreated macrophages was higher than that of the mutant at 4 h and 6 h post-infection (P<0.05) ( Figure 4A). Moreover, at the later time points (4 h, 6 h), there were no differences in the pH of vacuoles containing the wild type and mutant bacteria in untreated macrophages (P＞0.05) ( Figure 4B). Therefore, Eha may be involved in protecting the bacteria against the effects of acidification.

The eha mutants exhibited decreased survival under acidic conditions
The in vitro bacterial survival rates were observed after the bacteria were treated with different pH (7.2, 6.3, 5.5, 4.8, 4.3) LB for two hours. All strains (wild type, mutant and complemented) had reduced survival rates as the pH decreased ( Figure 5). However, at the lower pH values (6.3, 5.5, 4.8, 4.3), the survival rates of the mutant were lower than those of the wild type bacteria (P<0.05). The survival rates of the complemented strains were between those of the wild type and mutant bacteria. These results suggest that the eha mutant was more susceptible to acid stress than the wild type bacteria, indicating that Eha is required for free-living ET13 to resist acid stress.

Discussion
In the present study, we examined under what conditions the EHAMPE expressed the strongest Eha activities. On the condition, we compared the whole transcriptomic profiles of the wild type ET13 and its eha mutant, and found 147 DEGs to adapt to acid stress. Moreover, our experiments showed that Eha was involved in protecting the bacteria against acid stress in vitro and against the effects of acidification in vivo.
To survive and adapt to their environment, bacteria integrate multiple signal conduction pathways, such as the two-component system (TCS), to mediate an appropriate cellular response [26]. Our study found that 20 genes that were differentially expressed between the ET13 wild type bacteria and eha mutant are involved in the TCS pathways, seven of which were down-regulated and 13 of which were up-regulated under acid stress. Most interestingly, these down-   regulated DEGs (ETATCC_RS06125, ETATCC_RS06135, ETATCC_ RS06130, ETATCC_RS07080, ETATCC_RS06140, ETATCC_ RS06155, ETATCC_RS06145 and ETATCC_ RS06150) are all related to citric acid metabolism. The DEGs TATCC_RS06140, ETATCC_ RS06135 and ETATCC_RS06130 code for α, β and γ subunits of citric acid lyase C. The DEGs ETATCC_RS06145 and ETATCC_ RS06150 encode citric acid lyase G subunits. Citric acid and its lyase are the intermediate product and critical enzyme of the tricarboxylic acid (TCA) cycle, respectively. The TCA cycle is critical for carbon metabolism and energy generation [27]. In a previous study, a lyase gene mutant E. tarda was unable to proliferate and cause fatal infection in a fish model [28].
The DEGs ETATCC_RS14195 code for glutaminase, and both ETATCC_ RS00005 and ETATCC_RS12655 code for serine dehydratases. Glutaminases in Helicobacter pylori deaminate glutamine to glutamate, resulting in the production of ammonia, which may be of significance in the pathogenesis of H. pylori-associated diseases [29]. The gene ETATCC_RS06155 encodes and antiporter like a gastric proton pump. The pump in H. pylori is implicated in transient hypochlorhydria, which facilitates gastric colonization and potentially triggers the epithelial progression to neoplasia [30]. The ETATCC_ RS16530 gene codes for an outer membrane protein that transports capsule polysaccharides from the periplasmic space to the bacterial surface. Polysaccharides may serve to protect organisms from harsh environmental conditions and to increase their virulence [31].
Our present transcriptomic analysis also discovered that 17 of the DEGs were associated with the biosynthesis of secondary metabolites, and these were significantly up-regulated under acid stress. Secondary metabolites are synthesized using primary metabolites as precursors and regulators, and the process is easily affected by the bacterial environment [32]. ETATCC_ RS10420 and ETATCC_ RS10415 encode tryptophan synthases. The tryptophan synthases in Chlamydia trachomatis regulate the bacterial survival rates in Hela cells and the resistance to IFN-γ [33]. ETATCC_RS06795 encodes threonine dehydratase. The threonine dehydratase in Streptococcus pneumoniae contributes to the bacterial virulence and colonization in pneumococcal infection [34]. The findings suggest that these DEGs may regulate the virulence and survival of E. tarda.  Our transcriptomic analysis also found that Eha was associated with 12 up-regulated DEGs involved in ABC transporter systems, like oligopeptide permease (Opp) in the bacterial membrane that transports nickel, amino acids, thiamine, iron and so on [35]. Nickel enzymes have nickel in their active center, which is required for their biological function. Nickel enzymes include urease, hydrogenase, carbon monoxide dehydrogenase. The Opp1 ABC transporter in Staphylococcus aureus imports nickel and cobalt, and its mutant showed reduced mortality in terms of systemic infection and colonization of the bladder and kidneys in a urinary tract infection model [36]. Eha may regulate ABC transporter systems to control the entry of nickel and the other micronutrients into bacteria, and bacterial urease, hydrogenase and other nickel enzymes may affect the survival rates in different environments and may regulate the pathogenicity of E. tarda. ETATCC_RS11420, ETATCC_RS14855, and ETATCC_ RS14860 encode FliS of the chaperone of flagellin FliC [37] and FlgK and FlgL of the flagellar hook, respectively [38]. The FlgM-FliA regulatory circuit plays a central role in coordinating bacterial flagellar assembly. FliS modulates FlgM activity as a chaperone to control late flagellar gene expression, motility and biofilm formation in Yersinia pseudotuberculosis [39]. A bacterial biofilm is a structured community of bacteria in a self-produced extracellular matrix, which allows the bacteria to adhere to a solid surface or biological tissue. The formation of biofilms apparently improves the ability of the bacteria to resist adverse environmental conditions.
Our RNA-seq analysis revealed that 20 of the DEGs, including ETATCC_ RS13205, are directly involved in bacterial metabolism under different circumstances, and participate in material synthesis and degradation, as well as energy metabolism networks. ETATCC_ RS13205 encodes an acid phosphatase, an enzyme which is associated with the synthesis of bacterial capsule I and IV, which takes part in bacterial resistance [40]. These findings are not only important to better understand the roles of the various DEGs in bacterial physiology, but will help to generally understand the full potential and evolution of protein phosphorylation for signal transduction, protein modification and homeostasis in all bacteria [41].
In summary, the high-throughput RNA-seq was used to examine the whole transcriptomic profile of E. tarda for the first time. Overall, we observed that Eha regulates 147 genes that affect the bacterial energy, metabolism, and virulence, which allow E. tarda to adapt to an acid environment. Therefore, these results could be helpful for further investigations of the molecule mechanism(s) underlying how E. tarda survive in macrophages under acidification.