Shewanella inventionis ZQ14 and its phage vB_SInP-X14 were both isolated from the coastal waters of Qingdao, China (36.063°N, 120.320°E). Diluted 1 ml seawater to 10^–5 step by step and took 200 μL diluted seawater, inoculated it on 2216E solid medium by the plate coating method, and cultured it in an incubator at 28°C and 120 rpm. Shewanella inventionis ZQ14 was purified more than three times in the 2216E medium. Then, the 16s rRNA of the purified Shewanella inventionis ZQ14 was amplified by PCR and sequenced. After searching by BLASTn, the 16SrRNA similarity between the Shewanella inventionis ZQ14 and Shewanella inventionis KX27 (T) was 99.69%.

Phage was isolated by the double-layer plate method (Zhang W. et al., 2021). The same seawater sample was filtered through the 0.22 μm pore-size membranes (Millipore) to obtain viral water. The 200 μl host bacterial (Shewanella inventionis ZQ14) culture in the logarithmic growth stage was inoculated with 200 μl filtered seawater and incubated for 15 min. A single plaque was collected and placed into SM buffer (50 mM Tris-HCl, 0.1 M NaCl, and 8 mM MgSO4 [pH 7.5]) at 4°C. The viral mixture was passed through a 0.22-μm membrane filter and allowed to infect the host again. This step was repeated five times to obtain purified viral stock.

The purified virus was amplified to 500 ml using the reinfection of Shewanella inventionis ZQ14. Then viral culture was passed through a 0.22-μm membrane filter (Millipore) to remove the bacteria. The virus culture was concentrated to 3 ml with an ultrafiltration tube (Millipore Amicon Ultra-15) at 5,000 g (Zhang W. et al., 2021).

The vB_SInP-X14 morphology was examined by transmission electron microscopy (TEM). Briefly, 20 μl purified and concentrated phage solution was negatively stained with phosphotungstic acid (1%w/v, pH 7.0) for 5 min (Haq et al., 2012). The electron micrographs of vB_SInP-X14 were captured by JEOL Model JEM-1200EX TEM at 100 kV.

The 1 ml bacterial solution and 1 ml purified viral solution were incubated at 26°C for 15 min with MOI 0.1, the non-absorbed virus was removed by centrifugation, and the precipitate was washed three times with 1 ml fresh 2216E liquid to remove the remaining unabsorbed phage. Sampling was done every 5 min in the first hour and every 10 min from the second hour. Viral abundance was determined by the method of double-layer agar plaque assay. The number of plaques in different periods was counted to draw the growth curve. Three parallel experiments were conducted for this assay (Cai et al., 2019).

Liquid media of 2216E were prepared with pH of 3, 4, 5, 6, 7, 8, 9, 10, 11, and 12 and sterilized at 121°C for 20 min (Liu et al., 2019, 2021). The 100 μl phage solution (initial titer ∼ 10⁸ PFU/ml) was diluted with 1 ml medium of different pH and placed at room temperature for 2 h. The 200 μl phage solution treated with different pH was mixed with 200 μl of the host bacterial solution and the host was in the logarithmic growth stage, whose concentration was 2 × 10⁷ cfu/ml. After infection for 15 min, the double-layer plate (solid 2216E medium: agar 15 wt.%, semi-solid 2216E medium 7.5 wt.%, agar brand: Solarbio, Beijing) was poured and cultured at a constant temperature overnight at 26°C. Then the number of plaques was counted to draw the phage survival curve at different pH, and carried out three parallel experiments.

Several copies of the same phage solution (titer ∼ 10⁸PFU/mL) were treated at - 20°C, 4°C, 25°C, 35°C, 45°C, 55°C, 65°C, and 75°C, respectively, for 2 h (Zhang et al., 2020). The host bacterial solution in the logarithmic growth stage was vortex mixed with the treated phage solution. After infection for 15 min, the double-layer plate was poured and cultured at a constant temperature overnight. The plaques on the plate were counted and the phage growth curve related to temperature was drawn. We carried out three parallel experiments.

The DNA extraction of vB_SInP-X14 was performed by Virus DNA Kit (OMEGA).

