Twelve different rat stools and nine water samples contaminated by proximity sewage were collected. Rat stools were taken from a dry place one meter from the water sample area. Samples were harvested in November 2015 in La Ciotat, France, at the same GPS location where P. massiliensis LC8 (Levasseur et al., 2016a) and Faustovirus LC9 samples had been collected (Bou Khalil et al., 2016; Cherif Louazani et al., 2017) (N43.181834, E5.614423).

Vermamoeba vermiformis stain CDC19 was used as cell support. The amoebas were harvested after 48 h of culture in homemade peptone yeast extract glucose medium (PYG) when a concentration of 1.10⁶ amoebas/mL was reached. Cells were then rinsed twice in homemade page’s amoeba saline (PAS) and pelleted at 700 × g for 10 min. The amoebas were then re-suspended in the starvation medium (Bou Khalil et al., 2016) at a concentration of 1.10⁶ amoebas/mL. An antibiotic and antifungal mixture with vancomycin (10 μg/mL), ciprofloxacin (20 μg/mL), imipenem (10 μg/mL), and voriconazole (20 μg/mL) was added to the suspension in order to decrease or eliminate bacterial or fungal contamination. A cell suspension of 250 μL per well was then distributed onto a 48-well plate. The samples were then vortexed and 50 μL were added to each well. The rest of the wells served as negative controls by adding 50 μL of PAS. The plate was incubated at 30°C for 4 days in order to monitor any potential cytopathic effect. This co-culture was repeated twice in the same order. When confronted with a high degree of contamination detected in some wells, filtration using 1.2 μm syringe filter (Merck Millipore) was carried out and gentamycin (20 μg/mL) was added 24 h before the second plate of co-culture (sub-culture 1).

End-point dilution was performed in order to clone the virus before its production. To do so, we successively inoculated diluted viral supernatant on V. vermiformis at a dilution factor of 10. End point dilution was assessed for 5 days and the lysis was controlled by inverted microscopy.

For the production and purification processes, 14 infected flasks of 150 cm² (Corning^®, Corning, NY, United States) were pelleted using the Beckman coulter^® centrifuge Avanti^® J-26 XP (Beckman, France) at 14,000 ×g for 30 min (Andreani et al., 2016; Levasseur et al., 2016a). A 25% sucrose gradient was used for the final purification step. After finalizing production, we proceeded with genome sequencing.

Genomic DNA was sequenced on the MiSeq Technology (Illumina Inc., San Diego, CA, United States) using the paired end and mate pair applications. The DNA was barcoded in order to be mixed with 11 other projects for the Nextera Mate Pair sample prep kit (Illumina) and with 16 other projects for the Nextera XT DNA sample prep kit (Illumina).

gDNA was quantified using a Qubit assay with the high sensitivity kit (Life Technologies, Carlsbad, CA, United States) to 131.3 ng/μl.

For the paired end library, dilution was performed requiring 1ng of each genome as input. The “tagmentation” step fragmented and tagged the DNA. Limited cycle PCR amplification (12 cycles) then completed the tag adapters and introduced dual-index barcodes. The library profile was validated on an Agilent 2100 BioAnalyzer (Agilent Technologies Inc., Santa Clara, CA, United States) with a DNA High sensitivity labchip and the fragment size was estimated to 1.5 kb. After purification on AMPure XP beads (Beckman Coulter Inc., Fullerton, CA, United States), the libraries were then normalized on specific beads according to the Nextera XT protocol (Illumina). Normalized libraries were pooled for MiSeq sequencing. Automated cluster generation and paired end sequencing with dual index reads were performed in a single 39-h run in 2 × 250-bp.

A total of 6.6 Gb of information was obtained from a 697,000 per mm² for the density cluster with a cluster passing quality control filters of 94.6% (12,733,000 passed filtered clusters). Within this run, the index representation for Orpheovirus IHUMI-LCC2 was determined to 12.93%. The 1,942,146 paired end reads were trimmed and filtered according to the read qualities.

