Endogenization of ssDNA fragment from an infectious virus in the genome of the Chaetoceros tenuissimus

One of the smallest diatoms, Chaetoceros tenuissimus, maintains their population despite coexisting with infectious viruses during blooms. To further understand this relationship, here, we sequenced the C. tenuissimus NIES-3715 genome. A gene fragment of a replication-associated gene from its own infectious ssDNA virus (designated endogenous virus-like fragment, EVLF) was found to be integrated in a total of 41 Mbp of both haploid assemblies. In addition, the EVLF was transcriptionally active and conserved in nine other C. tenuissimus strains from different geographical areas, although the primary structures of their proteins varied. The phylogenetic tree further suggested that the EVLF was acquired by the ancestor of C. tenuissimus. A target site duplication, a hallmark for long interspersed nuclear element retrotransposons, anked the EVLF. Therefore, the EVLF was likely integrated by a retrotransposon during viral infection. The present study used genome information provides further insights into the diatom-virus evolutionary relationship.


Introduction
Diatoms (Bacillariophyta) are an important group of oceanic eukaryotic phytoplankton accounting for approximately 40% of primary marine production 1, 2 . The TARA Oceans project, a global plankton sampling campaign, has highlighted the signi cance of diatoms in global biogeochemical cycles 3,4 .
Research on the dynamics of diatoms is, therefore, important to understand global ecosystems. Generally, the growth and photosynthetic activity of diatoms are higher than those of other phytoplanktons and are often responsible for the blooms found in coastal and upwelling regions, thus playing a role as a major food resource for zooplankton, larvae, and lter feeders 5, 6 among others. Diatom dynamics, in addition to predation, are primarily controlled by environmental factors, such as water temperature, light, and nutrients 5,7 . Other than these abiotic factors, diatom populations are exposed to diverse biological stressors. Over the past three decades, numerous studies have suggested that viral infections are a major determinant of phytoplankton fate in aquatic environments 8 . Indeed, knowledge regarding diatom viruses has accumulated rapidly since their rst report in 2004. Roughly two diatom-virus groups have been identi ed, namely single-stranded (ss) DNA and ssRNA diatom viruses 9 . Diatom cell death due to viral infections can be readily observed in the laboratory, with infected cultures dying off in a few days 9 . In nature, however, the diatom population does not exhibit a rapid decrease in numbers, even in the presence of these infectious viruses; therefore, they seem to have the ability to coexist in the same region 10,11 . Hence, diatoms are thought to have evolved a mechanism to resist viral infections, as has been observed in other phytoplankton-virus systems.
The virus resistance mechanisms of eukaryotic microalgae have been reported by numerous researchers, e.g., variation in susceptibility in the host cell [12][13][14] , blockade of intracellular virus genome replication 15 , bacterially mediated virus resistance 16 , and host cell physiological barriers 17,18 . In addition, recent genomic and transcriptomic analyses have revealed virus resistance systems at higher resolutions. For example, variabilities in genomic islands containing the viral-attachment genes of the cyanobacteria Prochlorococcus facilitate host-virus coexistence 19 . Moreover, for the smallest eukaryotic phytoplankton, Ostreococcus tauri (Mamiellophyceae), the size of the hypervariable region in chromosome 19 differs greatly among strains, which is assumed to be closely related to its susceptibility to infection by its dsDNA virus, OtV 20,21 . Furthermore, in many host organism-virus relationships, whole-or partial-viral genome sequences present in the host DNA can act as a viral resistance factor, e.g., superinfection exclusion and RNA interference related mechanisms [22][23][24][25] . To gain a more in-depth understanding of the host-virus systems, recent studies have highlighted the importance of using host genome analyses.
The marine planktonic diatom, Chaetoceros tenuissimus Meunier (Bacillariophyta, Centrales), is rectangular in the girdle view and is one of the smallest (~ 5 µm) diatoms. This species is widely distributed and is observed in Japanese coastal waters 10 , the Narraganset Bay 26 , the Mediterranean Sea 27,28 , and the San Matıas Gulf 29 . A previous study showed that C. tenuissimus has a high growth rate of at least three divisions per day and blooms during spring and autumn to levels of ~ 10 7 cells/l 10 . To date, four different viruses capable of infecting C. tenuissimus have been isolated and characterised, two different ssDNA viruses (CtenDNAV type-I and type-II), and two different ssRNA viruses (CtenRNAV type-I and type-II) 9,30 . However, the C. tenuissimus population in natural environments sustains its bloom size even in the presence of these viruses 10,11 . The tolerance to viral infection, along with their fast growth rate, suggests that C. tenuissimus might play an important role in maintaining primary productions in coastal environments. Knowledge on the relationship between C. tenuissimus and its viruses has gradually accumulated from the viewpoint of growth-physiology studies based on traditional culture experiments 17,31 , however, studies focusing on the aspect of cell biology are lacking. Here, we have explored the utility of genomic sequencing for this diatom species to further the current understanding regarding host-virus interactions at the molecular level. We believe that novel genomic discoveries can provide critical insights into the evolution of the diatom-virus relationship.