The purified phage genomic DNA was sequenced by Novogene Company (Tianjin, China), using Illumina Miseq 2 × 300 paired-end method. The sequencing library used NEBNext Ultra DNA Library Prep Kit for Illumina (NEB, United States). In short, the 1 μg DNA sample was cut to 350 bp size for Illumina sequencing and further PCR amplification. To ensure the accuracy and reliability of subsequent information analysis results, the raw data was filtered to obtain effective data (clean data), and clean data produced 4,073,908 reads. Clean data were assembled using the software ABySS (v1.3.7) (Simpson et al., 2009) with multiple-Kmer parameters. The GapClose (v1.12) software was applied to fill the gap of preliminary assembly results with default parameters. The fragments shorter than 500 bp were filtered and the final results were obtained. BWA (bwa-0.7.8) software was used to compare the reads to the assembled genome sequence, and SAMTOOLS (v 0.1.19) software was used to interpret and count the coverage of the reads to the assembled genome sequence. Both BWA and SAMTOOLS software used default parameters. The average coverage of reads on the genome was 33,330. Orientation was obtained by using the program BreakDancer (Fan et al., 2014). The genome sequence Shewanella phage vB_SInP-X14 has been deposited in GenBank, whose accession number is MK796797. The cumulative GC skew¹ analysis was used to detect the origin and terminus of replication of the phage genome (Uchiyama et al., 2008; Zhang W. et al., 2021). The prediction of the putative open reading frames (ORFs) was combined with several different software, including prodigal (Hyatt et al., 2010), GeneMarkS sever (Besemer et al., 2001), and RAST (Brettin et al., 2015). All ORFs were annotated by Pfam with default parameters (Mistry et al., 2021), BLASTp (E-value < 1e-5) against by the Non-Redundant (nr) protein database of the GenBank, HHpred simultaneously based on four databases (PDB_mmCIF70_13_Jul, Pfam-A_v34, COG_KOG_v1.0, NCBI_Conserverd_Domains (CD)_v3.18) (Söding et al., 2005). The genome mapping was visualized by CLC Main Workbench 20. The tRNAscan-SE program predicted tRNA sequence (Lowe and Chan, 2016).

To expand the homologous sequences of vB_SInP-X14, genome sequence of vB_SInP-X14 was searched as a query against the NCBI database, searching through BLASTn (E-value < 1e -10, Query Cover > 25%), and against the Integrated Microbial Genomes/Virus (IMG/VR) through Diamond BLASTp (v.0.9.21) (E-value < 1e-10, Query Cover > 50%). Sequences shared more than 50% ORFs of vB_SInP-X14 were considered as candidate sequences (Buchfink et al., 2015).

A total of 65,899 Caudovirales phages were downloaded from NCBI-virus. The pre-experiment of protein sharing between phages was carried out by using vConTACT to find phages linked to vB_SInP-X14. Then, according to the pre-experiment results, only added vB_SInP-X14-related phages and NCBI-RefSeq for cluster analysis to facilitate the data visualization. Accordingly, 4,406 reference genomic sequences of isolated Caudovirales, were selected and combined with the above homologous sequences to vB_SInP-X14 in the IMG/VR dataset. All proteins were compared using all-versus-all DIAMOND BLASTp (E-value ≤ 1E-10, coverage ≥ 50%, amino acid identity ≥ 30%). Markov clustering algorithm (MCL) was used to divide protein clusters (PCs) (Enright et al., 2002). Then vConTACT 2.0 was used to calculate the similarity score between genomes. All predicted proteins were combined as a database with DIAMOND BLASTp; therefore, the database in vConTACT chose “none.” Viral clusters were defined using ClusterONE (Nepusz et al., 2012). To show the detailed cluster analysis result of vB_SInP-X14, the results of vConTACT analysis were screened. Only the classified phage genomes of Podoviridae in the International Committee on Taxonomy of Viruses (ICTV) were selected and combined with homologous sequences to vB_SInP-X14 as the source. The phages that generated weight with the source were selected and displayed in the detailed figure. The network diagram was visualized by Gephi.

All 23 isolated Shewanella phages in GenBank (Table 1) were used to build a whole-genome tree with Viptree² (version 2.0) and analyze their similarities. The Viptree calculates the sequence similarities of the whole genomes according to tblastx, and generates the “proteome tree” of the viral genome sequences (Nishimura et al., 2017). In addition, a separate phylogenetic analysis of vB_SInP-X14 was constructed using Viptree.

