Pseudomonas aeruginosa PA1 (Lu et al., 2015) and lytic phage PaP1 (Lu et al., 2013) were stored in our laboratory at -80°C in glycerol. P. aeruginosa was cultured aerobically in Luria-Bertani (LB) medium (5 g yeast extract, 10 g tryptone, 10 g NaCl per liter) at 37°C for all experiments. PaP1 particles were collected and purified using CsCl gradient ultracentrifugation (Lu et al., 2013).

For the one-step growth curve of PaP1, we used a modification of the methods of Lu et al. (2013). The early-logarithmic phase cultures of P. aeruginosa (OD600 = 0.2) were infected with phage (multiplicity of infection = 10). After incubation at 37°C for 5 min to allow adsorption, the mixture was centrifuged for 30 s at 13,000 × g. Un-adsorbed phage was removed from the supernatants by washing twice with LB medium. Pellets were resuspended in 5 ml LB, and the cultures were grown at 37°C with shaking at 160 rpm. The total time of adsorption and washing was about 10 min. A total of 50 μL of sample was taken every 10 min, and the number of PaP1 particles was determined using the double-layer agar plaque method. Three biological repeats were performed. The burst time and size were calculated based on the one-step growth curve.

For both microarray analysis and RT-qPCR, total RNA was isolated from four groups of bacterial cultures (0, 5, 15, and 45 min after phage infection) by SV Total RNA Isolation System (Promega, USA, Z3100). The uninfected samples (0 min) served as a control. Microarray hybridization was performed using 5665 custom-designed probes for the P. aeruginosa genome and 342 probes for phage PaP1 (Agilent 8 × 15 K array). The intergenic regions were not covered. Raw data from the microarray experiments was deposited in the GEO database with accession number GSE80117. In the RT-qPCR analysis, cDNA synthesis was prepared with PrimeScript RT reagent kit (TaKaRa Bio; Dalian, China) according to the manufacturer’s recommendations. Quantitative real-time PCR was performed using SYBR Premix Ex Taq II (TaKaRa Bio; Dalian, China). Primers used in this study are listed in Supplementary Table S1. The 16S rRNA gene was selected as the reference gene for normalization.

Phage-infected bacteria or phage-free bacteria were cultured for 40 min (Supplementary Figure S1). Briefly, 1 L of culture from each group was used for each NMR sample, and five replicate samples were prepared per group. The cultures were placed on ice for 20 min, centrifuged at 15, 000 × g for 10 min at 4°C, and the supernatant was removed. The cell pellets were subsequently washed thrice with pre-chilled PBS, followed by centrifugation at 15, 000 × g for 5 min at 4°C. All samples were snap frozen in liquid nitrogen and stored at -80°C until NMR analysis. Briefly, 50 mg of freeze-dried cells per sample was weighed and suspended in 1 mL ultrapure water. Cell suspensions were further sonicated on ice for eight cycles of a 4 s on/3 s off cycling program (Sonics VX-130, USA) and then centrifuged for 15 min at 13,000 rpm and 4°C. The cell lysate was filtered through a 3 KDa ultrafiltration membrane (Mimacon, USA) at 13,000 rpm for 45 min at 4°C. Briefly, 50 μl DSS standard solution (Anachro, Canada) was added to 450 μl filtrate per sample, and centrifuged at 13, 000 × g at 4°C for 2 min before being transferred to the NMR tubes (Norwell, USA).

Metabolic profiling and peak identification were performed by Anachro Technologies Inc. (Wuhan, China) using described methods (Weljie et al., 2006; Ouattara et al., 2012). All of the NMR experiments were performed on an Agilent DD2 600 MHz spectrometer equipped with a triple-resonance cryogenic probe at 298.15 K. Metabolites were identified and quantified using Chenomx NMR Suit (version 8.1, Chenomx, Edmonton, AB, Canada). A line broadening of 0.5 Hz was applied to the free induction decay prior to the Fourier transformation, phasing, and baseline correction. Data were then carefully phased and baseline corrected in Chenomx Processor. All spectra were referenced to the internal standard (DSS) and analyzed by experienced analysts against the Chenomx Compound Library. From the 10 spectra, 48 metabolites were identified and quantified. All information on metabolite concentration were exported to excel and normalized by weight across all parallel samples before use in the multivariable analysis. Principal component analysis (PCA) and partial least squares discriminant analysis (PLS-DA) were performed by pcaMethod Bioconductor and pls package in R. ggplot2 in R was used for data visualization. For quantitative metabolomics, raw data files were deposited to the EMBL-EBI MetaboLights database with identifier MTBLS431. The complete dataset can be accessed at the following url: http://www.ebi.ac.uk/metabolights/MTBLS431.