Results
Genome assembly and gene prediction A total of 16.6 Gigabases (Gb) of sequence, providing a 150-fold coverage of the genome sequence, was obtained using the Illumina Miseq and Nanopore MinION platforms and assembled into 41 megabases (Mb) in a total of 85 scaffolds, ranging in size from approximately 1 Kb to 4.46 Mb (Supplementary Table 1). Using the kmer-based statistic, the haploid genome size and its heterozygosity were estimated as 39.7 Mb and 1.56%, respectively (Supplementary Fig. 1), and the size corresponded to the assembly. The accuracy of the genome assembly was con rmed by comparing with the T. pseudonana and T. oceanica genomes using BUSCO software 32 . The 82 complete genes assigned by BUSCO were identi ed in the C. tenuissimus genome assembly, and the number was larger than that for the Thalassiosira species (Supplementary Table 2). Thus, a successful genome assembly of C. tenuissimus NIES-3715 was achieved. A summary of the genome structure is shown in Table 1. A total of 18,705 protein coding genes were predicted in the haploid nuclear genome. Of the predicted proteins, 14,860 had signi cant similarities (e-value < 1e-5) to protein sequences in the non-redundant proteins database (nr), and 9,544 had recognisable Interpro domains. Complete chloroplast and mitochondrial genomes were identi ed from the assembly, with sizes of 116 Kb and 36 Kb containing 131 and 33 predicted genes, respectively (  Supplementary Fig. 3a). In contrast, the mitochondrial genome in this study is the rst reported in a Chaetoceros species (Supplementary Fig. 2b) and showed similarity to the mitochondrial genome of T. pseudonana (83.3% identity) which is classi ed with Chaetoceros as being Centrales (Supplementary Fig. 3b). All sequence reads obtained in this study were deposited in the DDBJ Sequence Read Archive under accession number DRA009158, and the assembly scaffolds for the nuclear genome, as well as the chloroplast and mitochondrial genomes of C. tenuissimus NIES-3715, were also deposited under accession numbers BLLK01000001-BLLK01000085, LC537471, and LC537470, respectively. EVLF in the nuclear genome of C. tenuissimus strains An EVLF, which was similar to a replication-associated protein in C. tenuissimus DNA virus, was found to be integrated between the predicted proteins in the nuclear genome (Fig. 1a). The 181 amino acids were encoded by 546 bases of the sequence and partially aligned with the sequence of the replicationassociated protein in C. tenuissimus DNA virus SS12-43V (67% identity; Fig. 1b). Although EVLF translation was stopped in the middle by a termination codon (TGA), the amino acid sequence could be followed in the next translation frame, and was found to be similar in sequence up to amino acid residue 381 in the C. tenuissimus DNA virus (Fig. 1b), although a termination codon was also located in the middle of the 3′ end (Fig. 1b). A poly-adenine (A) like sequence was observed downstream of the fragment and common sequence "CATAAAA" anked the fragment (Fig. 1a).
PCR analysis using speci c primers designed to amplify the EVLF over its 1.5 Kbp length revealed that nine out of 10 other strains of C. tenuissimus also possessed the EVLF in their genomes (Fig. 2a). However, no target region in C. tenuissimus strain B was ampli ed in this study. In this case, it is possible that the primer set used to amplify the 1.5 kb region did not bind to the target site in strain B because other primer sets used for RT-PCR were able to amplify the EVLF region (Fig. 2a). As a result, strain B was also found to possess the EVLF in its genome. From the sequences of the cloned ampli ed gene fragments, a total of 15 clones with distinct sequences were obtained. The ML tree of these protein sequences formed two clades (Fig. 2b). Moreover, in clade I, three clades were divided with over 93% of bootstrap probability (BP, Fig. 2b). While the 10 sequences in clade I had stop codons and/or frame shifts (Fig. 2b), the translated amino acid sequences were aligned with the other sequences as well as the virus replication-associated protein from C. tenuissimus ( Supplementary Fig. 4). All the EVLFs from the C. tenuissimus strains were clustered with bacilladnavirus replication-associated proteins (Fig. 2c). In particular, the EVLFs were formed by the replication-associated proteins in the C. tenuissimus DNA virus type-V clade as a sister group, with moderate statistical support (BP = 85%; Fig. 2c).