A phylogenetic tree based on the whole genome, including at least one member of all genera within Podoviridae family in the ICTV, was constructed. The whole-genome phylogenetic tree including 82 viral isolates was constructed by the Virus Classification and Tree Building Online Resource (VICTOR) (Meier-Kolthoff and Göker, 2017) based on nucleotide, and the treefile of GBDP_Trimming_D6_FASTME was chosen and visualized by iTol (V.5) (Letunic and Bork, 2019). The whole-genome Average Nucleotide Identity (ANI) was computed by fastANI (Jain et al., 2018), then visualized by pheatmap in R. In addition, this result also verifies the results of VICTOR and fastANI. According to the results of the phylogenetic tree and ANI analysis, some viral isolates and representative UViGs similar to vB_SInP-X14 were selected to perform the whole genome gene comparison map with Viptree, and the protein clustering through All-to-all BLASTp analysis was performed by Orthofinder (version: 2.5.2) using the diamond blastp flag of E-value < 1e-5, Query Cover > 50%, and identity > 30% (Emms and Kelly, 2015).

The relative abundances of vB_SInP-X14 were expressed as tpm (transcripts per million) and calculated by CoverM (v0.3.1) (parameters: -m tpm –min-read-percent-identity 0.95 –min-read-aligned-percent 0.75) (Li B. et al., 2010; Li et al., 2021). The quality-controlled reads from the 154 Global Ocean Viromes (GOV 2.0) dataset were mapped to the vB_SInP-X14 contig by minimap2 (v2.17) program. Several other phages were selected as references, including the representative pelagibacter phages (HTVC010P, HTVC019P, HTVC011P), cyanophages (P-SSB7, P-SSM7) (Labrie et al., 2013), Puniceispirillum (SAR116) phage HMO-2011 (Zhao et al., 2013), Roseophage SIO1 (Angly et al., 2009), two phages similar to vB_SInP-X14 (Vibrio phage VvAW1 and Thalassomonas phage BA3) (Efrony et al., 2009; Nigro et al., 2012), and other Shewanella phages. Five marine viral ecological zones (VEZs) defined by the Global Ocean Viromes (GOV 2.0) dataset are Arctic (ARC), Antarctic (ANT), temperate and tropical epipelagic (EPI), temperate and tropical mesopelagic (MES), and bathypelagic (BATHY) (Gregory et al., 2019). The relative abundance of viruses in the five VEZs were added and transformed into log₁₀, and the subsequent results were visualized by heatmap and histogram.

Shewanella phage vB_SInP-X14 was isolated from surface coastal waters of Qingdao, China (36.063°N, 120.320°E), regarding Shewanella inventionis ZQ14 as the host. vB_SInP-X14 formed clear plaques with 2 mm diameter in double-layer agar (Figure 1A). The typical podoviral morphology of vB_SInP-X14 was observed by TEM, which had an icosahedron head (∼59 nm in diameter on average) and a short non-contractible tail (∼24 nm long on average) (Figure 1A).

Experiments of viral one-step growth, pH stability, and temperature tolerance were performed to characterize the biological features of vB_SInP-X14. From the result of viral one-step growth analysis (Figure 1B), the latent period is less than 5 min and the rise period is approximately 40 min. The short latent period of vB_SInP-X14 suggested that the viruses could be rapidly absorbed and injected into cells, and indicated the unique cellular virulence of vB_SInP-X14 (Ma et al., 2021). The growth curve of this trend has also been reported (Chan et al., 2014; Li et al., 2016). In the temperature stability experiment, the phage showed the maximum titer at 35°C. The titer of vB_SInP-X14 was at a high level when the temperature ranged from −20°C to 65°C (Figure 1C); however, the phage titer decreased linearly at more than 65°C, and completely lost its activity at 75°C. The optimum pH range of vB_SInP-X14 is between four and nine (Figure 1D). When the pH exceeds nine, the phage titer decreases sharply and until the pH is 12, the phage loses its activity.

Phage vB_SInP-X14 encapsulated a double-strand DNA genome with a length of 36,396 bp, G + C content of 44.1%. The result of the cumulative GC Skew analysis indicated that the origin and terminus of replication were at the regions 200 nt with the lowest value and 36,300 nt with the highest value, respectively (Supplementary Figure 3). A total of 40 ORFs were predicted, and 24 of which have orthologs with known functional genes in the public databases (Supplementary Table 1). The genome of vB_SInP-X14 is modularized by divergent functions and formed five gene clusters in the genomes. The replication/recombination and transcriptional regulation modules were viral early gene clusters, containing six (ORFs 3,5,11, and 14-17) and five (ORFs 1, 7-9, and 18) ORFs, respectively. The structural and packaging module, lytic module were viral late gene clusters, containing nine (ORFs 27, 28, 31, 32, 33, 35,36, 38, and 39) and three (ORFs 22, 24, and 25) ORFs, respectively (Figure 2A). Transcriptional regulation and replication/recombination modules accounted for large proportions in the genome of vB_SInP-X14. P22-homologous structural and packaging-related proteins were detected in the genome of vB_SInP-X14. Two P22-like packaging genes were annotated, including the P22 tail accessory factor (ORF 35) and the P22-like portal protein (ORF 31). In addition, vB_SInP-X14 contains a P22-like capsid protein (ORF 33).