Microarray data were normalized and analyzed by GeneSpring version GX10.0 (Agilent Technologies). Differentially expressed genes (DEGs; exhibiting at least twofold change in expression) were screened via one-way ANOVA (p < 0.05) with Benjamini–Hochberg correction for multiple comparisons (Hochberg and Benjamini, 1990). Gene co-expression analysis was performed with linear regression, which quantified the co-expression of interaction partners based on gene expression similarity of the microarray data sets. Linear correlation (Pearson’s correlation, cor) was first computed as described previously (Horvath and Dong, 2008). To make a visual representation, the significant interactions (P < 0.001) between phage and host genes were selected to construct the network by Cytoscape (Cline et al., 2007).

In NMR metabolomics, the software MestReNova 9.0 (MestreLab, Santiago de Compostela, Spain) was used to obtain peak areas from the raw spectrum. All information on metabolite concentration was exported to excel and normalized by weight across all parallel samples before use in the multivariable analysis. Principal Components Analysis (PCA) and Partial Least Squares Discriminant Analysis (PLS-DA) were performed by PCA Methods Bioconductor and pls package in R. For each metabolite, the Variable Importance in the Projection (VIP) score was calculated to determine the significantly changed metabolites (VIP score > 1) in phage-infected cells (Chong and Jun, 2005). ggplot2 in R was used for data visualization.

According to the one-step growth curve of PaP1, the latent time is approximately 20 min, and the burst time is approximately 40 min post-infection (Figure 1A). To track the dynamic interaction between phage and host throughout the infection cycle, we examined the transcriptional changes at three time points after phage infection (5, 15, and 40 min) using microarray analysis, with the phage-uninfected host cells (0 min) as the control. These time points correspond to early, middle, and late infection phases, respectively. Supplementary Figure S1 is a simplified scheme of the experimental design.

In three independent experiments, total cellular RNA was extracted from P. aeruginosa PA1 cultures uninfected (0 min) or infected with PaP1 for 5, 15, or 40 min and then used for microarray analysis. A total of 5,655 specific probes were designed for the analysis of PA1 gene expression. Compared with the uninfected host, 7.1% (398/5655) of genes of phage-infected host were identified as DEGs (fold change ≥ 2, p < 0.05), including 349 downregulated and 50 upregulated genes (Figure 1B and Supplementary Table S2). Approximately 85% (339/398) of DEGs were downregulated at the late infection phase, whereas fewer genes were differentially expressed at 5 min and 15 min. This suggests that gene expression of host P. aeruginosa was massively suppressed at 40 min after phage infection.

First, we specifically focused on the pathways inhibited by phage infection. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis was performed on the 349 downregulated genes at the three time points (Kanehisa and Goto, 2000). This pathway enrichment analysis based on kappa statistics and visualization was conducted through ClueGO/CluePedia plugin of the Cytoscape software with default parameters (kappa score threshold = 0.4, over view term = Smallest P-Value, sharing group percentage = 50.0) (Figure 2) (Bindea et al., 2009). The downregulated genes were concentrated in 13 pathways. The three most observed pathways (valine, leucine, and isoleucine degradation; cysteine and methionine metabolism; and arginine biosynthesis) are related to cellular amino acid metabolism. In addition, bacterial chemotaxis, two-component system, nitrogen metabolism, and oxidative phosphorylation were slightly repressed by phage infection. Thus, the inhibition of amino acid metabolism in the host may be important for phage PaP1 replication.

