Transcriptomic Analysis Reveals the Dependency of Pseudomonas aeruginosa Genes for Double-Stranded RNA Bacteriophage phiYY Infection Cycle

Summary Bacteriophage phiYY is currently the only double-stranded RNA (dsRNA) phage that infects Pseudomonas aeruginosa and is a potential candidate for phage therapy. Here we applied RNA-seq to investigate the lytic cycle of phiYY infecting P. aeruginosa strain PAO1r. About 12.45% (651/5,229) of the host genes were determined to be differentially expressed genes. Moreover, oxidative stress response genes katB and ahpB are upregulated 64- to 128-fold after phage infection, and the single deletion of each gene blocked phiYY infection, indicating that phiYY is extremely sensitive to oxidative stress. On the contrary, another upregulated gene PA0800 might constrain phage infection, because the deletion of PA0800 resulted in a 3.5-fold increase of the efficiency of plating. Our study highlights a complicated dsRNA phage-phage global interaction and raises new questions toward the host defense mechanisms against dsRNA phage and dsRNA phage-encoded hijacking mechanisms.


INTRODUCTION
Bacteriophage phiYY is a recently identified member of Cystoviridae , which have genomes consisting of three double-stranded RNA molecules, L, M, and S. Phage F6 was the first identified member of this family isolated in 1973 (Vidaver et al., 1973). Cystoviridae contains three dsRNA segments that are located inside a core particle composed of major structural protein, an RNA-dependent RNA polymerase, a hexameric NTPase, and an auxiliary protein (Nemecek et al., 2010;Mantynen et al., 2018Mantynen et al., , 2019. The core particle is enclosed within a lipid-containing membrane. However, until now, there are only seven sequenced dsRNA phages yet. And six of them primarily infect Pseudomonas syringae, whereas phiYY is the only member that primarily infects P. aeruginosa with rough LPS , which could be included into a phage cocktail to inhibit phage-resistant mutants during phage therapy and might be used to kill the rough LPS strains selected after chronic infection with antibiotic treatment (Hocquet et al., 2016). P. aeruginosa is a common opportunistic pathogen that causes infections of bloodstream, urinary tract, burn wound, and airways of patients with cystic fibrosis (Waters and Grimwood, 2018). Moreover, with the emergence of multidrug-resistant clinical isolates of P. aeruginosa, phage therapy has received renewed attention for treating multidrug-resistant bacterial infections (Smith et al., 2017;Waters et al., 2017;Forti et al., 2018;Jault et al., 2018;Bao et al., 2020). Thus, a solid understanding of the phage-host interaction at the molecular level is valuable for the regulation and legislation of phage therapy in the near future (Rohde et al., 2018).
Transcriptomic approach is a powerful tool to study phage-host interactions and has been widely applied in studying phage-host interactions (Lavigne et al., 2013;Chevallereau et al., 2016;De Smet et al., 2016;Zhao et al., 2017;Lood et al., 2020). It leads to a better understanding of the temporal transcriptional schemes, the impact of phage on host genes expression, and the regulatory mechanisms of phages. For example, the transcriptomic data of P. aeruginosa phage pap3 and its host revealed that phage early expressed protein Gp70.1 could affect host gene expression, and further experiment demonstrated that Gp70.1 binds to a global regulator Rpos to control host gene expression . Phages could also use phage encoded auxiliary metabolic genes to generate additional metabolites to support its replication (Zhao et al., 2017). P. aeruginosa phage PAK_P3 significantly depletes bacterial transcripts to facilitate phage replication (Chevallereau et al., 2016). Moreover, giant phage phiKZ could complete its infection cycle without the support from bacterial transcriptional machinery (Ceyssens et al., 2014), whereas another jumbo phage PA5oct relies on the host RNA polymerase to replicate (Lood et al., 2020).
In this study, we report the RNA sequencing (RNA-seq) analysis of dsRNA phage-P. aeruginosa interactions, reveal the gene expression patterns of both phage and host during infection cycle, and experimentally demonstrate that several deferentially expressed host genes are essential for phage life cycle. This study would contribute to the understanding of infection dynamics of dsRNA phage.

RNA-Seq Experiment Design of phiYY Infection
Phage phiYY is effective in infecting the rough strains and has been included in a phage cocktail to limit the emergence of phage-resistance mutant . PAO1r_8 is a phage-resistant mutant selected after dsDNA phage PaoP5 infection (Shen et al., 2018). PAO1r_8 deleted 341 genes (PA1880-PA2220) and encoded 5,229 ORFs. It lost O-antigen because of the deletion of galU gene. It was renamed PAO1r in short in this study and selected as the host strain for RNA-seq analysis (Table S1).
Next, a one-step growth curve of phiYY infecting PAO1r was measured ( Figure 1A). The latency period is approximately 11 min, and the infection cycle duration is approximately 18 min, after which majority of the host are lysed. Thus, we focused on 6 (early), 12 (middle), and 18 min (late) time points after phage was added as representative snapshots of the phage infection cycle.

dsRNA Phage Transcriptomic Profile
The proportion of reads mapped as phage genes was fluctuating around 0.2%-1% ( Figure 1B). This indicates that proportions of dsRNA phage transcripts are quite minor. On the contrary, the mRNA of dsDNA phage could reach a very high percentage. For example, the proportion of phage reads was over 80% and 60% for Acinetobacter baumannii phage phiAbp1  and P. aeruginosa Podovirus LUZ19 (B) Percentage of RNA-seq reads mapped to the phage phiYY genome. Phage reads account for 0.2%-1% of the total reads. S1, S2, S3 represent three biological repeats. (C) Transcriptomic profile of phiYY genes in the infected host. Genes from segment S and M are highly expressed early at 6 min, whereas genes from segment L are highly expressed later at 12 and 18 min. iScience Article (Lavigne et al., 2013), respectively. There are several explanations for the lower phage transcripts. First, dsRNA phages have a very small genome, for example, the sizes of the segments of phiYY were 3,004 (S), 3,862 (M), and 6,648 (L) bp . On the contrary, the genomes of dsDNA phages usually range from dozens to hundreds of kilo bases (Ceyssens and Lavigne, 2010;Shen et al., 2016). Moreover, dsRNA phages encoded proteins to block bacterial transcriptional machines are not reported yet, whereas these proteins are common in dsDNA phages (Roucourt and Lavigne, 2009). For example, bacteriophage Xp10 encoded protein P7 could bind to the host's polymerase to block host RNA transcription .
The phage gene expression pattern can be clustered into early and late expressed genes ( Figure 1C). Most genes from segment S and M are highly expressed at 6 min and decrease thereafter, whereas genes from segment L are highly expressed at 12 and 18 min. For phage phi6, RNA-dependent RNA polymerase within the phage particle transcribes S and M segments directly. However, the host protein YajQ is responsible for the activation of L transcription (Qiao et al., 2010b). Thus, we infer that segment S and M are immediately transcribed using RNA polymerase inside the phage particle and the initiation of L segment transcription might require a host protein(s), such as YajQ (PA4395), which results in delayed transcription of the L segment.

Bacterial Transcriptomic Profile
The bacterial gene expression level was measured by FPKM (Reads Per Kilobase Per Million Read), and genes with log2 fold change values of 1.5 and q values of <0.05 were defined as differentially expressed genes (DEGs) ( Table S2). Totals of 5.37% (281/5,229), 7.88% (412/5,229), and 8.89% (465/5,229) DEGs RNA-seq data were validated by RT-qPCR. Fourteen DEGs were selected and the expression patterns were validated by RT-qPCR. The primers are listed in Table S3, and the RT-qPCR results are consistent with the RNA-seq results (Table 1).
According to the gene ontology (GO) analysis, the host DEGs can be clustered into several groups (Figure 3). Most of the upregulated genes are involved in the pathways of transcription and translation, because these resources might be essential for dsRNA phage replication. Most downregulated genes are clustered as chemotaxis, transcriptional regulators, adaptation and protection, amino acid metabolism, energy metabolism, and carbon compound catabolism.

Determine DEG Genes that Are Essential for dsRNA Phage Infection
We wanted to test whether the DEGs are expressed passively or are functional for either supporting or blocking phiYY life cycle. Thus, we tested the effects of the six most upregulated and six most downregulated DEGs on phiYY infection, by making single insertional deletions of these genes in PAO1r and testing the efficiency of plating (EOP) of phiYY on each mutant. Among the 12 mutants, DPA4571, DPA3337 (DrfaD), and DPA0848 (DahpB) are completely resistant to phiYY, and the deletion of PA4613 (katB) decreased the EOP to 0.019 ( Table 2).
Four of the six most upregulated genes are anti-oxidative stress genes, including alkyl hydroperoxide reductase (PA0848 and PA0140), catalase (PA4613), and thioredoxin reductase (PA0849). And ahpB and katB are important for phiYY infection ( iScience Article more sensitive to oxidative damages than DNA molecules. For DNA phages, the deletion of anti-oxidative genes had a modest impact on phage infectivity. For example, Campylobacter jejuni phage NCTC 12673 has an EOP of 0.15 on a katA mutant compared with that of a wild-type bacterium (Sacher et al., 2018). Thus, we infer that host-encoded anti-oxidative stress machines are essential to protect phage RNA genomes from oxidative damage and are essential for dsRNA phage infection cycles.
Moreover, PA0800 is highly expressed during phage infection and the deletion of PA0800 increased the EOP to 3.5, indicating that this uncharacterized hypothetical protein might be involved in blocking phage infection. Although the function of PA0800 is unknown, its mechanism is a prospect for investigation in the near future.
For the most downregulated genes, four of six of the genes did not affect phage cycle, whereas DPA4571 and DPA3337 completely resisted phage infection. PA3337 encodes an ADP-L-glycero-D-mannoheptose epimerase (rfaD), which is involved in the biosynthesis of the lipopolysaccharide precursor (Drazek et al., 1995). We performed adsorption assay to test the binding of phiYY to DPA3337 and PAO1r. The adsorption rate to PAO1r was~90%, whereas it decreased to~10% for DPA3337 ( Figure S1). These data indicate that the deletion of PA3337 blocks phage adsorption and results in a complete loss of plaques. But its downregulation during phage infection is intriguing, as it may be unable to deter phage infection, seeing as the phage has already entered the bacteria. One possible explanation is to prevent further phage infection, after the first phage entering, which is similar to superinfection exclusion (Abedon, 2015;Labrie et al., 2010). PA4571 is annotated as a cytochrome c oxidase, and the mechanism of PA4571 on phage infection is unknown.
These data indicate that DEGs are not expressed passively but are the results of complicated bacterial defense mechanisms and dsRNA phage-encoded hijacking mechanisms. Some DEGs might be bacterial defense approaches to constrain phiYY infection, such as PA0800. On the contrary, some DEGs might be controlled by phages to assist phage infection cycles, such as the highly expressed anti-oxidative stress genes.
Phages are viruses that are dependent on bacteria for replication. Thus, the host-encoded genes that are essential for dsRNA phage phi6, phi8, or phi2954 infection had been investigated previously. Mindich et al. iScience Article have found that most of the cystoviridae family utilizes particular host proteins for the control of gene expression from the large genomic segment (Qiao et al., 2010b). Different dsRNA phages use YajQ, GrxC, GrxD, and Bcp for the induction of transcription of large L segments, which have guanine nucleotides missing at the first or second five prime positions relative to that found for segments S and M (Qiao et al., 2008(Qiao et al., , 2010a. Moreover, by screening the transposon mutant library, they found mutations affecting pilus formation or LPS composition to be of consequence for susceptibility to infection of different dsRNA phages, because pilus or LPS is the phage receptor (Qiao et al., 2010c).
However, in this study, we tested the effect of DEGs on phage infection cycle and found that some DEGs are essential for dsRNA phage infection, whereas PA0800 hinders phage infection. Since dsRNA phage phiYY only encodes 20 ORFs, it is reasonable to infer that dsRNA phage replication is massively dependent on host machinery.

Bacteria-dsRNA Phage Interaction Network Analysis
To predict the phage genes that might regulate host gene expression, the phage and bacterial gene coexpression analysis were conducted , which identified nine phage genes that have expression pattern correlations with host genes (Figure 4). The results indicate that phage genes, such as putative procapsid protein, RNA-dependent RNA polymerase, and packing NTPase, may play a central role in interacting with host genes.
Many dsDNA phages encode proteins that could regulate the host gene expression (Roucourt and Lavigne, 2009). For example, T4 phage uses Alc protein to cut off host transcription (Kashlev et al., 1993). However, dsRNA phage phiYY has a limited genome with only 20 annotated genes, and most of these genes are involved in phage replication, packaging, attachment, and lysis , with eight genes that have yet to be characterized. To our knowledge, the impact of dsRNA phage-encoded proteins on host expression has yet to be studied. Thus, further studies of the phage genes might be interesting to test whether dsRNA phages could control and manipulate host gene expression and which phage protein(s) could upregulate the expression of anti-oxidative stress genes.

Conclusion
This work represents the first RNA-seq analysis of the dsRNA phage infection processes, revealing dsRNA phage expression patterns and host responses. Furthermore, the phage-host interaction network analysis  iScience Article predicted some dsRNA phage-encoded proteins that may be involved in manipulating host gene expression. Finally, the host determinants of dsRNA phage susceptibility were investigated by studying the impact of significantly changed DEGs on phage infectivity, and the study reveals that phiYY is critically dependent on host machines to successfully replicate. These results raise new questions regarding host defense mechanisms against dsRNA phages and dsRNA phage-encoded hijacking mechanisms.

Limitations of the Study
The molecular mechanism of anti-oxidative genes on protecting phiYY infection has not been investigated yet.
Further research is needed to elucidate whether dsRNA phage could control and manipulate host gene expression.

Resource Availability Lead Contact
Further information should be directed and will be fulfilled by the Lead Contact, Shuai Le (leshuai2004@ tmmu.edu.cn).

Materials Availability
Available from the Lead Contact under the public/local regulation and Material Transfer Agreements (MTAs) for education and research purpose if they are used not for commercial purposes.

Data and Code Availability
The raw RNA-seq reads and analyzed files were deposited in the NCBI GEO: GSE128811.

METHODS
All methods can be found in the accompanying Transparent Methods supplemental file.

Bacterial strains, phages and culture conditions.
The bacterial strains and phages in this work are listed in Table S1. Phage phiYY was previously isolated from the sewage of Southwest Hospital in Chongqing, China. The accession numbers of phiYY genome segments were deposited in GenBank . (Segments L,M,and S are KX074201,KX074202,and KX074203,respectively). The previously sequenced PAO1r_8 was referred in short as PAO1r in this study (Shen et al., 2018). PAO1r lost 341 genes (PA1880-PA2220), including galU and therefore lost O-antigen. P. aeruginosa strains were grown in Luria-Bertani (LB) broth at 37°C. When required, gentamycin (25μg/ml) were used in LB broth.

One-step growth curve of phiYY on PAO1r
One-step growth curve of phiYY was performed as previously described . PAO1r was infected with phiYY with an MOI (multiplicity of infection) of 0.1.
After a 5 min adsorption, the mixed culture was centrifuged for 1 min at 12,000×g. It was then washed twice with LB medium to remove the unadsorbed phages. The pellets were then resuspended in 5ml of LB medium. After which the culture was diluted 10,000 fold to limit post-adsorption of uninfected bacteria. Then, the cultures were grown at 37°C with shaking. Samples were taken every 10 min until 80 min after phage infection. For each time point, the samples were pelleted immediately for 30 s at 12,000 ×g, and the phage titer in the supernatant was determined by double-layer agar plaque method. Three biological repeats were performed.

Preparing samples for RNA-seq and RT-qPCR
10 ml of bacterial PAO1r culture (OD600=0.5) was infected with phage phiYY at an MOI of 10, and 1 ml of uninfected sample was taken as the negative control (0 min sample). The infected culture was grown at 37°C with shaking. 1 ml culture was taken from the infected culture at three time points (6, 12 and 18min) for RNA extraction.
Three biological repeats were performed.

RNA extraction and RNA-seq analysis.
RNA extraction and RNA-seq analysis were performed as previously described . Briefly, total RNA was extracted from each sample (0, 6, 12, and 18 min after phage infection) by SV Total RNA Isolation System (Promega, USA). RNA quality and quantity were checked using Bioanalyzer (Anilent, USA) and RNA 6000 Nano kit (Anilent, USA) .For RNA-seq experiments, rRNA was removed using the Ribo-Zero rRNA removal kit. The cDNA libraries were constructed and sequenced on an Illumina HiSeq 2500 sequencer (Illumina, USA), using paired-end 2-by 150-bp reads.
Bowtie 2 was used to map the reads to PAO1r and phiYY genome, respectively (Langmead & Salzberg, 2012). Tophat2 and Cufflinks 2.2.1 were used to analyze the RNA-seq data and identify differential genes (Trapnell et al., 2012, Kim et al., 2013. The Reads Per Kilobase Per Million Read (FPKM) and the false discovery rate (FDR) (q value) were used to determined gene expression changes. Differentially expressed genes (DEGs) between the two groups were calculated by DESeq, and DEGs were determined with a log2 fold change value of 1.5 and a q value of 0.05.
GO term enrichment of DEGs was analyzed by Blast2GO software (BioBam).

RT-qPCR validation.
13 DEGs were selected for RT-qPCR validation. RT-qPCR was performed using SYBR Premix Ex Taq II (TaKaRa Bio, China). The primers used in this study are listed in Table S3. The 16S rRNA gene was used as the reference gene for normalization, and the expression of each gene was compared with that of the uninfected host using the delta-delta Ct method.

Insertion deletion of the 12 DEGs in PAO1r
12 insertional mutants were constructed by a previously described procedure (Le et al., 2014). Briefly, to knock out PA5471 in PAO1r, a small fragment within PA5471 was amplified by primers 5471-KO-F and 5471-KO-R (Table S3), using PAO1r genome as template. The PCR fragment was digested with EcoRI/BamHI, and ligated into the EcoRI/BamHI digested plasmid pEX18Gm, resulting in pEX-5471. Then, pEX-5471 was electroporated into PAO1r to generate an insertional mutant via a single crossover. The other 11 knock out mutants were generated with the primers listed in Table S3.

EOP experiment and statistical analysis
The EOP experiment of phiYY on 12 insertional deletion mutants was determined by double-layer agar plate assay. Briefly, 10 μl of serial 10-fold dilutions of phiYY was mixed with 200 μl host, and mixed with 5 ml of 0.7% LB agar, and poured on LB agar plate. The phage titer on each strain was calculated after observation of plaques after overnight incubation, and three biological repeats were performed. EOP was calculated as phage titer on host strain divided by phage titer on PAO1r. The statistical analysis was performed using student's t test, and a P value < 0.05 was considered as statistically significant.
Related to Table2. Percent adsorption of the phage was calculated as [(initial titer − residual titer) / initial titer] × 100%. There biological repeats were performed. The asterisks mark P-value of < 0.05 as calculated by Student's t-test.