In the structural and packaging module, ORFs 31, 32, and 33 (locus-tags: X14_000031, X14_000032, and X14_000033) encoded portal protein, scaffolding protein, and coat protein, respectively, which are necessary for procapsid formation of tailed phages. The bacteriophage scaffold protein (ORF32) is an essential protein for the formation of a procapsid (Black, 1995). The intermediate structure of the procapsid is formed by the assembly of scaffold protein and viral coat protein, in which the scaffold protein is located in the coat protein subunit during the arrangement of the icosahedron. In subsequent steps, the scaffold protein molecules are released, leaving shell proteins that cannot be self-assembled (Sun et al., 2000). Then the viral genome is packaged into a preformed procapsid. In the DNA phage, it is composed of two packaging proteins, named terminase large subunit and terminase small subunit (TerL and TerS) (Casjens, 2011). They are assembled into the dodecamer of portal protein and insert DNA into the empty capsid with high precision and efficiency, and the whole process needs to be driven by ATP (Roy and Cingolani, 2012). In the TerL, two or possibly three ATP-binding sites are expected, and the TerL usually shows a specific affinity for portal protein. In contrast, the TerS seems to regulate the activity of the large subunit (Rao and Black, 1988; Black, 1995). Besides, similar to many other phages, the TerS and TerL subunits are adjacent, which are ORFs 27 and 28 in the genome of vB_SInP-X14, respectively (Duhaime et al., 2011).

Three lytic-associated ORFs were encoded by the virus genome, suggesting that the life strategy of vB_SInP-X14 might be lytic, which was in accordance with the results of the one-step growth curve The lytic module includes two holin proteins (ORFs 24 and 25) and one endolysin N-acetylmuramoyl-L-alanine amidase (ORF 22, locus-tag: X14_000022). ORFs 24 and 25 belong to the holin family, and their potential functions appear to be the transport of murein hydrolases across the cytoplasmic membrane to the cell wall where these enzymes hydrolyze the cell wall polymer as a prelude to cell lysis (Schmidt et al., 1996). The potential function of ORF 22 is to dissolve the amide bond between N-acetylmuramoyl and L-amino acids in bacterial cell walls (Hynes et al., 2012).

Apparent modularization was observed for the genome of vB_SInP-X14 on the basis of the ORFs arrangement (Figure 2A).

Phage vB_SInP-X14 is the first isolated podovirus infecting Shewanella, distinct from other Shewanella phages. To identify the taxonomic position of the phage vB_SInP-X14, the phylogenetic tree based on the viral proteome of whole genomes was established by Viptree, including all isolated Shewanella phages (Table 1) and 100 viruses closely related to Shewanella phages (Wang et al., 2021). The result showed that vB_SInP-X14 was far away from other Shewanella phages (Figure 2B). Therefore, vB_SInP-X14 can be clustered with other podovirus within Podoviridae family (Figure 2B and Supplementary Figure 2). The hosts of these podoviruses are Betaproteobacteria and Gammaproteobacteria, while Shewanella is within the Gammaproteobacteria class.

Phage vB_SInP-X14 was also different from other isolated phages in the order Caudovirales. After searching against the entire NCBI database, we found four homologous phages to vB_SInP-X14, which were unclassified at the genus level and not included in the NCBI-ref database. And the genomes of three phages were incomplete. Through searching against the IMG/VR database (Zhan and Chen, 2019), we detected 14 non-redundant homologous UViGs closely related to vB_SInP-X14 (Supplementary Table 2). In the network diagram of the pre-experiment, only phages related to vB_SInP-X14 were presented (Supplementary Figure 1). The network diagram (Figure 3), vB_SInP-X14 has linkages between six isolated phages. Among them, vB_SInP-X14, four isolated phages without genus-level taxonomy information (Thalassomonas phage BA3, Vibrio phage 1.183.O, Vibrio phage 1.184.A, Vibrio phage 1.211.A), and 14 UViGs from IMG/VR were classified into the same viral cluster (VC 0_0). Although Vibrio phage VvAW1 and Rhodoferax phage P26218 have linkage with vB_SInP-X14, they were not classified under VC 0_0 and used as outlier references (Figure 3B). According to these results, we suggest that the members of VC 0_0 could represent a new genus, named here as Bocovirus. In addition, the VC 0_0 is a relatively independent cluster away from other defined viral genera (Figure 3A).