This pattern of bacterial response is similar to other phage–host studies (Lavigne et al., 2013; Leskinen et al., 2016; Zhao et al., 2016b) in which many host genes are repressed. Only a few DEGs (51, Supplementary Table S2) were upregulated, including membrane proteins involved in energy and small-molecule transport, such as narK1, narH, narJ, and fptA.

We further explored phage–host interactions on a genome-wide level using gene co-expression analysis of PaP1 and its host. According to the calculation of correlation (cor), 169 host genes had significant interactions (p < 0.001) with 39 PaP1 genes during PaP1 infection. Remarkably, all 215 gene pairs had cor values < 0, indicating a negative interaction. As shown in the gene interaction network (Supplementary Figure S2), the PaP1 genes occupied almost all centers of subnetworks. These results revealed that phage genes play a central and inhibitory role on host gene expression.

To further investigate the biological functions of the host that are modulated by phage infection, P. aeruginosa Genome Database (PseudoCAP) function analysis was performed on the 169 host genes involved in the network (Winsor et al., 2005). As shown in the circos plot (Figure 3), the top five most enriched PseudoCAP functions included transcriptional regulators; transport of small molecules; membrane proteins; energy metabolism; and translation, posttranslational modification, and degradation.

For the phage genes, both the network map and circos plot revealed that ORF12 had the most host gene partners, followed by ORF102 and ORF27 (Supplementary Figure S2). Gene co-expression network analysis is a powerful approach to detect functional relevance between genes. Genes with a common mediator tend to have a similar expression pattern (Horvath and Dong, 2008). Thus, we inferred that ORF12, ORF102, and ORF27 might be transcriptional regulators during phage infection. To further investigate the impact of these genes on host physiology, studies that include cross-expression of phage genes in the host should be conducted (Zhao et al., 2016a).

The roles of phage proteins on bacterial physiology has been extensively studied (Danovaro et al., 2008; Hurwitz et al., 2013; Ceyssens et al., 2014; Wagemans et al., 2014; Hood and Berger, 2016) and reviewed (Roucourt and Lavigne, 2009) because it is not only an interesting biological phenomenon, but also a rational strategy for identifying potential antimicrobial drug targets (Braff et al., 2016). Liu et al. (2004) found that Staphylococcus aureus DnaI is the target of phage protein 77ORF104. By screening small-molecule compounds that inhibit the binding of 77ORF104 to DnaI from commercially available libraries, Liu et al. (2004) identified two compounds that are active against DnaI and have a very low minimum inhibitory concentration against S. aureus in vitro. Our group also identified GP70.1 from transcriptomic analysis of phage–P. aeruginosa interactions and found that GP70.1 directly binds to the key regulator RpoS to redirect host metabolism (Zhao et al., 2016a). Therefore, further investigation of the interactions of PaP1-encoded proteins ORF102, ORF102, and ORF27 with host targets is necessary.

Metabolites are downstream products of transcriptional and translational processes and thus provide abundant information about the catalytic and regulatory properties of gene products. A ¹H-NMR-based metabolomics approach was applied to complement the microarray data and to investigate phage-induced changes of intracellular metabolites at the late infection phase (40 min) compared with the uninfected host (0 min) (Figure 4A).

Principal component analysis of the normalized NMR data from 10 samples (five replicate samples) shows a separation between control and phage-infected P. aeruginosa, indicating that there were significant changes in metabolites induced by phage infection. Overall, 48 metabolites were identified from the NMR spectra, including 21 amino acids and derivatives, 10 organic acids, 7 nucleic acid components, 3 alcohols, 1 sugar, and 6 other metabolites (Figure 4B). Comparison of the metabolite levels in phage-infected bacteria with uninfected samples based on variable importance in projection (VIP) analysis revealed 12 metabolites with significantly altered levels (VIP score > 1) (Figure 4C).

The majority of the altered metabolites were amino acids or derivatives, which is consistent with the microarray data. Five out of the 12 significantly altered metabolites, including thymidine, lysine, 4-aminobutyrate, ethanolamine, and N-acetylglucose, were increased in the phage-infected samples. The cellular concentrations of NAD⁺ and betaine were the most significantly decreased, with VIP scores >3.