Transcription of EVLF in C. tenuissimus NIES-3715
In the growth phase, the EVLF was found to be transcribed in three replicate cultures of C. tenuissimus NIES-3715; however, in two of the cultures, the EVLFs at 6 days were found to be more weakly transcribed than those at other days, although the actin gene was transcribed constantly during sampling (Fig. 3).
Moreover, as no ampli cation of the negative control was observed, all potentially contaminating DNA was completely digested in the samples (Fig. 3).
Comparative genes among organisms The orthologous genes among C. tenuissimus, T. pseudonana, P. tricornutum, and Cyanidioschyzon merolae were identi ed from a total of 18,705, 11,673, 10,408, and 4,803 protein sequences, respectively. The genes included in the 2,102 ortholog groups were common among these organisms (Fig. 4a). In the C. tenuissimus genes, 3,265 ortholog groups were common in both T. pseudonana and P. tricornutum, while in the other groups a total of 1,011 and 822 were common in T. pseudonana and P. tricornutum, respectively (Fig. 4a). The speci c paralog genes in C. tenuissimus also formed 29 groups (Fig. 4a).
In the ortholog group, OG0000, 385 paralog genes in C. tenuissimus were most abundant among these organisms and its 376 amino acid sequences possessed a leucine rich repeat 5 domain (IPR026906; Supplementary Table 3). Subsequently, seven ortholog groups, OG0004, OG0006, OG0021, OG0024, OG0065, OG0094, and OG0121, noticeably formed groups with genes possessing reverse transcriptase domains (IPR013103 and IPR000477), although the number of genes that contained these domains were not equal to that of the paralog genes (Fig. 4b). In the above ortholog groups, except for OG0006, the genes in C. tenuissimus were more abundant than those of other organisms (Fig. 4b).

Transposable elements
Although the total number of retroelements in C. tenuissimus was less than that in P. tricornutum, three retroelements of C. tenuissimus were detected at higher levels than those in other organisms: short interspersed nuclear elements (SINEs, 0.02%), long interspersed nuclear elements (LINEs, 0.64%), and Gypsy/DIRS1 of long terminal repeat elements (1.11%; Fig. 4d).