The phylogenetic tree established by VICTOR is an effective tool for viral classification. vB_SInP-X14 and its homologous sequences are novel in the VICTOR phylogenetic tree, after comparing with all the classified Podoviridae phages in the ICTV. The unclassified branch represented by vB_SInP-X14 is distant from other classified viruses (Figure 4). According to the whole-genome phylogeny and OPTSIL classification, 20 phages and UViGs represented by vB_SInP-X14 can be divided into the same subfamily, while 19 of them except Vibrio phage VvAW1 can be grouped into a novel genus (Zhang Z. et al., 2021). In addition, The FastANI heatmap showed that the ANI values of the 19 phages and UViGs ranged from 76.39 to 91.26% (Figure 4), which is higher than the threshold ANI values (> 70%) for the classification of the same genus by the ICTV (Turner et al., 2021).

Virus classification is a process of continuous improvement. At present, there is more and more software to help us classify viruses. Therefore, the combination usage of a variety of methods is necessary to achieve a persuasive result. This study used the method of Turner et al. (2021), combined the results at both nucleic and amino acid levels, and used multiple analysis methods to classify and analyze the vB_SInP-X14.

Exploring homologous genes between different genomes is a common research method for phylogenetic analysis (Brüssow and Hendrix, 2002; Hatfull and Hendrix, 2011). The 19 phages and UViGs represented by vB_SInP-X14 in the whole-genome phylogenetic tree were selected for the detailed genomic and comparative genomic analysis (Figure 5). The protein homology heatmap using all-versus-all BLASTp showed that their gene sharing rates were higher than 40% (Figure 5A). The gene sharing rate of vB_SInP-X14 with other genomes ranged from 54 to 76%. In the protein homology heatmap, two groups have a high gene sharing rate and selected one sequence from each of these two groups (Figure 5A). According to the above results, two incomplete genome sequences with shorter lengths were discarded, and the rest 11 genome sequences were selected for the comparative genomic analysis. The result showed that conserved genes are distributed in almost every module, except lytic (Figure 5B). In the packaging module, including terminase small subunit, portal protein, tail accessory factor, and stabilization protein are conservative in the 11 genome sequences. In structural genes, the 11 genome sequences shared a conserved coat protein. The methyltransferase subunit in the transcriptional regulation module, that is, the domain contained by the largest gene in the whole genome, except Vibrio_ phage_ 1.184. A, all 10 genomes are conserved. In addition, only one lytic-related gene was detected in several sequences, and there were no conserved cleavage genes in these 11 genome sequences, which also shows the particularity of vB_SInP-X14, suggesting their short lytic life cycle (Figure 1B).

The biogeographical distribution of vB_SInP-X14 was characterized in the 154 viral metagenomes from five VEZs of the Global Ocean Viromes (GOV2.0) dataset: Arctic (ARC), Antarctic (ANT), temperate and tropical epipelagic (EPI), temperate and tropical mesopelagic (MES), and bathypelagic (BATHY) (Gregory et al., 2019). The results confirmed that HTVC010P is one of the most abundant viruses in the ocean (Figure 6; Zhao et al., 2013). In addition, Prochlorococcus phage P-SSP7 and P-SSM7 and Pelagibacter phage HTVC019P, HTVC011P, and SAR116 phage HMO-2011 were also abundant (Kang et al., 2013; Zhao et al., 2013). From the perspective of relative abundance, Shewanella bacteriophages are generally not abundant. The relative abundance of vB_SInP-X14 was not high, mainly detected in the MES and ARC, and not detected in BATHY and ANT. Vibrio_phage_1.083.O was the most abundant Shewanella phage, and mainly detected from the MES, BATHY, and ARC, and followed by SppYZU01, SFCi1, 1/44, and vB_SInP-X14. The distribution patterns of SppYZU01 and vB_SInP-X14 were similar. In addition, some Shewanella phages were not detected in any VEZs of the GOV2.0 dataset, including Shewanella phage SppYZU05, Shewanella_phage_Spp001, or Shewanella_phage_S0112 (Figure 6). Although the relative abundances of Shewanella phages are not high in the ocean, considering the ecological significance and sedimental preference of Shewanella and their ability in the dissimilatory metal reduction and the heavy metal treatment (Yang et al., 2020), Shewanella phages should still play significant roles in the ocean through the mediation of their abundance and community structure, especially in the sediments and metal pollution environments.