The mate pair library was prepared with 1.5 μg of genomic DNA using the Nextera mate pair Illumina guide and two libraries were constructed. The genomic DNA sample was simultaneously fragmented and tagged with a mate pair junction adapter. The pattern of the fragmentation was validated on an Agilent 2100 BioAnalyzer (Agilent Technologies Inc., Santa Clara, CA, United States) with a DNA 7500 labchip. The DNA fragments ranged from 1.5 kb to 11 kb with an optimal size at 6.57 and 2.89 kb, respectively. No size selection was performed and 600 and 117 ng, respectively, of tagged fragments were circularized.

The circularized DNA was mechanically sheared to small fragments with an optimal size of 1029 and 1253 bp, respectively, on the Covaris device S2 in T6 tubes (Covaris, Woburn, MA, United States).

The library profile was visualized using a High Sensitivity Bioanalyzer LabChip (Agilent Technologies Inc., Santa Clara, CA, United States) and the final concentration libraries were measured at 5.13 and 5.4 nmol/l, respectively.

In each construction, the libraries were normalized at 2 nM and pooled. After a denaturation step and dilution at 15 pM, the pool of libraries was loaded onto the reagent cartridge and then onto the instrument along with the flow cell. Automated cluster generation and sequencing run were performed in a single 39-h run in a 2 × 151-bp.

Total information of the two flowcells at 6.2 and 7.9 Gb was obtained from a 648,000 and 863,000 cluster density per mm² with a cluster passing quality control filters of 96.1 and 94% (12,144,000 and 15,627,000 passing filter paired reads). Within these runs, the index representation for Orpheovirus IHUMI-LCC2 was determined at 3.16 and 12.43%. The 725,401 and 1,942,196 paired reads were trimmed and assembled with the paired end reads.

Mate pair and paired-end reads were trimmed using CLC Genomics Workbench v7.5¹. De novo assembly of all reads was conducted using 64-word size and 50 bubble size parameters. We obtained 20 scaffolds representing a total size 1,461,620 bp with an average coverage reads ranged from 423 to 551. In parallel, we used an A5 pipeline assembler (Tritt et al., 2012) with standard parameters on 3,884,384 raw reads (paired end reads) representing 621,103,741 nucleotides. We obtained one scaffold of 1,473,699 with a median coverage reads of 295 with a 10th percentile at a coverage of 226. However, two regions of repeats were not completely resolved. Blast alignments of the two different assembling strategies confirmed these two regions and also underlined a high degree of identity between the two methods of assembly (>99%). For these two regions on the A5 assembly, we used GapCloser (Luo et al., 2012) and GapFiller (Nadalin et al., 2012) to fill two gaps and obtained a final single scaffold of 1,473,573 base pairs.

The MAUVE program (Darling et al., 2004) was used to align and determine nucleotide divergence between genomes. BLAST nucleotide online was used to generate dot plots to explore large repeats in the whole genome and in all specific coding sequences. Emboss Explorer was used online using the following different software programs: palindrome of a 200 maximum length², an e-inverted program, an equicktandem for a fast detection.

Gene prediction was computed using Genemarks software (Besemer et al., 2001). We deleted predicted proteins having a size less than 50 amino acids, and 85 predicted protein from 50 to 99 amino acids were detected by Phyre2 (Kelley et al., 2015) as having abnormal tri-dimensional folding and finally were discarded from our dataset. A Blast protein was performed against the non-redundant (nr) protein database (June 19, 2017). Annotation was performed using a combination of Interpro³ version 63.0, a CD-search tool online (Marchler-Bauer and Bryant, 2004) and delta-blastp (Boratyn et al., 2012). Interpro detected 100 transmembrane domain-containing proteins, and with CD-search and delta-blastp they congruently identified domains in 443 proteins.