Discussion
Diatoms and their infectious viruses seem to be able to coexist in natural waters and are closely related to each other 10,11 . To understand the ecological relationships between host diatoms and their viruses, we sequenced the C. tenuissimus genome. The haploid genome size in C. tenuissimus NIES-3715 was estimated to be 39.7 Mb, and the sequences were assembled into 41 Mb with 85 scaffolds. Although the assembly is not yet at the chromosome level, a total of 18,705 genes were predicted from the haploid assembly. Surprisingly, a predicted gene having a similarity to a putative replication-associated protein of its infectious ssDNA virus was found to be integrated between two coding sites (Fig. 1a). The predicted EVLF protein sequence was highly similar to the mid-region of the replication-associated protein in C. tenuissimus ssDNA virus SS12-43V, which infects and lyses C. tenuissimus NIES-3715. Moreover, a frameshift caused by a two-base deletion, as well as nonsense mutations, was observed in the sequence (Fig. 1b). Therefore, the integration of the virus replication-associated gene appears to have not fully occurred. Like C. tenuissimus, an integrated virus fragment has been shown in the diatom P. tricornutum genome and its transcription has been detected by an EST analysis 34 . However, this fragment is not similar to that in C. tenuissimus but instead is similar to a viral replication gene identi ed only from viral metagenomic data (data not shown). Therefore, this is the rst case of the integration of an extant closely related infectious virus fragment into its host genome.
When a virus rst infects the host, its replication relies on the host cellular machinery 35 . In particular, the genes of DNA viruses must be transcribed as mRNA 35 . The ssDNA diatom viruses encode three putative proteins in their genomes 30 . To integrate the virus gene into the host genome, the nucleic acids of both the host and virus should be in the same cell during the infection. There are several key molecular processes required for the successful integration of a viral gene into the host genome. A Bornavirus-like nucleoprotein gene has been found to be integrated into mammalian genomes, and this integration is thought to be mediated through mRNAs by retrotransposons, which are mobile genetic elements with their own reverse transcriptase (RT), such as LINEs 36 . The banana genome also contains an integrated infectious virus gene, and it is thought to have been integrated by the Ty3/gypsy retrotransposon, which is a retroelement with long terminal repeats 37 . In the C. tenuissimus genome, the number of proteins that possessed a RT domain was higher than that in T. pseudonana, P. tricornutum, and F. solaris (Fig. 4c). Similarly, the percentage of predicted retroelements in C. tenuissimus was higher than that in T. pseudonana and F. solaris (Fig. 4d). Notably, the percentage of LINEs in the C. tenuissimus genome was the highest among the different species (Fig. 4d). In addition, considering that the number of paralog genes possessing or lacking an RT domain in C. tenuissimus was higher than that in other species, it is assumed that the C. tenuissimus genome contains more retrotransposons (Fig. 4b). These results indicated that more retrotransposons have been copied and integrated into the genome, which might be a factor in explaining the large genome size of C. tenuissimus. It is, therefore, possible that the retrotransposons possessing the RT domain converted the mRNA of the virus replication-associated gene into cDNA during the infection process. However, this then leads to the question of how the cDNA of the virus replication-associated gene integrates into C. tenuissimus genome. In general, once retroelements enter the nucleus, an endonuclease encoded by a LINE makes a single-stranded endonucleolytic nick in genomic DNA at a degenerate consensus sequence (5′-TTTT/A-3′, with "/" indicating the scissile phosphate), exposing a poly-(T) tail and a 3′ hydroxyl group that serves as a primer for the reverse transcription of retroelement mRNA, and subsequent integration 38, 39 . The integrated retroelements have a structural hallmark in that they are anked by variable size target site duplications (TSDs) 38, 39 . This structural hallmark was observed in the region of the integration of EVLF in the C. tenuissimus strains ( Supplementary Fig. 5). A poly-(A) like sequence was observed downstream of the stop codons, which corresponds to the terminus of the virus replication-associated gene ( Supplementary Fig. 5). In addition, the endonuclease cleavage sites "5′-TTTTATG-3′" anked the EVLF and were characterised as TSDs (Fig. 1a). In addition, a 5′ truncation and point mutations within the integrated retrotransposon are characteristic of integration by LINEs 39 . These structural hallmarks were also observed in the EVLFs of the C. tenuissimus strains (Fig. 2b) but not in the P. tricornutum genome 34 . In the case of C. tenuissimus, this evidence suggested that the virus replication-associated gene was possibly integrated into the C. tenuissimus genome by host LINEs.
In the culture experiments, CtenDNAV type-II was inoculated into the host cells at the log phase, and the cell density decreased after the cells reached the stationary phase 30 . In short, the inoculated cells maintained a high growth rate, suggesting that C. tenuissimus can resist virus infection 40 . Contrary to C. tenuissimus, in many marine microbial host-virus relationships, the viral mortality of the host increases when there is a high growth rate of host cells [41][42][43][44] . Considering this phenomenon together with the fact that EVLF is transcribed (Fig. 3), we speculate that the EVLF might act as an antiviral immunity-like RNA interference (RNAi). RNA silencing has been shown to act as a defence response against viral infections in plants [45][46][47] and mosquitos 48 , Moreover, this gene silencing mechanism functions in the diatom, P. tricornutum 49 . The C. tenuissimus genome also encoded genes important in RNAi function, such as genes encoding Dicer and Argonaute proteins (accession nos. GFH46084.1 and GFH61989.1, respectively). From the RT-PCR results, the EVLF expression levels for 6 days culture were lower than those for other culture periods. Therefore, there is a possibility that C. tenuissimus represses viral proliferation through the replication-speci c RNAi machinery during periods of high growth rate.
ssDNA viruses infecting C. tenuissimus are classi ed into two types (type-T and type-V) based on the sequences of replication-associated proteins 50 . The EVLFs were certainly derived from the type-V clade of C. tenuissimus DNA virus (Fig. 2c). Although the geographic distribution of these strains is different in Japan ( Supplementary Fig. 6), the EVLFs were highly conserved among strains and in both alleles (Fig. 2b). Moreover, nine out of 10 strains were found to have EVLF at the same locus (Fig. 2a). These results indicate that the ancestor of C. tenuissimus had acquired the replication-associated gene from one type of virus, and then might have acquired an antiviral immunity-like RNAi machinery to resist against infectious viruses. Consequently, the populations have survived to date, and the EVLFs remain encoded in their genome as a fossil. Although the enigma of the survival strategy against the infectious viruses used by C. tenuissimus still remains to be solved, it is hoped that the genome information from C. tenuissimus will shed further light on it.
In this study, we sequenced the C. tenuissimus genome and discovered that a replication-associated gene associated with its own infectious virus has been integrated into the genome. This discovery represents the rst case of an extant, closely related infectious virus fragment being integrated into a host genome; meanwhile, the EVLF may repress viral proliferation by RNA interference during periods of high growth rate. Finally, our analysis of their relationship suggests a close evolutionary relatedness.

