Persistent mirusvirus infection in the marine protist Aurantiochytrium

Domain: earth_science
Published: 2025-11-11
Editorial outcome: major_revision

Abstract
Mirusviruses are abundant and broadly distributed double-stranded (ds) DNA viruses recently discovered in marine metagenomic data. Their host range and ecological impact are unclear. The protist Aurantiochytrium limacinum possesses two mirusvirus-like genomic elements, one a circular episome (AurliV-1) and the other (AurliV-2) a chromosomal integrant. Here we show that genes in both genomes are expressed and viral particles containing mainly AurliV-1 DNA are produced under starvation conditions and when cells are cultured in standard growth medium. We detected viral particles of ~140 nm in the nucleus, in cytoplasmic vesicles, between the plasma membrane and cell wall, and in the extracellular environment. Of 67 AurliV-1-encoded proteins detected using proteomics, 45 are enriched under starvation conditions, including the structurally important major capsid and triplex proteins. Our results establish Aurantiochytrium as a model system for elucidating mirusvirus-host interactions and demonstrate persistent viral infection in a microbial eukaryote.

Introduction
Elucidating interactions between eukaryotic microbes and the viruses that infect them is critical for understanding microbial ecology on a global scale. The past twenty years has seen rapid growth in our knowledge of such viruses, in particular the nucleocytoplasmic large DNA viruses (NCLDVs). Often referred to as ‘giant’ viruses due to their atypically large virions (>300 nm in diameter) and genomes (0.3 to >2 megabase pairs) (e.g. refs. 1–3), NCLDVs belong to the phylum Nucleocytoviricota4 and have been shown to influence the biology and ecology of diverse algae and protists. For example, NCLDVs impact the persistence of harmful algal blooms by pelagophyte algae5,6. In Arctic green algae, co-infecting lytic prasinoviruses compete with one another in a temperature-dependent manner, an observation with implications for polar food webs in changing oceans7. The nuclear genomes of diverse eukaryotic taxa bear witness to past and present interactions with viruses. These genetic ‘footprints’ range from viral detritus in species of Acanthamoeba8,9 to giant endogenous viral elements (GEVEs) in green algae10, brown algae11 and fungi12, to domesticated mobile genetic elements such as classical retrotransposons and the recently discovered polintons13. The existence of stable, long-term associations between microbial eukaryotic hosts and viruses in nature remains an open question.

The 2023 discovery of mirusviruses underscores the power of metagenomics as an engine of exploration. These dsDNA viruses were inferred from metagenome-assembled genomes (MAGs) pieced together from Tara Oceans data14,15. Mirusvirus genes involved in DNA replication and gene expression share sequence similarity with viruses in the Nucleocytoviricota. In stark contrast, mirusvirus genes encoding structural proteins (including the distinctive HK97-fold major capsid protein (MCP)) are related to viruses in the Duplodnaviria, which includes tailed bacteriophages and herpesviruses. Mirusviruses are predicted to be abundant and active in ocean surface waters from pole to pole14 and mirusvirus genetic signatures have recently been detected in freshwater samples16. However, because of their MAG-based discovery, the physical nature of mirusviruses and the range of hosts they infect is unclear (beyond the fact that they include cells <5 μm in size). Mirusvirus-like MCP genes have been found in the genome of the phagotrophic green alga Cymbomonas17, and surveys of genomic data have identified mirusvirus-like sequences in diverse protists from across the eukaryotic tree of life18. Given this apparently broad distribution, elucidating the biology of mirusviruses is essential to understanding the evolution of microbial eukaryotes and their pivotal ecological roles in aquatic ecosystems.

Mirusvirus genomic elements were recently shown to be a prominent feature of thraustochytrids18,19, protists well known as aquatic decomposers and biotechnological workhorses used to produce polyunsaturated fatty acids for nutritional supplements20. The genome of the marine thraustochytrid Aurantiochytrium limacinum contains two ~300 kilobase-pair (kbp) mirusvirus-like elements19. The first is a circular episome initially dubbed CE1, and the second is integrated into a sub-telomeric region of chromosome 15 (LE-Chr15), reminiscent of the episomal and integrated versions of herpesviral genomes in metazoans21,22. Here we show that A. limacinum produces mirusvirus particles under starvation and when grown in the rich media used for routine culturing. Starvation induces a dramatic shift in mirusvirus protein expression from the episomal genome and triggers changes in host cell biology. Our results suggest that mirusviruses persistently infect A. limacinum and other thraustochytrids, and that dsDNA viruses can have long-term, intimate, hitherto unrecognized associations with protists in nature. Investigation of such interactions in laboratory cultures used for research and biotechnological applications promises new insights into their ecological and evolutionary impact.