tRNA prediction was computed online⁴ (Lowe and Chan, 2016) following different standard parameters successively with eukaryotes, archaea and bacteria. We identified orthologous and paralogous genes by using Proteinortho v5 (Lechner et al., 2011) with 60% coverage and 20% amino acid identity and an e-value of 10^-2 as significance thresholds. Moreover, we generated pan-genomic tree on GET_HOMOLOGUES package (Contreras-Moreira and Vinuesa, 2013) using OrthoMCL algorithm with the standard parameters expected for the coverage and e-value. We choose 60% as minimum coverage in Blastp pairwise alignments and 1 × 10^-2 as maximum e-value.

Orpheovirus IHUMI-LCC2 is available in the EMBL-EBI database under accession number LT906555.

All phylogenetic analyses were conducted using the following procedures. Blastp was used to find close homologous proteins. Then, the MUSCLE program (Edgar, 2004) was used to align amino acid sequences. The FastTree program (Price et al., 2009) was computed with standard parameters using the maximum likelihood method with 1,000 bootstrap replicates and the Jones–Taylor–Thornton (JTT) model for amino acid substitution. Phylogenetic trees were then visualized using iTOL v3 online (Letunic and Bork, 2016).

Bacterial contamination is common in viral co-cultures when using stool and sewage samples, with the frequent presence of resistance to the antibiotics and anti-fungal mixtures used. For this, we used a classic mix of antibiotics notably, vancomycin, ciprofloxacin, and imipenem, as previously reported. However, we added gentamicin to our first sub-culture plate, 24 h before inoculating the new plate in order to eliminate resistant bacteria from the stool samples. After three passages on V. vermiformis, lysis referring to the cytopathic effect was detected in some wells. We performed negative staining on the supernatant of a rat stool in well LCC2 and observed particles with an elongated aspect (Figure 1), some of them appear to be irregular, with a concave shape compared to Pandoravirus and Pithoviruses. In contrast, the apex appears to be more similar to Pandoravirus. We named it Orpheovirus.

The length of the Orpheovirus virions range from 900 to 1,100 nm (N = 10) with a maximum diameter of about 500 nm (N = 13). Some virions could reach 1,300 nm in length; this process was sometimes observed in the host cytoplasm. The cork does not seem to seal by a grid as opposed to Pithoviruses and Cedratviruses. We noticed shapes which were similar to the ostiole-like apex observed in Pandoraviruses (Philippe et al., 2013) with a diameter ranging from 70 to 80 nm (N = 8) obstructed by a thick membrane (Figures 2F,I). Nevertheless, in Pandoraviruses the tegument is composed of three layers, each measuring about 20 ∼ 25 nm. The Orpheovirus’ particles presented a dark dense outside layer coated with short, sparse fibrils on their external surface (black arrow). This dark layer is followed by a medium dense space (white arrow) which is in direct contact with the thin inner hyperdense membrane surrounding the viral core cavity containing the nucleic acid (Figure 2I). Altogether, these layers measure ≈40 nm. The replicative cycle of Orpheovirus showed classical stages of infection and replication in amoeba. Briefly, the virus entry by phagocytosis is the start of the cycle, where particles escape the phagosomal process. DNA delivery occurs in the amoeba cytoplasm via the ostiole-like apex (Figures 2A,B). An eclipse phase takes place at 4 h post entry. Functional viral factories (Figure 2C) are well installed and detected around 14–16 h post-infection. Similar forms corresponding to early virion synthesis are also observed, as it is the case for Pithoviruses and Cedratviruses (Figure 2D). At 20 h post-infection, the host cells’ cytoplasm is fully occupied by newly synthesized virions (Figures 2E,G,H). We were also able to detect viruses outside the amoeba due to cell burst or viruses exiting by exocytosis. Complete cell burst occurred between 24 and 38 h post-infection. This slow viral cycle is often observed in the case of V. vermiformis used as cell support, which is not the case when using Acanthamoeba spp. (Reteno et al., 2015; Andreani et al., 2017).