DNA extraction and DNA library construction
The C. tenuissimus strain NIES-3715 was isolated from Seto Inland Sea, Japan ( Supplementary Fig. 6)  For long-read sequencing of genomic DNA using MinION (Oxford Nanopore Technologies, Oxford, UK), total nucleic acid was extracted from the pellet using the DNAs-ici!-F (RIZO Inc., Tsukuba, Japan), according to the manufacturer's protocol. To extract genomic DNA from the total nucleic acid sample, the sample was treated with RNase A (Nippon Gene, Tokyo, Japan) and subsequently puri ed with phenol/chloroform prior to construction of a DNA sequencing library. The sequencing library was constructed using the ligation sequencing kit (SQK-LSK109, Oxford Nanopore Technologies) and sequenced by MinION, according to the respective manufacturer's instructions. After sequencing, basecalling was performed with Albacore (v2.3.1, Oxford Nanopore Technologies).

De novo genome assembly
To estimate the genome size and heterozygosity, k-mer counting was performed using the short pairedend reads and the Jelly sh programme 52 . The histogram of 21mer counts was visualised using GenomeScope 53 .
Hybrid assembly of all Illumina short reads and MinION reads was performed using MaSuRCA 54 (v3.3.0) with default parameters. The haploid genome sequence was constructed from the assembled genome using HaploMerger2 55 (v20180603) with default parameters. The assembly quality was evaluated by QUAST 56 . To evaluate the assembly accuracy, single-copy ortholog genes were searched using BUSCO 32 with alveolata_stramenophiles_ensembl datasets.