Pathway enrichment analysis was applied in order to investigate the altered metabolic pathways. The results revealed that 48 metabolites detected in the NMR spectra were enriched into 43 KEGG pathways. As shown in Supplementary Figure S3, the significantly changed pathways had a lower p-value and higher pathway impact value. The most enriched pathways were involved in the metabolism and biosynthesis of several amino acids, namely, cysteine, methionine, glycine, serine, threonine, phenylalanine, tyrosine, and tryptophan. This finding is consistent with the microarray dataset and indicates that the cellular amino acid metabolism was shifted after phage replication.

To investigate the cause of metabolic changes in the host, potential AMGs were predicted from the PaP1 genome (accession number HQ832595) (Table 1) using BLASTn (Altschul et al., 1990). Interestingly, the PaP1 gene ORF 110 encodes a thymidylate synthase, which is highly expressed during the infection cycle as detected by microarray data and RT-qPCR (Figure 5). Thymidine (p = 0.0059) was present in significantly higher levels after phage infection (Figure 5), which indicates that it may be an essential substrate for de novo PaP1 genome synthesis.

The effect of phage-encoded AMGs on host metabolism was suggested by two other metabolomics studies. Phage phiKZ encodes seven AMGs (Table 1), which are involved in the nucleotide synthesis, including a thymidylate synthase (De Smet et al., 2016) Consequently, a significant increase of the molecules of pyrimidine synthesis was detected, such as N-carbamoyl-L-aspartate and thymidine. Ankrah et al. (2014) identified viral AMGs from the Global Ocean Survey, which may partially explain the changed metabolites during phage infection, such as amino sugar and nucleotide sugar metabolism.

Although the expression of AMGs may explain some of the changes in host metabolite levels, many others could not be explained, such as the increase in the level of choline and the depletion of NAD⁺ and betaine. We further examined if the observed changes in metabolic pathways were a consequence of differential expression of host genes induced by phage infection. In P. aeruginosa, choline is oxidized to betaine through a two-step process catalyzed by choline dehydrogenase (CHDH) and NAD⁺-dependent betaine aldehyde dehydrogenase (BADH), which are encoded by betA and betB, respectively (Landfald and Strom, 1986; Velasco-Garcia et al., 2006a; Fitzsimmons et al., 2012) (Figure 5). As observed from the microarray and RT-qPCR, both betA and betB were downregulated at 40 min after phage infection, which may inhibit choline catabolism and lead to a decrease in betaine, potentially explaining the metabolic observations in which there was a significant increase in choline (p = 0.0004) and decrease in betaine (p = 0.008) in the phage-infected P. aeruginosa. NAD⁺ acts as a substrate in the reaction catalyzed by BADH, which constitutes the second step of the choline catabolic pathway.

Cellular NAD⁺ was significantly decreased after phage infection (Figure 4C), which further inhibits betaine synthesis. The decreased NAD⁺ might reflect downregulated energy and amino acid metabolism in the host; NAD⁺ is an important product of energy metabolic processes. According to the microarray results, a high percentage of downregulated host DEGs are related to energy metabolism. Therefore, phage-induced changes of host metabolic pathways are also mediated through the differential expression of host genes.

Interestingly, inhibition of the choline–glycine betaine pathway is highly toxic to P. aeruginosa. Choline is transported into cells through transporters BetT-I and BetT-III and oxidized by CHDH to betaine aldehyde, which is further oxidized to glycine betaine by BADH (Belda et al., 2016). Glycine betaine is an important osmoprotectant that helps P. aeruginosa survive during periods of osmotic stress, such as in infected human tissues. Furthermore, betaine aldehyde is highly toxic (Boch et al., 1996; Sage et al., 1997; Velasco-Garcia et al., 2003). Thus, Velasco-Garcia et al. (2003, 2006b) found that disulfiram, a drug used in aversion treatment of alcoholism, can inhibit BADH and could be a potential antimicrobial agent for P. aeruginosa infection. Therefore, the mechanism by which phage PaP1 inhibits the choline–glycine betaine pathway should be further investigated, and may provide a novel strategy for antimicrobial drug design (Liu et al., 2004).