Orpheovirus has a circular genome estimated at 1,473,573 base pair (including 100 N due to an incomplete elucidate region in its genome) with a GC%-content established at around 25% (Table 1). A megablast or a simple blastn against the nr/nt nucleotide collection database revealed no match with other known giant viruses. Dot plots show various areas of repeats (Supplementary Figures S1–S3). We found 57 palindromic sequences, 1,527 tandem repeats and 832 inverted sequence candidates. The number of repeats explains the complexity observed during the genome assembly steps. A comparison with other giant viruses (Supplementary Table S1) showed an extremely high number of tandem repeats and inverted repeats for Orpheovirus.

1,512 genes were predicted but, following our method, 313 genes with an abnormal conformation already cited in the material and methods section were discarded. We only retained 1,199 genes, resulting in a coding density of around 66.4% (979,005 base pairs). This value is close to that of Pithoviruses but lower than that of Cedratvirus A11. A Blast against the nr database retrieved 509 matched proteins with at least one known protein (≈42.5% of all predicted proteins), and 690 unmatched, which are classified like ORFans genes (≈57.5% of all predicted genes). Of the 509 proteins, two had a hit with unclassified sequences, 148 had a best hit with a virus (≈12.3% of all conserved proteins), 176 with eukaryotes (≈14.7%), and 183 with prokaryotes (≈15.3%) (Figure 3). Regarding the 148 best hits with viruses, we observed 27 best hits with P. sibericum, 11 with P. massiliensis, 15 with Cedratvirus A11, 24 with Mimivirus A, B, and C lineages, and 18 with Klosneuvirinae. Hence, the highest best hit viral was obtained with the putative family Pithoviridae with 53 best hits, although the value was also important with Mimiviridae and associated extend family.

Despite this, the 57.5% of Orpheovirus’ genes are ORFans, with an e-value cut-off of 10^-2. This could increase at ≈66% when we chose a more stringent cut-off value for the blastp at 10^-5. We found 343 genes that formed 167 clusters of paralogous genes in Proteinortho, regrouping few numbers of genes ranging from two to a maximum of four by cluster. Annotation of these clusters revealed predominant predicted proteins mainly as MORN-repeat (55 sequences), Ankyrin-repeat (37 sequences), and F-box domain-containing (81 sequences).

The Orpheovirus annotation presented translation system components as follows: eight aminoacyl tRNA synthetases (aaRS), four translation factors: three initiation factors and one release factor (Supplementary Table S2). Surprisingly, Orpheovirus didn’t present any tRNA. We used the aminoacyl tRNA synthetase, which appeared to be a good way to distinguish and classify some lineages and to describe hypothetical common ancestor (Abrahão et al., 2017; Schulz et al., 2017). The Glycyl-tRNA synthetase of Orpheovirus was found to branch with the Asgard Glycyl-tRNA synthetase, and not with Catovirus CTV1 nor Klosneuvirus KNV1 homologs. This Asgard superphylum described by metagenomic studies seems to be a controversial bridge between prokaryotes and eukaryotes (Spang et al., 2015; Da Cunha et al., 2017; Zaremba-Niedzwiedzka et al., 2017). Regarding the phylogenetic analysis of each tRNA synthetase (Supplementary Figures S4–S11), we observed different patterns for amino acyl tRNA synthetase. While some are monophyletic with other described giant viruses, others appear to be polyphyletic resulting from potential lateral gene transfer.