Gene prediction and annotation
To predict the gene regions in the genome sequence, we rst obtained the complete open reading frames (ORFs) from RNAseq (DNA Data Bank of Japan (DDBJ) Sequence Read Archive under accession number DRA011082). In RNAseq, the de novo assembly procedure was performed following that described by Hongo et al 57 . The complete ORFs were extracted from the assembled sequences and translated using TransDecoder 58 with default parameters. Next, the quality controlled paired-end reads of RNAseq were mapped to the assembled genome using TopHat2 59 with default parameters. Using the mapping data and the protein sequences of complete ORFs, the gene model was predicted using BRAKER2 60 (v2.1.0) with default parameters. Proteins predicted from the gene model were annotated based on their homology to sequences in the nr database from NCBI using the BLASTP programme with a threshold evalue of < 1e-5, and protein domains were found using Interproscan with a threshold e-value of < 1e-5.

Con rmation of an EVLF in the nuclear genome
To analyse the EVLF in the C. tenuissimus genome, we used nine other strains of this diatom species other than strain NIES-3715 ( Supplementary Fig. 6). One millilitre of a stationary growth phase C. tenuissimus culture was centrifuged at 17,400 × g for 3 min at 4°C. The resulting diatom cell pellets were then preserved at − 80°C until analysis. DNA was extracted from stored cell pellets using the DNeasy Plant Mini Kit (Qiagen), according to the manufacturer's instructions. The EVLF was ampli ed using a primer pair, ctEVLFout_v1_F: 5′-GCAAACACGTKTGTTGATATATCGG-3′ and ctEVLFout_v1_R: 5′-CGATCCTCTTGAAGACCCAGT-3′, (Fig. 1). PCR ampli cation was conducted in a reaction mixture with a 20 µl nal volume, containing 0.5 µl DNA, 1 × BlendTaq buffer (Toyobo, Japan), 200 nM dNTPs, 0.2 µM of each primer, and 1 U BlendTaq DNA polymerase. PCR was conducted using GeneAmp PCR system 9700 with the following cycle parameters: 30 cycles of denaturation at 94°C for 30 s, annealing at 55°C for 30 s, and extension at 72°C for 30 s. The PCR products were then electrophoresed on 1% (w/v) agarose ME gels (Wako Pure Chemical Industries, Osaka, Japan), and the nucleic acids were visualised using Midori green nucleic acid stain (Nippon Genetics, Tokyo, Japan). PCR amplicons of approximately 1.5 kb were excised, and their nucleic acids were extracted (NucleoSpin® Gel and PCR Clean-up; Macherey-Nagel GmbH and Co., KG, Düren, Germany). The PCR products were ligated into the pGEM-T Easy vector (Promega, Madison, WI, USA) and transformed into Escherichia coli DH5α-competent cells (Toyobo, Japan). Sequencing was conducted using the dideoxy method with ABI PRISM 3130 Genetic Analyzer (Thermo Fisher Scienti c).