Results
Mirusvirus genes are expressed and encode functionally diverse proteins
Initial characterization of the A. limacinum CE1 and LE-Chr15 mirusvirus elements, herein referred to as AurliV-1 and AurliV-2, respectively (for Aurantiochytrium limacinum viral elements 1 and 2), showed them to be extremely divergent from one another and from other known viral genomes, including Tara Oceans-derived MAGs14,19. Seventy-five % (128/170) of the AurliV-1 ORFs in the original Joint Genome Institute (JGI) short-read assembly had no hits in standard protein-based BLAST searches, and only 48 of the original 152 JGI-predicted AurliV-2 proteins retrieved an AurliV-1 homolog as a significant hit. To facilitate more detailed functional analysis of A. limacinum AurliV-1 and AurliV-2, we revisited the ORFs of both mirusvirus genomes based on new gap-free, nanopore-based, consensus sequences. Figure 1a shows a physical map of AurliV-1 with an expanded set of 218 ORFs cross-referenced against JGI gene models (when present). One hundred and seventy-three of these ORFs (~78%) were found to have significant functional annotations using BlastP, InterPro, and/or HMMER searches (Supplementary Data 1). For AurliV-2 (Fig. 1b), 215 ORFs (79%) had such annotations (Supplementary Data 2). All predicted proteins were analyzed using Chai-1 (and in some cases AlphaFold (Version 3) to inform downstream structure- / function-based investigations (see below). Combined with the fact that neither AurliV-1 or AurliV-2 show obvious signs of genome erosion19, the putative viral functions of many of these proteins suggest that both mirusvirus genomic elements represent current and possibly active infections of their host and could yield viral particles. A. limacinum could thus serve as an experimental system with which to explore mirusvirus-host interactions.

Transcriptomic data generated in support of the original A. limacinum genome sequencing project show that many predicted AurliV-1 genes and some AurliV-2 genes are transcribed when cells are grown in rich liquid medium (790 By+)19. Here we performed an independent RNA-seq experiment that yielded similar results, confirming gene expression from both mirusvirus elements when A. limacinum is grown on agar plates containing 790 By+ or mGPY media (Supplementary Fig. 1). AurliV-1 expression was significantly greater than AurliV-2 in both media and was comparable to that of the host genome (Supplementary Fig. 1a). All AurliV-1 genes in the JGI gene catalog showed evidence of expression, 49 of which had read counts above 500 arbitrary units (A.U.). In contrast, only 10 genes from AurliV-2 had read counts above 500 A.U. and some were not expressed at all. Furthermore, whereas genes for the AurliV-1 MCP, triplex proteins 1 and 2, and portal protein were found to be expressed in both 790 By+ and mGPY (MCP (i.e., Aurli1_48830) is the 13th most highly expressed gene, Supplementary Fig. 1b), transcripts for the AurliV-2 homologs of these key mirusviral structural proteins were negligible or completely absent. The expression of AurliV-2 may be suppressed under the conditions examined.

Of the top 10 most highly expressed genes from AurliV-1, five could not be annotated. Remote homology modeling using Phyre2, HHPred, and AlphaFold3/FoldSeek23–26 failed to provide statistically significant functional predictions for the top two most expressed genes (Supplementary Data 1; Supplementary Fig. 1b; both were detected by proteomics (below)). Of the five with confident annotations, none had potential virus-associated functions, with the exception of arylsulfatase Aurli1_41229 and protein kinase Aurli1_48848 (see below). In contrast, all but two of the top 10 most highly expressed AurliV-2 genes could be assigned virus-associated functions, most likely involved in manipulating host signal cascades, such as serine/threonine protein kinase Aurli1_142988. Viral protein kinases are a key feature of large DNA viruses such as herpesviruses, poxviruses, and baculoviruses, and are considered an important aspect of their evolutionary success27. The 5th most abundant AurliV-2 transcript, Aurli1_13274 (Supplementary Fig. S1c), was annotated as a RING-type E3 ubiquitin transferase. Some viruses, such as herpes simplex virus type 1 (HSV-1), encode RING-type E3 ubiquitin ligases, which they use to evade and manipulate host immune systems28. Overall, the functional ambiguity of many of the most highly expressed AurliV-1 and AurliV-2 genes underscore the novelty and divergence of mirusvirus proteins.

Aurantiochytrium limacinum produces mirusvirus particles
Under transmission electron microscopy (TEM), electron-dense virus-like particles 130–150 nm in diameter could readily be seen in the nucleus and cytoplasm of cells subjected to starvation (artificial seawater, ASW) (Fig. 2). The nuclear capsids, ~90 nm in diameter, were generally hexagonal and of variable electron density, consistent with the presence of different stages of virus maturation.

The nuclei of virus-producing cells (Fig. 2b, d, e) were ultrastructurally distinct from those lacking obvious viral particles (Fig. 2a, c); the former typically contained more peripherally located nucleoli and sometimes contained fibrous structures (Fig. 2b, e, Supplementary Fig. 2d, g). Unusual nucleolus-like structures with demarcations resembling procapsids were also observed (Supplementary Fig. 2b, c, e, f), as well as mature capsids very close to – and possibly in the process of budding from – the nuclear membrane (Supplementary Fig. 2h, i). Cytoplasmic virus particles appear to acquire a tegument- or envelope-like structure in discrete, electron dense regions of the cytoplasm, and accumulate in large numbers between the plasma membrane and the cell wall under starved conditions (Fig. 3a, f, g, Supplementary Fig. 3a–f). The presence of viral particles in extracellular material was also confirmed by TEM and negative staining of ultracentrifuged supernatant, where mature virions could be observed embedded within, and tethered to, layered organic material presumably derived from the A. limacinum cell wall (Fig. 4, Supplementary Fig. 4).

Collectively, these observations are remarkably like those made in the early 1970s by Kazama and Schornstein, who described “herpes-type” viruses associated with a thraustochytrid (at the time considered to be a fungus) grown under starvation conditions29,30. These authors studied a strain identified as Thraustochytrium sp., which is today considered to be a distinct genus within the same larger monophyletic clade to which Aurantiochytrium belongs. Kazama and Schornstein thus appear to have been the first to observe mirusviruses, at the time assumed to be herpesviruses (phylum Peploviricota), and long before the discovery of many members of the Nucleocytoviricota.

Unexpectedly, in addition to the starvation conditions studied by Kazama and Schornstein, we also observed virions in control cells grown in rich media (Fig. 2f, Supplementary Fig. 2, and Supplementary Fig. 3g–i). It is thus possible that virion production is associated not only with nutritional status but also specific cell types. The thraustochytrid life cycle is complex and dynamic – in Aurantiochytrium viewed under scanning electron microscopy (SEM), ovoid vegetative sporangia, limaciform (‘slug’) cells, and flagellated zoospores (produced from zoosporangia) can all be seen depending on culture conditions (Supplementary Fig. 5). Importantly, we observed viral particles in various cell types as well as embedded in the cell wall of a putative zoospore (Supplementary Fig. 2a, j, Supplementary Fig. 4d, e). This suggests a complex relationship between the virus and the life cycle of its host. Future work should explore the extent to which viral replication is linked to specific host cell types or states.

Nuclear capsids and extracellular virions contain mirusvirus DNA and proteins
We used quantitative PCR (qPCR) and droplet digital PCR (ddPCR) to determine whether the extracellular virions observed by TEM contain AurliV-1 and/or AurliV-2 DNA. Substrate for these experiments was prepared by (i) centrifugation of cultured biomass grown under fed (control) and starved conditions at low speed to sediment cells, (ii) passing the supernatant through a 0.45 μm filter to remove remaining cell-sized debris (and possibly clumped virions), and (iii) DNase treatment to remove contaminating DNA not protected inside viral particles. qPCR and ddPCR primers targeted unique regions of the AurliV-1 and AurliV-2 mirusvirus DNA polymerase genes as well as two non-viral genes on chromosome 15, i.e., Ubco (encoding a putative ubiquitin-conjugating enzyme) and LeuA (which encodes an alpha-isopropylmalate synthase) (Supplementary Table 1). Ubco and LeuA are not part of the endogenized AurliV-2 region of chromosome 15 and thus serve as a control for host gDNA contamination. The qPCR results suggest that both mirusvirus genomic elements are encapsidated, and that the AurliV-1 DNA polymerase gene is approximately three orders of magnitude (103) more abundant than that of AurliV-2 (Supplementary Fig. 6a).

We next used qPCR to examine whether the inferred relative abundance of AurliV-1- and AurliV-2-derived virions varied across treatments, again using the DNA polymerase gene as a molecular marker for mirusvirus. The relative abundance of the DNA polymerase gene in the fed and starved samples, with and without DNase treatment, was normalized against a spiked-in luciferase-containing plasmid (to control for DNA extraction variation). The results demonstrate that levels of AurliV-1 and AurliV-2 DNA polymerase genes remained relatively consistent at the 24 h time point (Supplementary Fig. 6b). These results suggest that a larger number of extracellular virions contain DNA from the AurliV-1 episome compared to AurliV-2.

Unexpectedly, the fed (control) and starved cells did not differ significantly in terms of extracellular mirusvirus abundance. A more sensitive and robust ddPCR experiment confirmed the presence of AurliV-1 and V-2 DNA polymerase genes in both fed (790 By+) and starved (artificial seawater, ASW) conditions. As with qPCR, extracellular AurliV-1 DNA was found to be much more abundant than AurliV-2 (e.g., under fed conditions at 72 h, ~105 copies of AurliV-1 DNA per μL of supernatant compared to <102 copies per μL of AurliV-2; Fig. 5a, b). AurliV-2 DNA was less abundant than that of AurliV-1 under both conditions. The amount of host genomic DNA contamination, as measured using primers that target LeuA, was minimal in all samples, and DNase treatment removed contamination from unprotected DNA.

Finally, proteomics was used to positively identify AurliV-1- and AurliV-2-encoded proteins in whole-cell pellets, nuclear-enriched pellets, and concentrated supernatants. From a total of 6,101 identified proteins in the whole-cell fraction, 67 were AurliV-1-encoded, as opposed to only nine from genes on the chromosomally integrated AurliV-2 mirusviral element (Fig. 6a, Supplementary Data 1 and 2). Some AurliV-1-derived proteins were detected in all three fractions (e.g., MCP, triplex proteins, a protein kinase, and a peptidase), others in two of the three fractions, and eight were found only in the supernatant (Fig. 6b, Supplementary Data 1). The supernatant fraction was dominated by MCP followed by triplex 2 and 1 proteins (Fig. 6b), consistent with the herpesvirus capsid hexon31, and an AlphaFold3 structural prediction of the AurliV-1 MCP supports the existence of a hexameric structure in vivo (Fig. 6c–f). Notably, AurliV-2 proteins were not detected in the supernatant fraction as seen for AurliV-1 proteins. Given our ddPCR results (Fig. 5), it is conceivable that AurliV-2 DNA is co-packaged at a low level with the more abundant AurliV-1 in capsids composed entirely of proteins encoded by AurliV-1.

We also compared the proteomes of A. limacinum cultures grown in triplicate under control (fed) conditions (790 By+ medium) to those of cultures subjected to starvation (i.e., grown in ASW). Of the 67 detected AurliV-1-derived proteins, 45 were found to be differentially expressed under starvation conditions (Fig. 6a). Most striking among those enriched in starved cells were MCP, triplex proteins 1 and 2, a putative chromosome partitioning protein (ParA), and two distinct arylsulfatase proteins (Fig. 6g; see below). AurliV-1 protein abundances were consistently higher in the starvation treatment compared to control cells, supporting the hypothesis that nutrient stress triggers an increase in production of mirusvirus particles from the AurliV-1 episomal element. The exception among the proteins shown in Fig. 6g is kinase ORF729, which is not differentially enriched between the two conditions and may not be a component of viral particles. Two of nine AurliV-2-encoded proteins were found to be differentially expressed, in this case depleted in ASW relative to standard 790 By+ medium. Neither of these AurliV-2 proteins could be assigned a predicted function.

Persistent mirusvirus infection in Aurantiochytrium limacinum and other thraustochytrids
We have shown that A. limacinum cells grown in replete media and under starvation conditions express AurliV-1 and AurliV-2 genes and produce virions. Our data suggest that A. limacinum AurliV-1 is the predominant mirusvirus genomic element yielding mature virions under starvation and ‘standard’ growth conditions. It is nevertheless possible that AurliV-2 is active to a much lower extent and/or under different conditions. To explore the possibility that the host is limiting the expression of the endogenized AurliV-2 epigenetically, we performed a genome-wide screen for 5-methylcytosine in the Oxford Nanopore sequence data generated herein. We found no obvious differential methylation patterns between AurliV-2 and AurliV-1 under standard growth conditions, and in fact no clear differences in methylation levels between these two mirusvirus genomes and the rest of the A. limacinum nuclear genome (Supplementary Fig. 7). Viral DNAs recently integrated into the nuclear genomes of the amoebozoan Acanthamoeba8 and Amoebidium (a single-celled relative of animals) are silenced by heavy methylation, a mechanism suggested to enable controlled coexistence of viral DNA in eukaryotic nuclear genomes32.

In terms of base composition, the AurliV-1 episome has the same GC content as the A. limacinum nuclear genome (45.0%), whereas the endogenized AurliV-2 is slightly lower at 41.6%19. This is consistent with the hypothesis that AurliV-2 became associated with A. limacinum more recently than AurliV-1, and yet is somehow under tighter host control, perhaps due to its chromosomal location.

Why have mirusviruses gone largely unnoticed in A. limacinum and other thraustochytrids? In control (fed) cells, only ~3% of 397 cell sections observed under TEM contained at least one obvious virus-like particle, and only a single cell was observed with >5 such particles (Fig. 2f). Under starvation conditions, the number of observed virus-producing cells were similarly low: ~3% and ~4% of intact, starved cells in TEM sections were found to have at least one viral particle at 19 h (n=693) and 26 h (n=514), respectively, and between 1% and 1.5% of the sampled population had >5 obvious viral particles (Fig. 2f). Consideration of TEM sections alone undoubtedly underestimates the number of viral particles per cell, but given the large number of cells examined (~400–700 cells), these data are consistent with the hypothesis that in culture a relatively small number of A. limacinum cells produce mirusvirus virions in substantial quantities at any given time. This represents an interesting parallel to the recent discovery of an endogenous giant (NCLDV-type) viral element in lab strains and natural isolates of the green alga Chlamydomonas, which produce mature virions in late-stage growth in only a small number (up to 3%) of cells even during peak viral activity33.

While the endogenous AurliV-2 genomic element is present in all cells, it is formally possible that AurliV-1, whose genome is represented in the original A. limacinum JGI assembly at ~9X higher coverage than AurliV-219, is present in only a subset of cultured cells but as a high-copy episome. The high coverage of AurliV-1 was confirmed here by re-sequencing and assembly of the A. limacinum genome using the latest nanopore sequencing chemistry. The same ~9:1 AurliV-1:AurliV-2 depth of coverage ratio was again recovered, showing that the stoichiometry of the two mirusvirus genomic elements is stable across different sequencing platforms (454 pyrosequencing and nanopore sequencing) and between different laboratory cultures grown years apart. At present, the diversity of A. limacinum cell types in culture is a barrier to further investigation of mirusvirus genome biology at the population level, and long-term synchronization of specific cell types is presently not possible. Vegetative cells grow and divide within sporangia, and even zoospores can exist as clusters (Supplementary Fig. 5)34. Assuming these technical challenges can be addressed, application of single cell ‘omics approaches to the study of A. limacinum should be very informative. Interestingly, putative circular mirusviral genomes similar in size to AurliV-1 were recently inferred from draft genomic data for two other thraustochytrids assigned to Schizochytrium spp. and the chlorophyte green alga Ostreobium quekettii18. This suggests that episomal mirusvirus genomes may be common in nature and may impact the growth and production of strains used commercially to produce fatty acids20.

Mirusvirus appears to represent a persistent infection of A. limacinum, and perhaps other thraustochytrids, similar to, but distinct from the phenomenon of latency in animal-infecting herpesviruses22; mirusviruses replicate constitutively in cells grown under standard conditions and become a more prominent, if sporadic, feature of cells subjected to starvation. How mature virions are dispersed – and how (and if) they infect naïve hosts – is unclear. In virus-containing starved cells (but less so in fed cells), viral particles were commonly seen accumulating in large numbers in the space between the plasma membrane and cell wall (Fig. 3a, f, Supplementary Fig. 3), and only rarely in isolation between cell wall layers and in the extracellular milieu (Fig. 4a, b, Supplementary Fig. 4). A potential clue to escape comes from the AurliV-1 genome, which encodes three distinct arylsulfatases (Supplementary Data 1), one of which contains a proline- and threonine-rich, BLLF1-type/herpesvirus-like outer envelope glycoprotein extension on its amino-terminus (Supplementary Fig. 8). While this BLLF1 domain-containing arylsulfatase was not detected in our proteomic data, the other two (which lack obvious extensions) were enriched under starvation (Fig. 6g). Three AurliV-1-encoded alpha-galactosidase-like proteins also have BLLF1 domains, but on their carboxy-termini (Supplementary Fig. 9).

Bioinformatic analysis revealed complex evolutionary histories for both sets of enzymes. Genes encoding arylsulfatase and arylsulfatase-like proteins also reside within AurliV-2 and several other host nuclear chromosomes (1, 2, and 11), as well as in the thraustochytrid Hondaea fermentalgiana; these proteins belong to a monophyletic clade and appear to be related to bacterial alkaline phosphatases (Supplementary Fig. 8). In the case of the AurliV-1 alpha-galactosidases, homologs with and without obvious glycoprotein domains can be found on A. limacinum chromosome 6, in H. fermentalgiana and other thraustochytrids, and in select fungi and protists (Supplementary Fig. 9, Supplementary Information File 1). Regardless of the origins of these AurliV-1-encoded proteins, given that thraustochytrid cell walls are composed of sulphated polysaccharides rich in galactose, xylose, mannose, and other sugars35, surface-localized arylsulfatase and alpha-galactosidase activities could serve to degrade the A. limacinum cell wall, thus allowing virions to be released from sporangia. They could also facilitate mirusviral infection by digesting cell walls from the outside in. The AurliV-1 genome is predicted to encode a variety of other carbohydrate-degrading enzymes as well, including a putative alpha-amylase and a pectinesterase (Supplementary Data 1). There is clearly much to learn about mirusvirus-host ‘cross-talk’ at the biochemical and cellular level.

The persistent and productive nature of mirusvirus ‘infection’ in A. limacinum suggests a stable, asynchronous virus-host relationship, one with significant implications for understanding viral evolution and ecology in marine environments. The relative impact of the AurliV-1 and AurliV-2 genomic elements on host cell biology is unclear. As recently shown for Chlamydomonas spp.33, concerted exploration of thraustochytrid strain diversity in nature should provide insight into when and how mirusviruses first became woven into their biology, and the extent to which other viruses play a role. In that vein, it is important to note that although not identified in the A. limacinum strain examined herein, lytic single-stranded RNA viruses have been shown to infect the thraustochytrid Schizochytrium sp. strain NIBH N1-2736,37 and a lytic large DNA virus (SmDNAV) infects Sicyoidochytrium minutum38. We thus cannot discount the possibility that multiple distinct DNA and RNA viruses co-occur in at least some thraustochytrid species and genera. This may complicate efforts to disentangle mirusvirus-host interactions but also offer deeper insight into viral ecology and co-evolution.

Our work on A. limacinum has revealed intriguing parallels between the biology of mirusviruses and herpesviruses. In addition to the presence of episomal and sub-telomerically-integrated viral genomes in both mirusviruses and latent herpesviruses19,21,22, this includes (i) angularized capsids in the nucleus with and without electron-dense centres, (ii) margination of chromatin and nucleoli in virus-producing nuclei, and (iii) virus ‘budding’ at the nuclear envelope. Elucidating the mechanisms by which mirusvirus capsids are packaged in the nucleus, acquire a lipid envelope, proceed to the cell surface, and exit the cell is an important task for the future. Equally important will be functional characterization of the mirusvirus-encoded proteins that underpin the viral ‘lifecycle’ and, in the case of A. limacinum at least, what environmental conditions trigger competition and/or synergy between AurliV-1 and AurliV-2.

Metagenomics data predict that mirusviruses are abundant in the sunlit ocean14 and their genomes have been detected in freshwater ecosystems in Japan16; genomic data also suggest that they may be broadly distributed among thraustochytrids18,19. Now that we know how to look, the recognized geographical distribution and host range of mirusviruses will no doubt continue to grow. Establishment of A. limacinum as a model lab system for studying mirusvirus-host interactions provides a foundation for future research with evolutionary, ecological, cell biological, and biotechnological implications.

Methods
Cell culture
Aurantiochytrium limacinum (ATCC MYA-1381) was maintained in ATCC 790 By+ medium (1.0 g/L yeast extract, 1.0 g/L peptone, and 5.0 g/L D-glucose dissolved in filtered artificial seawater (ASW; 423 mM NaCl, 8.99 mM KCl, 9.28 mM CaCl2·2H2O, 22.9 mM MgCl2·6H2O, 22.5 mM MgSO4·7H2O, 2.14 mM NaHCO3)) either as liquid culture grown on an orbital shaker (VWR 12620-926 at 125 rpm), or on solid medium with the addition of 1% (w/v) agar. Cultures were incubated at 27°C, with routine subculturing every 10–14 days (liquid media) or every 5–6 weeks (solid media). For RNA-seq (below), cells were also grown on mGPY (30 g/L glucose, 10 g/L peptone, 5 g/L yeast extract, 26 g/L sea salt) agar plates.

Starvation and fed treatments
To simulate starvation conditions, two-day-old cultures of A. limacinum grown in nutrient-rich 790 By+ were pelleted (5000 × g for 5 min), washed once with ASW, and cultured in ASW on an orbital shaker at 27°C until harvest. Control (fed) cultures were maintained as liquid cultures grown in 790 By+ medium on an orbital shaker at 27°C. Harvest timepoints were as follows: 21 h for SEM; 0 h, 19 h, and 26 h for TEM; 24 h supernatant harvest for qPCR; 0 h, 6 h, 18 h, 24 h, and 72 h supernatant harvests for digital droplet PCR; 72 h for proteomics.

Electron microscopy
A. limacinum cells were pelleted (5000 × g for 5 min) and fixed in 2.5% glutaraldehyde prepared in 0.1 M sodium cacodylate buffer at 4°C for 90 h. Samples were rinsed three times with 0.1 M sodium cacodylate buffer, fixed for 2 h with 1% osmium tetroxide, then rinsed with distilled water. For transmission electron microscopy, samples were stained overnight at 4°C with 0.25% uranyl acetate, dehydrated using a graduated series of acetone (50%, 70% twice, 95% twice, 100% twice), embedded with Epon Araldite resin (3:1 dried 100% acetone to resin for 3 h, 1:3 dried 100% acetone to resin overnight, 2 × 3 h 100% resin), and cured at 60°C for 48 h. Thin (80–100 nm) sections, obtained using a Reichert-Jung Ultracut E ultramicrotome, were placed on 300 mesh copper grids and further stained with 2% aqueous uranyl acetate and lead citrate. Samples were imaged with a Hamamatsu ORCA-HR digital camera attached to a JEOL JEM 1230 microscope at 80 kV, or with a Gatan 832 ORIUS SC1000 digital camera attached to a FEI Tecnai-12 microscope at 80 kV. For scanning electron microscopy, samples prepared as described above were dehydrated in a graded series of ethanol (50%, 70% twice, 95% twice, 100% twice). Drying was performed using the LEICA CPD300 critical point dryer, following the manufacturer’s “cell protocol”. The dried samples were then attached to SEM stubs using double-sided carbon tabs. Samples were sputter coated with 5 nm of gold/palladium (80/20) using the LEICA EM Ace 600. Imaging was performed using the Zeiss SIGMA 300 SEM with an SE2 detector operating at 5.0 kV. For negative stain electron microscopy of mirusvirus particles, samples collected from sucrose (20%) gradient ultracentrifugation (103,538 × g (average), 29,000 rpm) for one hour at 8°C in a Beckman Coulter SW32 Ti rotor) were stained with aqueous 0.25% uranyl acetate.