First of all, we searched for five ancestral genes of NCLDV (Colson et al., 2013) encoding the major capsid protein, the helicase-primase (D5), the DNA polymerase elongation subunit family B, the DNA-packaging ATPase (A32), and the viral late transcription factor 3 VLTF3. Four of these genes were found with the exception of the A32-like packaging ATPase, which was absent in all four viruses (Legendre et al., 2014; Andreani et al., 2016). As reported for Cedratvirus A11, P. massiliensis, and P. sibericum, Orpheovirus presented two distinct RNA polymerase II subunit 1. Multiple ribonucleases such as Ribonuclease R, two Ribonuclease III and one ribonuclease HI were detected in Orpheovirus. Orpheovirus presented glycosyltransferase and numerous proteins involved in lipid pathways. We also identified two proteins presenting multiple fusion bacteria domains involved in Riboflavin (Vitamin B2) biosynthesis. Indeed, we observed in ORPV_596 Tri-functional domains of “Di-Hydro-Folate-Reductase/deoxycytidylate deaminase/Riboflavin biosynthesis protein RibD” presenting a homology with Indivirus ILV1. And the second protein is ORPV_666, annotated like Tri-functional domains “3,4 dihydroxy-2-butanone 4-Phosphate synthase/GTP cyclohydrolase II/Lumazine synthase (RibA+RibB+RibH)” presenting homologies with Indivirus ILV1, Bacillus subtilis, and Acanthamoeba castellanii strain Neff. Vitreschak (2002) demonstrated that Riboflavin operon gene fusion is frequently found in bacteria (Vitreschak, 2002).

After that, phylogenetic analysis based on the DNA polymerase B protein, VLTF3 and RNA polymerase II subunit 1 showed deep branching with Cedratvirus and Pithoviruses (Figure 4 and Supplementary Figures S12, S13). Moreover, 58 reciprocal best hit proteins were only shared between Orpheovirus IHUMI-LCC2, Cedratvirus A11, P. sibericum P1084-T and P. massiliensis LC8. In addition, 14 reciprocal best hit proteins were found to be shared between Orpheovirus and Pithoviruses (14+58) and 15 between Cedratvirus and Orpheovirus (15+58) (Figure 5), while Cedratvirus shared 151 proteins (58+93) with Pithoviruses. Meanwhile, 946 of 1,034 protein clusters (≈91.4%) are unique to Orpheovirus, 319 of 497 clusters (≈64.2%) to Cedratvirus A11, and 114 clusters of 543 (≈21%) to Pithoviruses. There were only two colinearity blocks and nine lines connecting Orpheovirus to other viruses (Supplementary Figure S14).

Following the discovery of this divergence between the four viruses, we decided to investigate Orpheovirus position in the “Megavirales” order further with the help of a parsimonious pan-genomic tree (Supplementary Figure S15). The long branch length observed for Klosneuvirus, Pandoravirus inopinatum and Orpheovirus is explained by the large genome size and by the number of predicted proteins compared to the close relative strains in the tree. These long branches could be a new common marker to explain the emergence of new viral family or lineages in the proposed “Megavirales” order. In the case of Orpheovirus, the pan-genomic analysis confirms this distant relation with the proposed Pithoviridae family.

Orpheovirus has a predicted gene a 434 amino acids protein that we called V21-like protein. This protein had homologs in Blastp, respectively, at 86% coverage, 21% identity with Sputnik virophage V21 protein (La Scola et al., 2008), and 85% coverage, 27% identity with Zamilon (Gaia et al., 2014). These two homologous proteins are annotated as hypothetical proteins, and showed no other homology using the blast strategy. However, HHpred online (Supplementary Data Sheet 3) and Phyre2 (Supplementary Figure S16) detected homology between the V21-like Orpheovirus sequence, Sputnik virophage V21 protein, the Zamilon protein, and a putative transferase present in the genomes of Mimivirus lineages A, B, and C. This V21-like protein also shared a common ancestor with all Sputnik virophages, and Zamilon virophage (Figure 6). No transposase or other mobile elements could be detected, no other special interest homology with other proteins was detected although a Ribonuclease III such as that in MIMIVIRE (Levasseur et al., 2016b) was present near this V21-like sequence in the genome of Orpheovirus.