RT-PCR analysis
The axenic algal cultures of C. tenuissimus strain NIES-3715 were grown in a modi ed SWM3 medium under a 12/12 h light-dark cycle of ca. 500 to 600 µM of photons m − 2 s − 1 , using cool white, uorescent illumination at 25°C for 3 days. For the RT-PCR analysis, preconditioned cultures were inoculated into 1-l of fresh SWM3 medium at a nal density of 2.5 × 10 3 cells/ml using a 2-l polycarbonate Erlenmeyer ask (431255; Corning Inc, Glendale, AZ, USA). This experiment was performed in triplicate. The cultures were subsampled at the early logarithmic growth phase (day 1 and day 2) and at late logarithmic growth phase (day 4 and day 7). One each sampling day, diatom cells in the cultures were retained on 0.4-µm polycarbonate membrane lters (Nuclepore membrane; Cytiva). The number of diatom cells on the lters ranged from 10 7 to 10 8 cells per lter, which were frozen in liquid nitrogen and stored at − 80 ºC until analysis.
The retained lters containing the diatom cell samples were cut into small pieces in the TRIzol reagent (Thermo Fisher Scienti c), and total RNA was extracted using a TRIzol Plus RNA Puri cation Kit (Thermo Fisher Scienti c), with PureLink DNase (Thermo Fisher Scienti c) digesting any contaminating DNA, in accordance with the manufacturer's instructions. Moreover, to completely digest any contaminating DNA, DNase treatment was performed using TURBO DNase free kit (Thermo Fisher Scienti c). The quantity of the total RNA was measured using a Qubit RNA HS assay kit (Thermo Fisher Scienti c). cDNA was constructed from 1 µg of total RNA using an oligo(dT) 15 primer and SuperScript IV Reverse transcriptase (Thermo Fisher Scienti c), in accordance with the manufacturer's instructions.
The EVLF sequence was ampli ed from the constructed cDNA using Ex Taq hot start version (Takara, Shiga, Japan) using the following conditions: initial denaturation phase of 98 ºC for 1 min, followed by 30 cycles of 98 ºC for 10 s, 60 ºC for 30 s, and 72 ºC for 40 s. Actin was also ampli ed using the cDNA and the total RNA to be used as a positive and a negative control, respectively. The primers used for this analysis were as follows: ctEVLFin_v1_F, 5′-AAGAAGAAGAGTCGACTGGATCAAC-3′; ctEVLFin_v1_R, 5′-ACAATAACGGTCTCATGATTGAGC-3′; ctActin_F, 5′-CTGGATGTGTTCTTGATTCTGGAG-3′; ctActin_R, 5′-CTTAGACATACGCTCACTGATTCCTG-3′. The amplicon lengths using these primers pairs were 456 bp for the EVLF and 500 bp for actin.

Phylogenetic analysis of EVLFs and virus replication-associated genes
The genome sequence of EVLF was obtained for nine strains of C. tenuissimus ( Supplementary Fig. 6) using the above genome sequence con rmation. To clarify the evolutionary relationship of the EVLFs, a maximum-likelihood (ML) tree analysis was conducted. First, comparing within the strains, the EVLF protein sequences, including all alleles, were aligned with a replication-associated protein from C.
tenuissimus DNA virus SS12-43V (accession no. BBE21064.1) using MAFFT 61 (v7.212) with default parameters. All stop codons and the frame-shifted amino acids in the alignment were removed by manual editing. The best-t evolutionary model for the optimum alignment was calculated using ModelFinder 62 and the Akaike information criterion. The ML tree was inferred from an evolutionary model using RAxML 63 (v8.2.4) with 100 bootstrap replicates. Second, to compare with virus replicationassociated proteins, the virus protein sequences in the NCBI database were retrieved based on similarity to the EVLF using the BLASTP programme. These accession numbers are shown in Supplementary    the "S" and "T" on the boxes refer to the start and terminal codons, respectively. A yellow box indicates a poly-(A) like sequence, and green boxes indicate a "CATAAAA" sequence. Black and gray arrowheads indicate the ampli cation primer sets, ctEVLFout_v1 and ctEVLFin¬_v1¬, which ampli ed the 1508 bp and 456 bp PCR products, respectively. (b) Alignment of amino acid sequences between the EVLF of C.
tenuissimus NIES-3715 and a replication-associated protein in its infectious DNA virus. The alignment was constructed by MAFFT61. Pink boxes indicate consensus amino acid residues. Hyphens and asterisks indicate gaps and termination codons, respectively. A two-nucleotide deletion led to a separation into two open reading frames.

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download.