RNA-seq and whole genome re-sequencing
For RNA-Seq, sample preparation was performed as described by Horowitz et al.39. A. limacinum cultures grown for 24 h in liquid culture were plated onto 790 By+ or mGPY agar plates and incubated at 28°C. After three and four days of growth, agar plates were scraped with a plate-spreader, collected in a 1.5 mL microcentrifuge tube, and frozen. Three biological replicates from each treatment and time-point were collected. RNA was extracted from frozen samples (see above). Quality control was performed using an Agilent Tapestation and Thermo Fisher Qubit 2.0 fluorometer. Library preparation, poly(A) selection, and Illumina HiSeq sequencing were performed by GENEWIZ LLC (NJ). Each sample yielded 812,500 (median) reads from each end. Cutadapt typically removed 8.6% (median) of total bases and trimming by fastQC typically left 712,500 147 base reads from each end.

For long read nanopore sequencing, A. limacinum was cultured for three days in 50 mL 790 By+ medium. Genomic DNA (gDNA) was extracted using a DNeasy Blood & Tissue Kit Quick-Start Protocol (QIAGEN, #69506) with an overnight Proteinase K digestion. Agarose gel electrophoresis (1%) was used to visually assess and confirm the integrity of high molecular weight (20+ kb) DNA. The DNA quantity was assessed using a Qubit 2.0 Fluorometer with the dsDNA broad range assay kit. A nanopore (PromethION, Oxford Nanopore Technologies (ONT)) sequencing library was prepared using the ONT ligation sequencing kit (SQK-LSK114) and ~1 μg of gDNA. A total of 350 ng of prepared library was loaded on a single PromethION R0.4.1 chemistry flow cell (FLO-PRO114M) and sequenced via Oxford Nanopore Technology’s MinKNOW software (v24.06.15). Basecalling was performed with Dorado (v0.8.0) using super accurate model v5.0.0, 400 bps and trimming adaptors. Middle adaptors were removed by Porechop ab initio (v0.5.0) and data filtered using Filtlong (v0.2.1) for 20.6 Gb of reads above q20 and longer than 1 kb. The resulting data were assembled with Canu v2.2 and improved by Medaka (v1.8.0). The long-read assembly was polished using Pilon (v1.22) with Illumina short-read sequence data generated using genomic DNA from cultures grown in liquid 790 By+ medium. The Illumina HiSeq 1000 sequencing (2 × 150 bp) performed by GENEWIZ LLC (NJ) generated 16.5 Gb of paired-end reads available at the NCBI database as BioProject PRJNA1242385. The Illumina reads were trimmed with Trimmomatic (v0.39) to trim 15 bases from the start of read, to trim bases below q20 at either end, and to cut the read when the sliding window average quality dropped below q25. The cleaned reads were then assessed with FastQC (v0.12.1) before being used for Pilon polishing. The AurliV-1 and AurliV-2 sequences were manually checked for sequence errors and submitted to NCBI (PX058987 (AurliV-1) and PX058988 (AurliV-2)).

Viral DNA isolation
A. limacinum cultures grown in 790 By+ for 24 h were kept at 27°C for an additional 24 h or subjected to starvation. Samples were then centrifuged (5000 × g for 5 min) and the supernatant was transferred to a new tube and centrifuged again (21,000 × g for 5 min). The supernatant was then passed through a 0.45 μm filter. Viral DNA extraction was performed as described by Pringle et al.40. Briefly, 180 μL of the filtered supernatant from each treatment was mixed with 20 μL of 3 mg/mL DNase I (Sigma-Aldrich, 69182) and incubated at 37°C for 60 min to remove unprotected DNA. Viral DNA extraction was performed using the DNeasy Blood and Tissue Mini Kit (QIAGEN, 69504) according to the manufacturer’s instructions but with the following modifications: Lysis Buffer AL was spiked with 5 μg/sample of salmon sperm DNA (Invitrogen, 15632011) to improve DNA binding to the column, and 100 pg/sample of a luciferase (luc2)-containing pGL4.26 plasmid (Clontech) to monitor variation during DNA extraction.

Quantitative PCR (qPCR)
Quantitative PCR primers targeting the mirusviral DNA polymerase gene on both AurliV-1 and AurliV-2 as well as two non-viral Chromosome 15 genes (a ubiquitin-conjugating enzyme (Ubco) and alpha-isopropylmalate synthase (LeuA)) were designed using the PrimerQuestTM tool and synthesized by IDT. Primer sequences are provided in Supplementary Table 1. PCR amplification was conducted using the GoTaq qPCR Master Mix (Promega, A6001), following the manufacturer’s protocol across all treatments, with and without DNase treatment, in two technical replicates. The entire setup was repeated four times to ensure consistency and reliability of the results. Additionally, qPCR was conducted on the spiked-in luciferase-containing plasmid to monitor for any inconsistency or variation in the extraction process40. The relative abundance of each gene, across treatments and DNase conditions, was normalized against the luciferase control and the Fed-DNase+ reference using the ΔΔCT method.

Droplet digital PCR (ddPCR)
Aurantiochytrium limacinum was seeded into 30 mL cultures containing 790 By+ and grown to a density between 1 × 106 and 1.5 × 106 cells/mL. Cultures were pelleted (3220 × g for 5 min) and resuspended with 30 mL of either fresh 790 By+ (fed, control medium) or ASW (starvation medium), then maintained until the supernatants were collected (0 h, 6 h, 18 h, 24 h, and 72 h). The collected supernatants were processed for viral DNA isolation as described earlier. The copy number of extracellular mirusvirus was quantified by performing ddPCR on the extracted DNA using the Bio-Rad QX200 system with EvaGreen Supermix (Bio-Rad, 1864034). Each 20 μL reaction contained 10 μL of 2× EvaGreen Supermix, 5.5 μL of template, and 250 nM of the forward and reverse primer. Nuclease-free water was used as the no template control. Droplets were generated using QX200 Droplet Generation Oil for EvaGreen (Bio-Rad, 1864006) and loaded onto a ddPCR 96-well plate (Bio-Rad, 12001925). Amplification was done with a ramp rate of 2°C/second: 95°C for 5 min., then 50 cycles of 95°C for 30 s, 59°C for 30 s, and 72°C for 60 s. Primers targeting the DNA polymerase gene in both AurliV-1 and AurliV-2 were used as a proxy for measuring viral copy number, and the non-viral chromosome 15 host gene LeuA was used as a proxy for measuring host DNA contamination. Primers were designed using Primer3 and synthesized by IDT. Primer sequences are provided in Supplementary Table 1.

Proteomics
Cell culture. A whole cell, data-independent acquisition (DIA), proteomics experiment was performed in triplicate as follows. A two-day-old A. limacinum culture (grown in 790 By+) was used to seed three cultures in 790 By+ (control) and three cultures in ASW (starvation treatment). All six cultures were grown for three days at 27°C and spun down at 3220 × g for 5 min. A viral enrichment sample was made from an ASW-grown culture supernatant that was filtered through a 0.22 μm filter and concentrated on a 50 kDa MWCO centrifuge filter (Millipore Sigma, UFC9050). A pellet from ASW-grown A. limacinum was also enriched for cell nuclei by incubating in nuclear lysis buffer (250 mM sucrose, 10 mM NaCl, 10 mM KCl, 10 mM MES-KOH (pH 5.4), 0.5% NP-40)41 for 10 min, aided by gentle sonication. After 10 min, cells were spun down and the process was repeated until supernatant ran clear. Both viral and nuclear enriched samples were used for data-dependent acquisition proteomics.

Sample preparation for data-independent acquisition (DIA) analysis. Our protocol was adapted from Hughes et al.42. Protein was extracted from cell pellets using 150 μL of a solution of 20 mM Tris-Cl pH 7.5, 150 mM KCl, 5 mM MgCl2, 0.5% (v/v) NP-40, 0.5% (w/v) sodium deoxycholate, and 0.5X cOmplete protease inhibitor (Sigma). A 10 μL nuclease solution (1 U/μL benzonase) was added to the lysate and incubated at 24°C for 10 min. After nuclease digestion, 50 μL of a solution of 600 mM Tris Cl pH 7.5, 8% (w/v) SDS, and 10 mM dithiothreitol (DTT) was added to the lysate and incubated at 95°C for 5 min in a heating block. After cooling to 24°C, lysates were alkylated using 20 mM chloroacetamide (final concentration) and subsequently quenched with DTT. Protein concentration of the resulting lysate was measured using an EZQ assay (Thermo Scientific). Protein cleanup prior to digestion was carried out using an adapted version of the SP3 protocol43. Briefly, 10 μL of a 50 mg/mL carboxylate-coated bead stock (Cytiva) was added to 40 μg of protein in 100 μL of solution followed by 400 μL of acetone. After a 5-min incubation at 24°C, mixtures were centrifuged for 5 min at 5000 × g. Resulting supernatants were discarded and pelleted material was resuspended in 800 μL of 80% ethanol and pelleted again as above. The ethanol supernatant was removed and 1 μg of trypsin/LysC (Promega, V5071) in 100 μL of 100 mM Tris-Cl pH 7.5 with 1 mM CaCl2 was added. After digestion, the supernatant was collected after centrifugation at 12,000 × g for 2 min and subsequently desalted using C18 HLB SPE (Waters) micro-elution plates. SPE wells were initially conditioned with 100% (v/v) acetonitrile (100 μL) followed by equilibration with 0.1% (v/v) trifluoroacetic acid (TFA, Thermo Scientific) (100 μL). Loaded peptides were rinsed using 4% (v/v) acetonitrile in 0.1% (v/v) formic acid (Thermo Scientific) (200 μL) and eluted using 60% (v/v) acetonitrile in 0.1% (v/v) formic acid (100 μL). All SPE steps were carried out using a vacuum manifold. Eluted peptides were evaporated to dryness in a SpeedVac centrifuge and reconstituted in a solution of 1% dimethyl sulfoxide (v/v) in 1% (v/v) formic acid. Reconstituted peptide solutions were quantified on a NanoDrop spectrophotometer using measured signal relative to a standard of known concentration.

Data acquisition (DIA). Data-independent acquisition analysis was carried out on an Orbitrap Lumos instrument interfaced to a Waters M-Class liquid chromatography (LC) system. For analysis, peptide samples were initially trapped on a pre-column (ACQUITY UPLC M-Class Symmetry C18 Trap, 180 µm × 2 cm, 5 µm beads, Waters) and then separated by an analytical column (nanoEase M/Z Peptide BEH C18, 75 µm × 25 cm, 1.7 µm beads, Waters). Samples were initially trapped for 3 min at a flow rate of 20 μL/min of 1% mobile phase B (0.1% formic acid in acetonitrile). Gradient separation of trapped peptides started at an initial condition of 3% mobile phase B, ramped to 7% B in 1 min, to 28% B in 25 min, to 80% B in 0.5 min, hold at 80% B for 2 min, ramp to 3% B in 0.5 min, and a final hold at 3% B for 11 min for a total run time of 40 min (mobile phase A = 0.1% formic acid in water, flow rate = 300 nL/min). For analysis on the Orbitrap Lumos, each peptide sample was injected three individual times (40-min acquisition time for each), with each acquisition covering analysis of a separate mass range (injection 1 = 430–550, 2 = 550–690, 3 = 690–930 m/z). Specifically, the Orbitrap Lumos MS was globally set to use a positive ion spray voltage of 2200 V, an ion transfer tube temperature of 275°C, a default charge state of 3, and an RF Lens setting of 45%. The complete duty cycle of the acquisition method consisted of two MS1 scans and two sets of windowed MS2 scans (order MS1-MS2-MS1-MS2). For the first, low mass range injection, the initial MS1 scan covered a mass range of 425–555 m/z at a resolution of 60,000 with an AGC target of 4e5 (100%) and a max injection time set to ‘Auto’. The following set of DIA MS2 scans covered a precursor range of 430–550 m/z with an isolation window size of 4 m/z (0 m/z overlap) for a total of 30 scan events. Each scan used an HCD energy of 30% and covered a defined mass range of 200–1800 m/z at a resolution of 30,000 with a normalized AGC target of 1000%, and the maximum injection time set to ‘Auto’. Loop control was set to ‘N’ with a value of 30 spectra. The next MS1 scan and following DIA MS2 used the same settings as the previous exception that the MS2 precursor mass range was set to 428–552 m/z to give a 2 m/z stagger with the previous scan windows. For the second, medium mass range injection, the scan settings were the same as for low mass injection with the following changes: (i) MS1 scan range of 545–695 m/z; (ii) DIA-MS2 precursor scan ranges of 550–690 m/z and 548–692 m/z; (iii) Loop control number of spectra = 35. For the third, high mass range injection, the scan settings were the same as the others with the following changes: (i) MS1 scan range of 685–935 m/z; (ii) DIA-MS2 precursor scan ranges of 690–930 m/z and 686–934 m/z; (iii) DIA-MS2 isolation window of 8 m/z; 4. Loop control number of spectra = 30. All scan data were acquired in centroid mode.

Sample preparation for data-dependent acquisition (DDA) analysis. As with DIA, protein was extracted from cell pellets using 150 μL of solution (20 mM Tris-Cl pH 7.5, 150 mM KCl, 5 mM MgCl2, 0.5% (v/v) NP-40, 0.5% (w/v) sodium deoxycholate, and 0.5X cOmplete protease inhibitor). DNA/RNA was supported using 10 μL of a nuclease solution (1 U/μL benzonase) added to the solution and incubated at 24°C for 10 min. After nuclease digestion, 50 μL of a solution of 600 mM Tris Cl pH 7.5, 8% (w/v) SDS, and 10 mM DTT was added to the lysate and incubated at 95°C for 5 min. After cooling to 24°C, lysates were alkylated using 20 mM chloroacetamide (final concentration) and quenched with DTT. Protein concentration of the lysate was measured using a BCA assay (Thermo Scientific). Protein cleanup prior to digestion was performed using a version of the SP3 protocol43. 10 μL of a 50 mg/mL carboxylate-coated bead stock (Cytiva) was added to the protein solution followed by acetone (80% v/v final) to support precipitation44. After a 5-min, 24°C incubation, mixtures were centrifuged for 5 min at 5000 × g. Supernatants were discarded and pelleted material resuspended in 800 μL of 80% ethanol and pelleted again. Ethanol supernatant was removed and 1 μg of trypsin/LysC (Promega, V5071) in 100 μL of 4 M urea was added to each sample and incubated at 30°C for 4 h with shaking (800 rpm). Samples were then diluted to 1 M urea by addition of 300 μL of 100 mM Tris-Cl pH 7.5 and digested for an additional 16 h under the same conditions. Supernatant was then collected after 2 min centrifugation at 12,000 × g and desalted using Strata-X SPE (Phenomenex, 8M-S100-4GA) micro-elution plates. SPE wells were initially conditioned with 100% (v/v) methanol (100 μL) followed by equilibration with 0.1% (v/v) trifluoroacetic acid (TFA, Thermo Scientific) (100 μL). Loaded peptides were rinsed using 4% (v/v) methanol in 0.1% (v/v) formic acid (Thermo Scientific) (200 μL) and eluted using 60% (v/v) methanol in 0.1% (v/v) formic acid (100 μL). SPE steps were all carried out using a centrifuge (250 × g for 3 min). Eluted peptides were evaporated as above for DIA and reconstituted in 1% dimethyl sulfoxide (v/v) in 1% (v/v) formic acid.

Data acquisition (DDA). As with DIA, DDA-MS analysis was performed on an Orbitrap Lumos connected to an UltiMate 3000 LC. Peptide samples were trapped on a pre-column (75 µm × 3 cm, 1.9 µm C18 Reprosil-Pur Basic beads from Dr. Maisch) and separated by an analytical column (75 µm × 25 cm with an integrated nanospray tip, 1.9 µm C18 Reprosil-Pur Basic beads). Data acquisition utilized a standard DDA scheme. The Orbitrap Lumos was globally set to use a positive ion spray voltage of 2200 V, an ion transfer tube temperature of 275°C, a default charge state of 2, advanced precursor determination activated, and an RF Lens setting of 45%. MS1 survey scans covered a mass range of 375–1400 m/z at a resolution of 60,000 with an automatic gain control (AGC) target of 4e5 (100%) and the injection time set to ‘Auto’. Precursors for tandem MS/MS (MS2) analysis were selected using monoisotopic precursor selection (‘Peptide’ mode), charge state filtering (2-4 z), dynamic exclusion (20 ppm low and high, 60 s duration, exclude isotopes = TRUE, exclude within cycle = TRUE), and an intensity threshold of 2.5e4. MS2 scans were carried out in the Orbitrap at a resolution of 30,000 using a higher-energy collision dissociation setting of 30%, an AGC target of 2e5 (400%), injection time set to ‘Auto’, and an isolation window setting of 1.6 m/z. The total allowable MS2 cycle time was set to 3 s. MS1 data were acquired in profile mode, and MS2 in centroid.

Data analysis. MS data were analyzed against a fasta reference genomic/transcriptomic database generated from (i) transcriptome assemblies stemming from the Marine Microbial Eukaryote Transcriptome Sequencing Project45, (ii) JGI-based gene models for A. limacinum, and (iii) ORFs equal or greater than 150 nucleotides generated from AurliV-1 and AurliV-2 consensus sequences from the nanopore-derived, Illumina-polished A. limacinum genome assembly. DDA-MS raw data were analyzed using FragPipe software (version 22.0) with MSFragger (version 4.1) using the default workflow settings46. Specifically, precursor and fragment tolerances of 20 parts-per million (PPM), a fixed modification of 57.02146 at C, and variable modifications of 15.9949 at M and 42.0106 at protein N-termini were used. Matching spectra were rescored using MSBooster and peptide spectral match (PSM) validation utilized Percolator prior to ProteinProphet scoring to support filtering at <1% FDR at the PSM, peptide, and protein levels. Peptides were quantified using the IonQuant package (version 1.10.27) built into FragPipe with default settings. Only protein groups with two or more razor peptides were considered as having been identified by proteomics and used for subsequent analyses.

DIA raw data were de-multiplexed using msconvert (ProteoWizard, version 3.0, peakPicking=‘vendor’, demultiplex=overlaponly) and analyzed using FragPipe software with DIA-NN (version 1.8.2) for quantification. The three injection intensities were summed for each sample. Only protein groups with intensities >0 in at least three of six samples were kept for subsequent analysis. Statistical analysis was performed through the Fragpipe analyst package47, where a cutoff of the p adjusted value of 0.05 and a log2 fold change >1 was used to determine differentially expressed proteins.

Protein annotation and structural predictions. Predicted proteins equal or greater than 150 nucleotides from nanopore assembled, Illumina-polished AurliV-1 and AurliV-2 contigs were annotated using InterPro, BlastP and phmmer (UniRef50). Top hits with the smallest e-values were reported. Further assessments of potential protein function were conducted for select poorly annotated genes using remote homology modeling performed with Phyre2, AlphaFold, FoldSeek, Chai-1, and HHPred23–26,48. The MCP, MCP hexon, triplex and arylsulfatase protein structures were predicted with AlphaFold3 and searched using FoldSeek. VIRALpro was used to confirm prediction of capsid proteins.

Bioinformatic analysis of arylsulfatase and alpha-galactosidase proteins
Arylsulfatase and alpha-galactosidase homologs with similar architectures were identified using NCBI’s Conserved Domain Architecture Retrieval Tool (CDART). Additional homologs were identified using BLASTP, TBLASTN, and RPS-BLAST searches against (i) a custom database representing genomic data from multiple protist lineages49, (ii) a long-read-based assembly and annotation of the Aurantiochytrium limacinum genome (Aurli2)19, NCBI’s (iii) non-redundant protein database, (iv) a metagenomic protein database (env_nr), and (v) transcriptome shotgun assemblies (tsa_nr). RPS-BLAST hits with an e-value <0.005 were retained. Alignments of arylsulfatase and alpha-galactosidase homologs were generated using MUSCLE. Alignments were trimmed with TrimAl. Phylogenetic trees were constructed using IQ-TREE, with the models Q.pfam+R5 (arylsulfatase) or Q.pfam+I+G4 (alpha-galactosidase), and midpoint rooted.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Chung et al present an impressive study that demonstrates virion production in the marine protist Aurantiochytrium limacinum, which has two endogenous mirusviruses in its genome (one as an episome, one as an integrated copy). Previous work from these authors described the endogenous viruses, but it was unclear if they were viable and could initiate a successful infection. Here, they use qPCR, RNA-seq, TEM, and proteomics to show that a small number of viruses are indeed produced in different culture conditions, and they provide evidence for a kind of low-level persistent infection. No other model system for mirusviruses has yet been developed, so this is novel and significant. Overall the experiments are well-done and the results are quite timely. The study of large DNA viruses in protists is surging in recent years, so many readers will be interested in these findings. I have some suggestions to help clarify and expand upon some parts of the manuscript.\n\nMajor comments\n-Some additional details about the hypothesized infection program within the cell would be highly valuable. The asynchronized cultures may seem like a liability, but it may actually be useful for investigating the infection cycle of the mirusvirus. In the giant virus community, infection programs are usually investigated through asynchronous infections so that single TEM experiments can yield different infection stages. In particular, more details about the “peripherally located nucleoli” in Figure 3 would be useful. Is there any evidence of a structure similar to the “virus factory” seen in NCLDVs, or do the TEM results suggest that the entire infection cycle plays out in the nucleus? Based on the TEM it seems like there is something interesting happening with cytoplasmic membranes. This is a quite fascinating part of the infection cycle, so more details would be welcome. A tool was recently developed to look for the scaffolding protein that NCLDVs use to form their virus factory – because these proteins have conserved energetic/chemical properties, it may be possible to see if these mirusviruses here have these proteins too, even if they are divergent from NCLDV homologs (https://www.biorxiv.org/content/10.1101/2024.09.01.610734v2.abstract).\n-Based on the qPCR experiments it seems that some of the sequences (especially those of AurliV-2) are amplified only at very low concentrations. I am concerned that this may be nontemplate amplification and it would be useful to get some information on the noise level in the qPCR results. In Fig 5A, if cycle ~32 is already noise amplification in the assay for the viral genes, the assay should also amplify the luciferase spike-in sequence at a similar cycle when the sequence is not in the sample. It is possible to add in this information? In general I think the qPCR assay is useful for assessing relative abundance shifts, which is the major discussion point of the authors, but eventually the authors may wish to use a plasmid construct with the viral genes so that absolute abundances can be ascertained using a standard curve. This would also allow for the the host/virus ratio to be quantified, which would be highly valuable.\n-Some additional details on the proteomics data analysis would be welcome. Was there a cutoff of the number of peptides matching to a protein used? Many studies only report proteins with either 2 unique peptides matching, or one peptide above a certain score. I was a bit unclear on this in the Methods, and I think it is relevant to the conclusions because the proteomics is probably the most sensitive method used to detect virions. Only 8 genes from AurliV-2 were found in the proteome, which could perhaps be due to false-positive spectral mapping or peptides matching to areas of the protein that are identical between AurliV-1 and 2? I am not opposed to the hypothesis that AurliV-2 is revived and low levels – this would indeed be fascinating. But because detection of this element is a bit at the margins some additional discussion and description of the methods may be useful.\n\nMinor comments:\n-I am curious why the authors chose cold shock as a treatment as opposed to heat shock or some other stressor (e.g. UV or a mutagen). This is not a critique – I am just curious if there is some rationale for this that may bear relevance to the expected natural conditions experienced by the thraustochytrid host.\n-In the abstract and Fig 5 legend the authors refer to some genes detected in the proteome as “up-regulated” but I this term is typically used for mRNA expression levels. Maybe clarify that these are “enriched”?\n-It is not entirely clear what is meant by the phrase: “The putative functions of these proteins suggest that both 5 mirusviral genomic elements have the potential to yield viral particles and greatly enhance our ability to infer mirusvirus-host interactions.” Are the authors indicating that there is no obvious evidence that these elements are relicts or show signs of genomic erosion? It may be useful to state that directly and describe what signatures would indicate erosion (probably duplications, lack of some key marker genes, gene truncations, etc).\n-“the average RNA-seq read count (arbitrary units, AU) for AurIiV-1 was 1218.95 (standard error 180.68), compared to 175.04 (standard error 12.26) for AurIiV-2.” – how does this compare to host gene expression? This would be an interesting comparison. Also, to be clear, not all genes are equal, so it would be useful to specifically discuss the expression of the MCP, triplex, and portal protein genes here, because those would be the signature genes needed for virion formation. Were these genes found for AurliV-2?\n-For the highest expressed gene AurIiV-1 JGI 48845 some annotation information is provided with pyre2. Wouldn’t it be better to use foldseek results from the alphafold predicted structure? In my experience, this yields much better results than pyre2, and it seems that the alphafold results are already available?\n-“The most abundant AurIiV-2 transcript stemmed from JGI annotation Aurli1 142988” – how is this different than the most highly expressed gene discussed in the previous paragraph? I’m not clear on what is being used to quantify gene expression here. For purposes of confirming virion presence, the structural proteins should be far more abundant than the others, so perhaps some discussion of the relative abundance of the proteomic data would be warranted too (especially for AurliV-2)?\n-“Of relevance here is the fact that these experiments were carried out using RNA extracted from unsynchronized bulk cultures” – this is interesting – across a bulk culture there will be some cells that may experience starvation conditions even if the culture as a whole is nutrient replete.\n-“This represents an interesting parallel to the recent discovery of an endogenous giant (NCLDV-type) viral element in lab strains and natural 12 isolates of the green alga Chlamydomonas” – it is also worth noting that in chlamy the NCLDV virions are produced in later-stage growth, which is similar to the starvation conditions described in this study. So there are actually quite a few interesting similarities here.\n-At several points in the manuscript the authors refer to data not shown, but I think it would be easy enough to include supplementary figures that demonstrate these points.\n-“it replicates constitutively in healthy cells” – this isn’t necessarily true, though, right? It could be dormant most of the time and only activated under certain conditions, and it just looks persistent because of population heterogeneity? There could be a lot of hidden complexity at the single-cell level.\nFrank Aylward", "reviewer_2": "In this manuscript, Chung et al., describes persistent infection of two mirusviruses in heterotrophic protists, Aurantiochytrium sp., by following up on the previous publication (https://doi.org/10.1016/j.cub.2023.10.009). Using transcriptome, genomic and proteomic analyses, and electron micrographic imaging, the authors demonstrated persistent and continuous viral expressions and observed viral particles in fed, starved, and cold-shock treatments. Additionally, the authors attempted correlated cold and starved conditions to the re-activation of persistent viruses from its 'low-expressing' stage. The noteworthy finding of their work is; 1) the first time demonstrating the infection strategy of mirusviruses in lab experiments with the only known host so far; 2) potential non-metazoan system for studying persistent and slow infection of dsDNA viruses. I recommend this article to be published upon clarifying some comments here and in the manuscript attached (minor revision?).\n\nDespite the challenges in synchronizing the cultures for differential gene expression experiments, and dynamic life stage of A. limacium, the authors used appropriate method to address what would trigger high infection, as observed in herpesviruses. As I commented in the manuscript (see attachment), some clarification may improve the manuscript should the authors agree - 1) how different life stage of the host may affect persistent infection or reactivation of viral lysis? - furthermore, impacting the interpretation of the data; 2) are there any antagonistic or synergistic gene expression between Aurli-1 and Aurli-2? The result of the second point is something that author never addressed in the manuscript (or I couldn't find other than the volcano plot).\n\nMy takeaway from the experiments regarding cold vs. starved treatment are still unclear - there are differential genes expressed between fed vs. starved (with starved hosts had more virions observed) but for the cold-treatment, only few genes show differences (no imaging due to the reasons described by the authors). What does it mean? Does it mean that the reactivation is not temperature-dependent but nutrient? Based on sequencing data alone, are there evidences or lack their of, that the cold treatment (doesn't) trigger the high virion production? There's no conclusion about this..\n\nLastly, the authors don't highlight enough how unique this finding is, even in the conclusion. Persistent infection is known in herpesviruses and some similar types of infection have been observed in several phycodnaviride. This has a huge evolutionary implications in terms of the viral infection strategy and add to the chimeric characteristic of the mirusviruses. This comment is also made in the main manuscript.", "reviewer_3": "We thank the editor for this opportunity to review the paper “Persistent mirusvirus infection in the marine protist Aurantiochytrium”. The paper presents an intriguing study of a virus infecting a single celled organism of ecological and biotechnological importance. Overall the paper is well written, addresses a clearly defined gap in the literature and presents some interesting conclusions. While we are not experts in the field it appears that it will open new avenues of research. Our expertise is primarily in biological mass spectrometry and consequently our comments will focus on those aspects of the paper. There are some issues with the clarity of the presented methods that make it difficult to assess the quality of the proteomics data presented.\n\nMethods - Sample preparaton:\nStarvation, Cold Shock and Fed Treatments\nCan the authors please provide a table or more clarification on the time points selected for starvation conditions and the time used in the orbital shaker? We do not see these time points presented in the figures (I’m presuming figure 5 would be the natural place to show these).\n\nProteomics/Sample preparation:\n‘nuclease solution’ is unclear: what is its composition and/or source?\n\nData acquisition:\n‘a standard data-dependent acquisition scheme’ is not sufficient detail to reproduce the work: we would expect at least the top-N details, mass range selected from, charge states etc, and this would take up minimal word count.\n\nFrom my understanding, there are three proteomics experiments provided in the paper: a quantitative experiment comparing A. limacinum under rich vs starvation conditions, and two expression cataloguing experiments on nuclear-enriched pellets and supernatants. The first was performed using Data Independent Acquisition and the latter two using Data Dependent Acquisition. This set of experiments is rather hard to parse from the existing methods and results and it would be helpful to provide some sort of flow diagram that summarises the experimental procedures performed. Can the authors also provide a figure (perhaps supplementary) demonstrating the confidence in identification in the cataloguing experiments? E.g. what was the total number of protein groups? How were protein groups annotated? Were a minimum number of peptides required to confidently annotate a protein group?\n\nAdditionally, the experimental plan shows that the starvation experiment was performed in biological triplicate and then each sample was injected three times. How was this data handled statistically? And in figure 5D, are the error bars indicative of biological or technical variation?\n\nFinally, how was normalisation performed on the data? Were the samples normalised by dry weight/protein content/etc before digestion? Or was this done at the data analysis level. If so by what method?\n\nIn all, we believe that these clarifications will increase the impact of the paper and allow other researchers to expand on these interesting findings.", "reviewer_4": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": "We thank Reviewer 1 for their time and constructive feedback on our manuscript. We are pleased that the novelty and impact of these first ‘wet lab’ studies of mirusviruses are recognized.\n\nMajor comments\n-Some additional details about the hypothesized infection program within the cell would be highly valuable. The asynchronized cultures may seem like a liability, but it may actually be useful for investigating the infection cycle of the mirusvirus. In the giant virus community, infection programs are usually investigated through asynchronous infections so that single TEM experiments can yield different infection stages. In particular, more details about the “peripherally located nucleoli” in Figure 3 would be useful. Is there any evidence of a structure similar to the “virus factory” seen in NCLDVs, or do the TEM results suggest that the entire infection cycle plays out in the nucleus? Based on the TEM it seems like there is something interesting happening with cytoplasmic membranes. This is a quite fascinating part of the infection cycle, so more details would be welcome. A tool was recently developed to look for the scaffolding protein that NCLDVs use to form their virus factory – because these proteins have conserved energetic/chemical properties, it may be possible to see if these mirusviruses here have these proteins too, even if they are divergent from NCLDV homologs (https://www.biorxiv.org/content/10.1101/2024.09.01.610734v2.abstract).\nThank you for these insightful comments. It is good to know that our studies of mirusvirus using asynchronous thraustochytrid cultures are comparable to those used by giant virus researchers for identifying different stages of infection in a single experiment. Establishing methods for de novo synchronous infection of naïve thraustochytrid hosts and methods for single-cell analysis are important future goals, which will hopefully bring the mirusvirus infectious cycle into sharper focus. As discussed below, based on the information currently in hand – including new data generated during revision – we believe that in Aurantiochytrium at least, mirusvirus gene expression, genome replication, and assembly all take place in the nucleus.\n\nIn our mixed cell populations, some cells are full of herpesvirus-like capsids, consistent with our genomic analysis, RNA-seq data, and proteomics experiments that indicate accumulation of viral gene products, including proteins related to herpesvirus structural proteins. In cell nuclei, we observe angularized capsids with electron dense centers that resemble herpesvirus ‘type C’ capsids that are the product of proper capsid assembly and DNA packaging [Tandon et al. 2015. Viruses. 7, 899-], which in animals are normally selectively exported to the cytoplasm via a nuclear egress mechanism. We also see what appear to be ‘empty’ capsids in the nucleus that could represent newly assembled capsids prior to DNA packaging, or could be the product of failed DNA packaging, both of which are commonplace in nuclei of herpesvirus-infected cells. Similarly, a hallmark of herpesvirus infected nuclei is margination of nucleoli and chromatin [Rosemarie et al. 2023. mBio. E0268622; Rosemarie and Sugden 2023. PLoS Biol. 21, e3002347], which is something that originally caught our attention as well as Reviewer 1 (our original Fig. 3).\n\nIn our revised manuscript we present an additional panel of TEMs (Fig. S2 of revised manuscript), which shows unusual nucleolus-like structures with demarcations resembling procapsids (Fig. S2b, c, e and f), as well as mature capsids very close to – and possibly in the process of budding from – the nuclear membrane (Fig. S2h-m). The observation of ‘budding’ events into deformed vesicles (e.g., Fig. 3d and e, Fig. S2h, i, l and m) is reminiscent of herpesvirus envelopment by budding at the trans Golgi network and thus may give some clues about how mirusviruses acquire a lipid envelope. It will take substantially more research to fully elucidate the mechanism(s) involved, but it is worth noting that we observe very few apparently empty capsids in the cytoplasm, which is consistent with a selective nuclear egress mechanism for capsids with newly packaged mirusviral DNA.\n\nOverall, our original and new TEM images support the idea that mirusvirus capsid assembly and DNA packaging takes place in the nucleus, as in herpesviruses. We have added more of this informed speculation about mirusvirus assembly and egress mechanisms on page 6 of our revised manuscript, while being cautious not to overinterpret our data based on well characterized herpesvirus mechanisms in the literature.\n\nAs for identification of virus factory ‘scaffolding proteins’ in mirusviruses, we thank Reviewer 1 for bringing the Rigou et al. preprint to our attention. This is something well worth exploring soon. However, given the anticipated complexity of the search process for such proteins (machine learning, etc.) and the time it has taken for us to address the most substantive concerns of the reviewers, it is not possible for us to include the results of such an investigation in the present manuscript.\n\n-Based on the qPCR experiments it seems that some of the sequences (especially those of AurliV-2) are amplified only at very low concentrations. I am concerned that this may be nontemplate amplification and it would be useful to get some information on the noise level in the qPCR results. In Fig 5A, if cycle ~32 is already noise amplification in the assay for the viral genes, the assay should also amplify the luciferase spike-in sequence at a similar cycle when the sequence is not in the sample. It is possible to add in this information? In general I think the qPCR assay is useful for assessing relative abundance shifts, which is the major discussion point of the authors, but eventually the authors may wish to use a plasmid construct with the viral genes so that absolute abundances can be ascertained using a standard curve. This would also allow for the the host/virus ratio to be quantified, which would be highly valuable.\nWe appreciate Reviewer 1’s comments about the qPCR results included in our original submission. We shared those concerns and considered the suggestion to use a plasmid construct to generate standard curves for absolute quantification. In the end we decided to go a big step further and perform ultrasensitive and quantitative droplet digital PCR (ddPCR) analysis of capsid-protected AurliV-1 and Aurli-V2 genomes in cell supernatants. These new results, along with luciferase controls, are presented in Figure 5a and b and Figure S6 of our revised manuscript. The original qPCR data have now been moved entirely to Supplementary Figure S6, as we believe it still provides some value in helping us piece together the host-virus dynamics in this complex system. Note too that we have removed our RT-PCR data; these experiments were exploratory in nature and not sufficiently replicated to allow us to come to any meaningful conclusions (see below).\n\nOur new ddPCR analysis is generally consistent with our previous data acquired via PCR and qPCR methods and has the advantage of being able to give us an absolute number of viral genomes per unit volume of cell supernatant, enabling robust quantitative comparisons. As described on pages 7 and 8 of our revised manuscript, our ddPCR experiments (performed in triplicate) demonstrate conclusively that both AurliV-1 and Aurli-V2 genomes are released into the extracellular environment, presumably encased in viral capsids. These results are consistent with our original qPCR data in showing that Aurli-V1 is much more abundant than AurliV-2, but that the latter is also present.\n\nWe have experience in the capture of large viral genomes for study and mutagenesis [Duguay, et al. J. Virol. 2025. E01671-24] and have plans to capture and assemble the AurliV-1 and AurliV-2 genomes independently using transformation-associated recombination methods in yeast. These genomes can be used for applications like Reviewer 1 suggests, as well as for transfection into naïve hosts to study viral genome retention and whether the AurliV-1 episome is sufficient to launch an infectious cycle. This kind of tool building will take time but will hopefully allow us to dissect relationships between the AurliV-1 and AurliV-2 viruses and their protist host.\n\n-Some additional details on the proteomics data analysis would be welcome. Was there a cutoff of the number of peptides matching to a protein used? Many studies only report proteins with either 2 unique peptides matching, or one peptide above a certain score. I was a bit unclear on this in the Methods, and I think it is relevant to the conclusions because the proteomics is probably the most sensitive method used to detect virions. Only 8 genes from AurliV-2 were found in the proteome, which could perhaps be due to false-positive spectral mapping or peptides matching to areas of the protein that are identical between AurliV-1 and 2? I am not opposed to the hypothesis that AurliV-2 is revived and low levels – this would indeed be fascinating. But because detection of this element is a bit at the margins some additional discussion and description of the methods may be useful.\nWe thank Reviewer 1 for raising this issue, which was also brought up by Reviewer 3 (see below). As requested, we have greatly expanded the text devoted to describing our proteomics analysis (both in the main text and methods section). The concern that the low number of AurliV-2 proteins detected in our proteomics data could be artifactual is one that we had initially considered. We can say that it is not due to false-positive spectral mapping and/or the possibility of peptides mapping to regions of proteins encoded by, and conserved in, both the AurliV-1 and AurliV-2 genomes. It is important to note that we have reanalyzed our proteomics data using our new gap-free AurliV-1 and AurliV-2 genomes, newly produced on a PromethION instrument using the latest Oxford Nanopore chemistry. This allowed us to detect an additional AurliV-2 protein (bringing the total to 9). Recall that the AurliV-1 and AurliV-2 genomes are highly divergent from one another and their protein-coding genes sets are very different [only 48 of 177 originally predicted Aurli-V1 proteins have a reciprocal best BLAST match to a Aurli-V2 homolog; Collier et al. 2023, Curr. Biol. 33, 5199-]. And to Reviewer 1’s specific concern, there is zero overlap between the 9 AurliV-2 proteins detected using proteomics (none of which correspond to viral ‘marker genes’) and the 67 detected AurliV-1 proteins (one more than originally reported). All this is to say that we are confident that a small but non-zero number of proteins encoded by AurliV-2 genes can be detected by proteomics. The biological significance of this observation is an important question for the future.\n\nMinor comments:\n-I am curious why the authors chose cold shock as a treatment as opposed to heat shock or some other stressor (e.g. UV or a mutagen). This is not a critique – I am just curious if there is some rationale for this that may bear relevance to the expected natural conditions experienced by the thraustochytrid host.\nAurantiochytrium is a tropical organism, and we have observed that it does not do well under the common ‘throw the plates in the fridge’ microbial storage conditions. When considering ways to stress Aurantiochytrium, cold shock thus seemed like a good place to start. Similarly, zoospore formation is known to be associated with starvation, and in fact Kazama and Schornstein (1972/73) suggested that the “herpesvirus” they observed in a thraustochytrid (then considered a fungus) could be specifically activated in starved zoospores (we make this point clearer in our revised manuscript). Starvation was thus chosen as a treatment as well.\n\nIn the future, we plan to study different stressors and their effects on both Aurantiochytrium mirusviruses, with the goal of elucidating relationships between host stress response pathways and viral replication. But for now, please note that we have decided to remove mention of our cold-shock treatment in the revised manuscript. The reason is that it was by far the weakest piece of the puzzle – cold-shocked Aurantiochytrium cells were obviously sick, to the point that they did not yield interpretable data in TEM experiments, and we could not proceed with them for qPCR. They also gave inconsistent results in our RT-PCR experiments (Figure S2 of our original manuscript). We felt it was best to present only data from control (i.e., ‘fed’) and starved cells.\n\n-In the abstract and Fig 5 legend the authors refer to some genes detected in the proteome as “up-regulated” but I this term is typically used for mRNA expression levels. Maybe clarify that these are “enriched”?\nAs suggested, we have changed “up-regulated” to “enriched” in the abstract and all other relevant locations. We have also changed “down-regulated” to “depleted” on page 9.\n\n-It is not entirely clear what is meant by the phrase: “The putative functions of these proteins suggest that both 5 mirusviral genomic elements have the potential to yield viral particles and greatly enhance our ability to infer mirusvirus-host interactions.” Are the authors indicating that there is no obvious evidence that these elements are relicts or show signs of genomic erosion? It may be useful to state that directly and describe what signatures would indicate erosion (probably duplications, lack of some key marker genes, gene truncations, etc).\nThere are indeed no obvious signs of genome erosion in the AurliV-1 and V-2 genomes, and we are even more confident in saying this having produced high-quality, gap-free consensus sequences for both genomes using the latest ONT chemistry. We have modified the sentence as follows:\n“Combined with the fact that neither AurliV-1 or AurliV-2 show obvious signs of genome erosion19, the putative viral functions of many of these proteins suggest that both mirusviral genomic elements represent current and possibly active infections of their hosts and could yield viral particles.”\n\n-“the average RNA-seq read count (arbitrary units, AU) for AurIiV-1 was 1218.95 (standard error 180.68), compared to 175.04 (standard error 12.26) for AurIiV-2.” – how does this compare to host gene expression? This would be an interesting comparison. Also, to be clear, not all genes are equal, so it would be useful to specifically discuss the expression of the MCP, triplex, and portal protein genes here, because those would be the signature genes needed for virion formation. Were these genes found for AurliV-2?\nAs requested, we have included RNA-seq read counts for the host to better contextualize the observed levels of AurliV-1 and V-2 transcription (recognizing that based on sequence data, AurliV-1 copy number is ~9X higher than the nuclear chromosomes and the integrated AurliV-2, which somewhat complicates interpretation). The new data, presented in Figure S1a as a whisker plot to better illustrate data dispersion, show that AurliV-2 gene expression is significantly different than the rest of the host nuclear genome, whereas the episomal AurliV-1 is not.\n\nRegarding the presence and expression of MCP, triplex, and portal protein genes, AurliV-1 and AurliV-2 both contain homologs of these signature genes (see Table 1 in Collier et al. 2023 Current Biology). Their RNAseq expression varies from ~70 to ~2000 A.U. for AurliV-1 and all are detected in the AurliV-1 proteomic data (See Table S1; MCP (Aurli1_48830) appears in Fig. S1b as the 12th most highly expressed AurliV-1 gene by RNA-seq). Interestingly, the AurliV-2 versions of these genes exhibit ~100-fold lower expression in RNA-seq (MCP and portal protein transcripts were not detected) and were not detected by proteomics (Figure 6 of revised manuscript). Naturally our manuscript focused more on the proteomics results when it comes to the MCP, triplex and portal proteins, but more information on this topic is now provided in the RNA-seq section of the revised manuscript (page 5/6) as a nod forward to the data that follow.\n\n-For the highest expressed gene AurIiV-1 JGI 48845 some annotation information is provided with pyre2. Wouldn’t it be better to use foldseek results from the alphafold predicted structure? In my experience, this yields much better results than pyre2, and it seems that the alphafold results are already available?\nWe have now modeled all AurliV-1 and V-2 predicted proteins using Chai-1 and, in the case of proteins detected using proteomics, Alphafold3 / FoldSeq. In select cases, such as predicted proteins stemming from the three genes with the highest level of RNA-seq, we have also used Phyre2 and HHPred. These details are now discussed on pages 5 and 6 of the revised manuscript. The take-home message is that structure-based predictions somewhat improved our ability to assign putative functions to the most highly expressed mirusvirus genes.\n\n-“The most abundant AurIiV-2 transcript stemmed from JGI annotation Aurli1 142988” – how is this different than the most highly expressed gene discussed in the previous paragraph? I’m not clear on what is being used to quantify gene expression here. For purposes of confirming virion presence, the structural proteins should be far more abundant than the others, so perhaps some discussion of the relative abundance of the proteomic data would be warranted too (especially for AurliV-2)?\nWe apologize for the confusion over this statement. Here we are referring to an AurliV-2 transcript, whereas the previous paragraph is devoted to AurliV-1 transcripts. We have rectified this in the revised text on pages 5 and 6. And as noted above, we have specifically discussed the expression of the hallmark / structural proteins on page 5.\n\n-“Of relevance here is the fact that these experiments were carried out using RNA extracted from unsynchronized bulk cultures” – this is interesting – across a bulk culture there will be some cells that may experience starvation conditions even if the culture as a whole is nutrient replete.\nThe link to starvation is indeed interesting, and of course there are other features of unsynchronized cultures that probably contribute to differences in expression / packaging / release of AurliV-1 and/or V-2, such as lifecycle stages (some of which themselves can be linked to starvation). As emphasized elsewhere in this document and in the manuscript itself, we simply do not have enough information to make robust inferences about stress-induced viral replication. Combined with continued investigation of bulk culture, single-cell approaches to this question will be very important moving forward.\n\n-“This represents an interesting parallel to the recent discovery of an endogenous giant (NCLDV-type) viral element in lab strains and natural 12 isolates of the green alga Chlamydomonas” – it is also worth noting that in chlamy the NCLDV virions are produced in later-stage growth, which is similar to the starvation conditions described in this study. So there are actually quite a few interesting similarities here.\nThe parallels between what we see in Aurantiochytrium and what Erazo-Garcia et al. described in Chlamydomonas is indeed remarkable. In addition to discussing and citing that paper on page 10, we have added mention of ‘late-stage growth’ in our revised manuscript.\n\n-At several points in the manuscript the authors refer to data not shown, but I think it would be easy enough to include supplementary figures that demonstrate these points.\nWe appreciate this comment and address it in several ways. In the case of the results of our genome-wide scan for 5-methylcytosine, which reveals no obvious epigenetic explanation for the observed difference in expression between AurliV-1 and AurliV-2, we now present the data as Figure S7. The only other instances of ‘data not shown’ in our original manuscript was mention of our nanopore re-sequencing of the genome (now formally described in our revised manuscript) and reference to the lack of quality SEM and TEM data stemming from cold-shocked samples. These images are of very poor quality and, as noted above, we have removed all of data associated with our preliminary investigation of cold shock. There are now no mentions of ‘data not shown’.\n\n-“it replicates constitutively in healthy cells” – this isn’t necessarily true, though, right? It could be dormant most of the time and only activated under certain conditions, and it just looks persistent because of population heterogeneity? There could be a lot of hidden complexity at the single-cell level.\nPoint well-taken. In this instance, we have removed “constitutively” and changed “healthy” to “cells grown under standard growth conditions” (Page 13). The reviewer is correct that there is likely hidden complexity at the single-cell level, and we have done our best to emphasize that point.\n\nFrank Aylward\nThank you Frank for your very helpful critique of our work.\n\nReviewer #2 (Remarks to the Author):\nIn this manuscript, Chung et al., describes persistent infection of two mirusviruses in heterotrophic protists, Aurantiochytrium sp., by following up on the previous publication (https://doi.org/10.1016/j.cub.2023.10.009). Using transcriptome, genomic and proteomic analyses, and electron micrographic imaging, the authors demonstrated persistent and continuous viral expressions and observed viral particles in fed, starved, and cold-shock treatments. Additionally, the authors attempted correlated cold and starved conditions to the re-activation of persistent viruses from its 'low-expressing' stage. The noteworthy finding of their work is; 1) the first time demonstrating the infection strategy of mirusviruses in lab experiments with the only known host so far; 2) potential non-metazoan system for studying persistent and slow infection of dsDNA viruses. I recommend this article to be published upon clarifying some comments here and in the manuscript attached (minor revision?).\nWe appreciate Reviewer 2’s summary of the significance of our manuscript and efforts to provide constructive feedback (including provision of a marked-up pdf file with specific comments). We have addressed all these comments to the best of our ability.\n\nDespite the challenges in synchronizing the cultures for differential gene expression experiments, and dynamic life stage of A. limacium, the authors used appropriate method to address what would trigger high infection, as observed in herpesviruses. As I commented in the manuscript (see attachment), some clarification may improve the manuscript should the authors agree - 1) how different life stage of the host may affect persistent infection or reactivation of viral lysis? - furthermore, impacting the interpretation of the data; 2) are there any antagonistic or synergistic gene expression between Aurli-1 and Aurli-2? The result of the second point is something that author never addressed in the manuscript (or I couldn't find other than the volcano plot).\nWe appreciate the thoughtful questions raised by Reviewer 2. We address them in turn below.\n\n1) how different life stage of the host may affect persistent infection or reactivation of viral lysis? - furthermore, impacting the interpretation of the data;\nExcellent question, to which we do not have a definitive answer. As noted above in response to Reviewer 1, the data of Kazama and Schornstein (1972/73) suggested viral activation in starved zoospores (which is one of the reasons why starvation was one of our target treatments). But we are not confident moving beyond speculation at this point, as all our work was done on bulk (non-synchronized) cultures. This limitation is emphasized on page 11 of the revised manuscript. Single-cell approaches are sorely needed here and are a focus of ongoing research in our labs.\n\n2) are there any antagonistic or synergistic gene expression between Aurli-1 and Aurli-2? The result of the second point is something that author never addressed in the manuscript (or I couldn't find other than the volcano plot).\nAgain, a great question but one without an answer. In analyzing RNA-seq data, we detected no obvious antagonistic / synergistic expression patterns. The RNA-seq profiles of AurliV-1 and V-2 are dramatically different in overall transcript abundance and the specific genes they express. Given this complexity and given that the two mirusviral genomes do not share very many genes in the first place, it is not possible to assess competition / synergy between the two. For such correlations to be meaningful in terms of host-virus interactions, they would need to be defined down to the level of single cell types / lifecycle stages, and we are not yet able to do that. Strictly speaking, we do not know whether mosaic viral particles exist – e.g., whether the AurliV-2 genome can be packaged into an AurliV-1-derived viral particle, or whether virions containing proteins encoded by both AurliV-1 and AurliV-2 are possible. If so, our data suggest that such mosaicism would be minimal, as RNA-seq and proteomic data indicate that AurliV-2 genes and proteins are expressed at very low levels relative to AurliV-1 (under the conditions we have examined at least).\n\nMy takeaway from the experiments regarding cold vs. starved treatment are still unclear - there are differential genes expressed between fed vs. starved (with starved hosts had more virions observed) but for the cold-treatment, only few genes show differences (no imaging due to the reasons described by the authors). What does it mean? Does it mean that the reactivation is not temperature-dependent but nutrient? Based on sequencing data alone, are there evidences or lack their of, that the cold treatment (doesn't) trigger the high virion production? There's no conclusion about this..\nAs noted in response to Reviewer #1, we elected to remove our cold-treatment data in the revised manuscript. Cold-shocked Aurantiochytrium cells were not conducive to proper TEM experiments and yielded inconsistent results in our RT-PCR experiments, which is why we quickly moved on to qPCR and, during the revision of this manuscript, droplet digital PCR.\n\nIn any case, with this comment, we take “reactivation” to mean transition from latency to productive lytic replication (see next comment). We do not know the precise links between nutritional status (starvation versus fed), temperature, and induction of AurliV-1 and V-2. As noted, this uncertainty is confounded by the fact that we analyzed bulk (asynchronous) cultures, but also by the nature of our now removed data on cold-shocked cells (limited to a handful of genes by RT-PCR, which showed variation from gene to gene and under different time points). Our cold shock was in fact more of an extended treatment, which may have triggered a variety of cellular responses, only some of which may be linked to viral production. We could not use sequencing to explore a possible link between virion production and cold treatment and, as noted, such cells are very ‘sick’ and do not lend themselves to scale up for TEM or bulk sequencing.\n\nMoving forward, we are particularly keen to learn about the fate of the cells that are chock-full of viral particles (found primarily in cultures grown under starvation conditions). Are they destined to die and release particles to the medium? Is viral infection non-cytolytic, and release occurs upon cell wall degradation during the natural growth/differentiation of Aurantiochytrium? Again, single cell analysis will be key here.\n\nLastly, the authors don't highlight enough how unique this finding is, even in the conclusion. Persistent infection is known in herpesviruses and some similar types of infection have been observed in several phycodnaviride. This has a huge evolutionary implications in terms of the viral infection strategy and add to the chimeric characteristic of the mirusviruses. This comment is also made in the main manuscript.\nWe are grateful for the reviewer’s thoughts and enthusiasm here. We are happy to have the chance to provide a more fulsome discussion of what we have described as persistent infection — what it is, similarities to herpesviruses, episomal maintenance, etc. — and have done so on page 13 of our revised manuscript. We do, however, want to note our desire to not overinterpret the similarities between mirusvirus-protist interactions on one hand, and herpesviruses and their interactions with animal hosts on the other. The mirusvirus - herpesvirus connection is very intriguing, but to our knowledge the term ‘persistent infection’ has not been applied to the biology of herpesviruses and is not, strictly speaking, synonymous with latency. We make this point clear in our revised manuscript.\n\nResponses to Reviewer #2 comments provided directly on manuscript PDF\n[note — page numbers refer to the location of Reviewer 2 comments on the original manuscript]\nPage 3 comments\nThis [demonstration of stable, long-term association between Aurantiochytrium and mirusvirus] is the major finding of this paper but the delivery is not powerful enough. This needs to be linked the typical persistent infection of herpesvirales - by defining what persistent infection is and how it's usually observed (e.g., stay dormant for years without typical lytic events, low to zero transcription activities while dormant, episome maintenance ....). This is not something that's usually observed in NCLDVs but some phycodnavridiae-infecting NCLDVs show this \"stable association\" (e.g., OtVs and OmVs Yau et al, 2016; 2020, PgV Stevens et al 2019; and other Mamielles-infecting viruses)\nFurther to above, we are delighted that Reviewer 2 sees the importance of framing our results through the lens of persistent infection. Please note that at this stage of the text (i.e., the first paragraph of the introduction), our goal was to set the stage for what is to come. As summarized above, we have expanded our discussion of these issues on pages 11-13 of our revised manuscript.\n\n“However, because of their MAG-based discovery, the physical nature of mirusviruses and the range of hosts they infect is unclear.” Agreed but we know that they likely infect euks that are less than ~5um\nWe have added a short parenthetical on page 3 to acknowledge the apparent small cell size of the hosts based on the size fractions where mirusvirus MAGs were detected in the Gaia et al. 2023 paper. We do note, however, that thraustochytrid vegetative cells (including in Aurantiochytrium analyzed herein) are generally larger than 5 micrometers, and the same applies to many of the other putative mirusvirus hosts described by Zhao et al. 2024, Current Biology), including Paulinella, Acanthamoeba, and cryptophye algae.\n\nPage 4 comments\nWhat about the integrant (Aurli-2?)\nAt the tail end of the Introduction, we prefer not to ‘unpack’ the complexity of AurliV-2 because we do not have a definitive answer. We address these uncertainties in the results and discussion sections.\n\nSpell out when the number starts the sentence\nWe thank the reviewer for pointing this out. We have fixed this on page 4 and elsewhere throughout the manuscript.\n\nPage 5 comments\n“did you check GOEV database from the mirusvirus paper to possibly link whether other mirusviruses have similar persistent infection strategy”\nAs presented, we do not understand the meaning of this question in the context of our description of the most highly expressed mirusviral genes in Aurantiochytrium. Is the idea that we could use metatranscriptomics to assess expression of particular AurliV-1 and V-2 genes in the data of Gaia et al.? If so, it’s hard to know how we could assess “persistent infection” in the absence of cultured hosts. Unfortunately, without more information to go on we are unable to address this question.\n\nPage 6 comments\n“this [AurliV-1 JGI = putative kinase] might be relevant to episomal maintenance and replication”\nSuch proteins could indeed play a role in episomal maintenance and replication, though there are many possible paths through which that role could be manifested. Note that we expanded our structural modeling analyses by applying Chai 1 analysis on all predicted V-1 and V-2 proteins and, in the case of proteins detected by proteomics, Alphafold3 / FoldSeq. Using these data, we have expanded our discussion of the most highly expressed Aurli-V-1 and V-2 genes and their potential functions, including putative kinases. Please refer to pages 5 and 6 of the revised manuscript.\n\n“…in cells subject to starvation”. Different nutrient conditions affect life cycle stage of the host. This may impact viral latency or transcriptome expression level more than the stress arising from the starvation.\nThe authors mentioned below that there were infrequent observation of zoospores during starvation but I'm curious what life stages the host were predominant during different treatment and its duration (e.g., rich media mostly vegetative sporangia but no zoospores vs. more zoospores in more starved host?) and whether they relates to patterns in expression or latency - this can be elaborated in Discussion.\nThis is a great question, but unfortunately it is not possible to provide a definitive answer. What we can say is that in the bulk cultures we examined in this study, zoospores were only observed under starvation conditions and even then, only occasionally. It is thus not possible to quantify the proportion of zoospores relative to vegetative cells. As noted elsewhere, the development of single-cell approaches to the study of thraustochytrids and their viruses will hopefully provide more clarity.\n\nPage 7 comment\n“What are they? Didn’t find it in the sections below unless I am missing”\nThe AurliV-1 and AurliV-2 orthologs in question are those examined by RT-PCR and presented in our original Figure S2. These genes were listed in the previous sentence (“…replication genes as well as structural protein genes (RNA polymerase, helicase, DNA polymerase, MCP and terminase). As described above, we have removed our RT-PCR data and thus the relevant text no longer appears in our revised manuscript.\n\nPage 8 comments\n“This is cool, this fits the definition of slow persistent infection (not latency).”\n“This sort of rediscovery is always exciting”.\nWe agree on both points!\n\nPage 10 comment\n“We also compared the proteomes of A. limacinum cultures grown in triplicate under control (fed) conditions … to those of cultures subjected to starvation …” Apart from the transcriptome expression between different treatments, did you see any differences between Aurli-1 and 2, regardless of the treatment? I am curious if there's any competitive/synergistic infection going between the two.\nAs noted above, RNA-seq analysis did not reveal any obvious antagonistic / synergistic expression patterns; given that the two mirusviral genomes have relatively few genes in common and that relatively few mirusviral genes are properly annotated, it is hard to know how we could get a reliable foothold on this interesting problem. What we can say is that between treatments, the gene expression profiles of AurliV-1 and V-2 are dramatically different in overall transcript abundance and in terms of the specific genes they express. With the data in hand, it is not possible to come to any meaningful conclusions about competition / synergy between the two.\n\nPage 12 comment\n“At present, the diversity of A. limacinum cell types in culture is a barrier to further investigation of mirusvirus genome biology at the population level, …” At what life stage did you observe the VLP the most?\nVirus-like particles were observed primarily in solitary vegetative cells (Figures 2 and 3) and occasionally in multicelled sporangia (on this point we have added new data in Figure S2 added at revision, see S2j-m). We have not yet directly observed VLPs in zoospores.\n\nPage 14 comments\n“This [discussion of other DNA and RNA viruses in thraustochytrids] feels out of place”\n“this [text related to SmDNAV and RNA virus] sounds more like random repetitions of a bit of intro.. can you summarize the experimental results here and what that signify?”\nWe apologize for this ill-formed section of text. In our revised manuscript we have revised and expanded it to make clearer the significance of a lytic large DNA virus (SmDNAV) in at least one thraustochytrid, as well as an RNA virus in a Schizochytrium sp. strain. The revised text is as follows:\n“In that vein, it is important to note that although not identified in the A. limacinum strain examined herein, lytic single-stranded RNA viruses have been shown to infect the thraustochytrid Schizochytrium sp. strain NIBH N1-2736,37 and a lytic large DNA virus (SmDNAV) infects Sicyoidochytrium minutum38. We thus cannot discount the possibility that multiple distinct DNA and RNA viruses co-occur in at least some thraustochytrid species and genera, which would make dissecting mirusvirus-host interactions even more difficult.” (Page 13 of revised manuscript)\n\nProteomics section comment\nIn our expanded Materials and Methods section, we have defined all proteomics-related abbreviations.\n\nFigure S5 comment\nWe apologize for the difficulty reading the taxon labels in this tree. The version rendered into the reviewer PDF file defaulted to a low-res image. We have uploaded high-quality images for Figures S8 and S9, fingers crossed these can be accessed appropriately by the reviewer.\n\nReviewer #3 (Remarks to the Author)\nWe thank the editor for this opportunity to review the paper “Persistent mirusvirus infection in the marine protist Aurantiochytrium”. The paper presents an intriguing study of a virus infecting a single celled organism of ecological and biotechnological importance. Overall the paper is well written, addresses a clearly defined gap in the literature and presents some interesting conclusions. While we are not experts in the field it appears that it will open new avenues of research.\nWe thank Reviewer 3 for their thorough and useful feedback on our manuscript. We are pleased that the quality, novelty and impact of these first ‘wet lab’ studies of mirusviruses are appreciated by the reviewer.\n\nOur expertise is primarily in biological mass spectrometry and consequently our comments will focus on those aspects of the paper. There are some issues with the clarity of the presented methods that make it difficult to assess the quality of the proteomics data presented.\n\nMethods - Sample preparaton:\nStarvation, Cold Shock and Fed Treatments\nCan the authors please provide a table or more clarification on the time points selected for starvation conditions and the time used in the orbital shaker? We do not see these time points presented in the figures (I’m presuming figure 5 would be the natural place to show these).\nAs noted above in response to Reviewer 1, we have greatly expanded our description of the proteomics methods to make it clear how our samples were prepared and how MS data were acquired and analyzed for both DIA and DDA (Please refer to pages 19-25 of revised manuscript).\n\nRegarding culture treatments, we have added a detailed list of the analyzed timepoints associated with our different treatments and analyses (SEM, TEM, qPCR, and ddPCR) on page 15 of the revised manuscript. Note that we have removed all our RT-PCR and cold-shock data from the paper (see comments to Reviewers #1 and 2 above).\n\nProteomics/Sample preparation:\n‘nuclease solution’ is unclear: what is its composition and/or source?\nAs described in our expanded methods, the nuclease was benzonase at 1U /μL.\n\nData acquisition:\n‘a standard data-dependent acquisition scheme’ is not sufficient detail to reproduce the work: we would expect at least the top-N details, mass range selected from, charge states etc, and this would take up minimal word count.\nWe appreciate the reviewer’s concerns here, and as noted have expanded our proteomics methods section accordingly (Pages 19-25). We hope that this information is now sufficient to allow our experiments to be reproduced.\n\nFrom my understanding, there are three proteomics experiments provided in the paper: a quantitative experiment comparing A. limacinum under rich vs starvation conditions, and two expression cataloguing experiments on nuclear-enriched pellets and supernatants. The first was performed using Data Independent Acquisition and the latter two using Data Dependent Acquisition. This set of experiments is rather hard to parse from the existing methods and results and it would be helpful to provide some sort of flow diagram that summarises the experimental procedures performed. Can the authors also provide a figure (perhaps supplementary) demonstrating the confidence in identification in the cataloguing experiments? E.g. what was the total number of protein groups? How were protein groups annotated? Were a minimum number of peptides required to confidently annotate a protein group?\nWe appreciate the challenges associated with parsing the materials and methods associated with our different proteomics experiments. We have greatly expanded the materials and methods and relevant sections of the main text to address the issues as best we can, and hopefully conveyed our workflow in a more linear manner (this is helped by addition of several sub-sections of the proteomics methods) Our hope is that a summary flow diagram is no longer needed.\n\nAdditionally, the experimental plan shows that the starvation experiment was performed in biological triplicate and then each sample was injected three times. How was this data handled statistically? And in figure 5D, are the error bars indicative of biological or technical variation?\nWe apologize for the lack of clarity on these points in our original submission. Our expanded materials and methods and revised legend to Figure 6 provide information about statistical methods and significance. All three injections were summed for each sample. The error bars are indicative of biological variation in the Figure 6D.\n\nFinally, how was normalisation performed on the data? Were the samples normalised by dry weight/protein content/etc before digestion? Or was this done at the data analysis level. If so by what method?\nOur methods now clarify this matter. Protein concentrations of the lysates were measured using a BCA assay.\n\nIn all, we believe that these clarifications will increase the impact of the paper and allow other researchers to expand on these interesting findings.\nWe sincerely thank Reviewer 3 for these helpful suggestions to improve the readability and reproducibility of our study.\n\nReviewer #4 (Remarks to the Author)\nI co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts.\nThank you for making this important contribution to training early-career researchers."}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have substantially revised the manuscript to address my comments, and they now include digital PCR quantification of viral particles, new TEMs, and additional text and analysis of proteomic data. It is quite challenging to quantify the dynamics of these latent viruses when their activity is so subtle, so these experiments are nontrivial. The manuscript was already quite impressive and described important findings, and these changes have improved it further. It is ready for publication and will make an important mark on the field - congrats to all authors.\nFrank Aylward", "reviewer_2": "I have no further suggestions for changes. The authors removed the data for the cold-shock treatment which previously complicated the story and further clarified the results from starved vs. enriched media. I particularly appreciate their further elaboration on proteomic data (per Rev #1 comments) and additional images showing how they the virions may replicate (i.e., the budding virus) - as well as addressing some of the questions that required speculations.", "reviewer_3": "The authors have made an excellent attempt to address our concerns in the initial review. The methods are now highly comprehensive and link well to the revised figure and methods for the proteomics so that the experimental plan can be followed effectively and supports the data well.\n\nA point of clarification regarding the selected 2 peptides (row 732) required for identification of a protein - are these selected on the basis of unique peptides or razor peptides? This should be explained in the methods (and ideally should be unique peptides).\n\nSecondly, (in row 739) there is now a fairly terse description of what I think is a minimum detections filter: \"protein groups with intensities <3 but >0 among all six samples were removed from subsequent analysis.\". Do the authors mean that the protein groups with intensities greater than the limit of detection found in 2 or fewer samples out of 6 were removed from the analysis? Or is it, as literally written, peptides with a raw intensity of <3 and >0 were removed? I would think the latter is unlikely on a Lumos instrument. Assuming it's the former, it just needs to be described more clearly.", "reviewer_4": "I co-reviewed this manuscript with one of the reviewers who provided the listed reports. This is part of the Nature Communications initiative to facilitate training in peer review and to provide appropriate recognition for Early Career Researchers who co-review manuscripts."}, "author_response": "We thank Reviewer 3 for their time and constructive feedback on our original manuscript, and these final comments related to our proteomics methodologies.\n\nRegarding the text on our original line 732, these were in fact razer peptides. This is because we used multiple gene models from genomic and transcriptomic data to generate our in silico proteome so as to maximize protein identification. The text has been revised as follows:\n“Only protein groups with two or more razer peptides were considered as having been identified by proteomics and used for subsequent analyses.”\n\nSecondly, (in row 739) there is now a fairly terse description of what I think is a minimum detections filter: \"protein groups with intensities <3 but >0 among all six samples were removed from subsequent analysis.\". Do the authors mean that the protein groups with intensities greater than the limit of detection found in 2 or fewer samples out of 6 were removed from the analysis? Or is it, as literally written, peptides with a raw intensity of <3 and >0 were removed? I would think the latter is unlikely on a Lumos instrument. Assuming it's the former, it just needs to be described more clearly.\nThe reviewer is correct, there was indeed an error in our description of our minimum detection filter. We apologize for the confusion. Only protein groups with intensities >0 in at least three of the six samples collected were retained for downstream analysis. We have corrected the relevant text in the final revised manuscript as highlighted in red below.\n“DIA raw data were de-multiplexed using msconvert (ProteoWizard, version 3.0, peakPicking = ‘vendor’, demultiplex = overlap only) and analyzed using FragPipe software with DIA-NN (version 1.8.2) for quantification. The three injection intensities were summed for each sample. Only protein groups with intensities >0 in at least three of six samples were kept for subsequent analysis. Statistical analysis was performed through the Fragpipe analyst package47, where a cutoff of the p adjusted value of 0.05 and a log2 fold change greater than 1 was used to determine differentially expressed proteins.” Page 23."}]

Decision letter