1. Genetics and Genomics
  2. Microbiology and Infectious Disease
Download icon

Virophages and retrotransposons colonize the genomes of a heterotrophic flagellate

  1. Thomas Hackl
  2. Sarah Duponchel
  3. Karina Barenhoff
  4. Alexa Weinmann
  5. Matthias G Fischer  Is a corresponding author
  1. Max Planck Institute for Medical Research, Department of Biomolecular Mechanisms, Germany
Research Article
  • Cited 0
  • Views 1,155
  • Annotations
Cite this article as: eLife 2021;10:e72674 doi: 10.7554/eLife.72674

Abstract

Virophages can parasitize giant DNA viruses and may provide adaptive anti-giant virus defense in unicellular eukaryotes. Under laboratory conditions, the virophage mavirus integrates into the nuclear genome of the marine flagellate Cafeteria burkhardae and reactivates upon superinfection with the giant virus CroV. In natural systems, however, the prevalence and diversity of host-virophage associations has not been systematically explored. Here, we report dozens of integrated virophages in four globally sampled C. burkhardae strains that constitute up to 2% of their host genomes. These endogenous mavirus-like elements (EMALEs) separated into eight types based on GC-content, nucleotide similarity, and coding potential and carried diverse promoter motifs implicating interactions with different giant viruses. Between host strains, some EMALE insertion loci were conserved indicating ancient integration events, whereas the majority of insertion sites were unique to a given host strain suggesting that EMALEs are active and mobile. Furthermore, we uncovered a unique association between EMALEs and a group of tyrosine recombinase retrotransposons, revealing yet another layer of parasitism in this nested microbial system. Our findings show that virophages are widespread and dynamic in wild Cafeteria populations, supporting their potential role in antiviral defense in protists.

eLife digest

Viruses exist in all ecosystems in vast numbers and infect many organisms. Some of them are harmful but others can protect the organisms they infect. For example, a group of viruses called virophages protect microscopic sea creatures called plankton from deadly infections by so-called giant viruses. In fact, virophages need plankton infected with giant viruses to survive because they use enzymes from the giant viruses to turn on their own genes.

A virophage called mavirus integrates its genes into the DNA of a type of plankton called Cafeteria. It lays dormant in the DNA until a giant virus called CroV infects the plankton. This suggests that the mavirus may be a built-in defense against CroV infections and laboratory studies seem to confirm this. But whether wild Cafeteria also use these defenses is unknown.

Hackl et al. show that virophages are common in the DNA of wild Cafeteria and that the two appear to have a mutually beneficial relationship. In the experiments, the researchers sequenced the genomes of four Cafeteria populations from the Atlantic and Pacific Oceans and looked for virophages in their DNA. Each of the four Cafeteria genomes contained dozens of virophages, which suggests that virophages are important to these plankton. This included several relatives of the mavirus and seven new virophages. Virophage genes were often interrupted by so called jumping genes, which may take advantage of the virophages the way the virophages use giant viruses to meet their own needs.

The experiments show that virophages often co-exist with marine plankton from around the world and these relationships are likely beneficial. In fact, the experiments suggest that the virophages may have played an important role in the evolution of these plankton. Further studies may help scientists learn more about virus ecology and how viruses have shaped the evolution of other creatures.

Introduction

Many eukaryotic genomes harbor endogenous viral elements (EVEs) (Feschotte and Gilbert, 2012). For retroviruses, integration as a provirus is an essential part of their replication cycles, but other viruses also occasionally endogenize, for instance with the help of cellular retroelements (Holmes, 2011). Some green algal genomes even contain giant EVEs of several hundred kilobase pairs (kbp) in length (Moniruzzaman et al., 2020), but unlike prophages in bacteria and archaea, most eukaryotic EVEs are thought to be ‘genomic fossils’ and incapable of virion formation and horizontal transmission. However, some viral genes may be co-opted for various host functions (Frank and Feschotte, 2017; Aswad and Katzourakis, 2012). In recent years, the exploration of protist-infecting giant viruses has uncovered a novel class of associated smaller DNA viruses with diverse and unprecedented genome integration capabilities.

Viruses of the family Lavidaviridae, commonly known as virophages, depend for their replication on giant DNA viruses of the family Mimiviridae and can parasitize them during coinfection of a suitable protist host (La Scola et al., 2008; Krupovic et al., 2016a; Duponchel and Fischer, 2019). A striking example is the virophage mavirus, which strongly inhibits virion synthesis of the lytic giant virus CroV during coinfection of the marine heterotrophic nanoflagellate Cafeteria sp. (Stramenopiles; Bicosoecida) (Fischer and Suttle, 2011; Fischer and Hackl, 2016). Virophages possess 15–30 kbp long double-stranded (ds) DNA genomes of circular or linear topology that tend to have low GC-contents (27–51%) (Fischer, 2020). A typical virophage genome encodes 20–30 proteins, including a major capsid protein (MCP), a minor capsid or penton protein (PEN), a DNA packaging ATPase, and a maturation cysteine protease (PRO) (Krupovic et al., 2016a). In addition to this conserved morphogenesis module, virophages encode DNA replication and integration proteins that were likely acquired independently in different virophage lineages (Yutin et al., 2013). Viruses in the genus Mavirus contain a rve-family integrase (rve-INT) that is also found in retrotransposons and retroviruses, with close homologs among the eukaryotic Maverick/Polinton elements (MPEs) (Fischer and Suttle, 2011). MPEs were initially described as DNA transposons (Pritham et al., 2007; Kapitonov and Jurka, 2006), but many of them carry the morphogenesis gene module and thus qualify as endogenous viruses (Krupovic et al., 2014).

Phylogenetic analysis suggests that mavirus-type virophages share a common ancestry with MPEs and the related Polinton-like viruses (PLVs) (Fischer and Suttle, 2011; Yutin et al., 2013). We therefore tested the integration capacity of mavirus using the cultured protist Cafeteria burkhardae (formerly Cafeteria roenbergensis; Fenchel and Patterson, 1988; Schoenle et al., 2020) and found that mavirus integrates efficiently into the nuclear host genome (Fischer and Hackl, 2016). The resulting mavirus provirophages are transcriptionally silent unless the host cell is infected with CroV, which leads to reactivation and virion formation of mavirus. Newly produced virophage particles then inhibit CroV replication and increase host population survival during subsequent rounds of coinfection (Fischer and Hackl, 2016). The mutualistic Cafeteria-mavirus symbiosis may thus act as an adaptive defense system against lytic giant viruses (Fischer and Hackl, 2016; Koonin and Krupovic, 2016). The integrated state of mavirus is pivotal to the proposed defense scheme as it represents the host’s indirect antigenic memory of CroV (Koonin and Krupovic, 2017). We therefore investigated endogenous virophages to assess the prevalence and potential significance of virophage-mediated defense systems in natural protist populations.

Here, we report that mavirus-like EVEs are common, diverse, and most likely active mobile genetic elements (MGEs) of C. burkhardae. Our results suggest an influential role of these viruses on the ecology and evolution of their bicosoecid hosts.

Results

Endogenous virophages are abundant in Cafeteria genomes

In preparation of screening for endogenous virophages, we generated high-quality de novo genome assemblies of four cultured C. burkhardae strains (Hackl et al., 2020). These strains, designated BVI, Cflag, E4-10P (E4-10), and RCC970-E3 (RCC970), were isolated from the Caribbean Sea in 2012, the Northwest Atlantic in 1986, the Northeast Pacific in 1989, and the Southeast Pacific in 2004, respectively. We sequenced their genomes using both short-read (Illumina MiSeq) and long-read (Pacific Biosciences RSII) technologies in order to produce assemblies that would resolve 20–30 kb long repetitive elements within the host genomic context. Each C. burkhardae genome assembly comprised of 34–36 megabase pairs with an average GC-content of 70% (Hackl et al., 2020).

To identify endogenous virophages, we combined sequence similarity searches against known virophage genomes with genomic screening for GC-content anomalies. The two approaches yielded redundant results and virophage elements were clearly discernible from eukaryotic genome regions based on their low (30–50%) GC-content (Figure 1A). Each element had at least one open reading frame (ORF) with a top blastp hit to a mavirus protein, with no elements bearing close resemblance to Sputnik or other virophages outside the genus Mavirus. In the four Cafeteria genomes combined, we found 138 endogenous mavirus-like elements (EMALEs, Figure 1B and C; Figure 1—figure supplement 1, Supplementary file 1). Thirty-three of these elements were flanked by terminal inverted repeats (TIRs) and host DNA and can thus be considered full-length viral genomes.

Figure 1 with 1 supplement see all
Endogenous virophages in Cafeteria burkhardae.

(A) GC-content graph signature of a virophage element embedded in a high-GC host genome. Shown is a region of contig BVI_c002 featuring an integrated virophage (pink box) flanked by host sequences. (B) Location of partial or complete virophage genomes and Ngaro retrotransposons in the genome assemblies of C. burkhardae strain BVI (see Figure 1—figure supplement 1 for all four strains). Horizontal lines represent contigs of decreasing length ordered from left to right and from top to bottom, with numbers shown for the first contig of each line; colored boxes indicate endogenous mavirus-like elements (EMALEs). Fully assembled elements are framed in black. Ngaro retrotransposon positions are marked by black symbols; open symbols indicate Ngaros integrated inside a virophage element. (C) Graphic summary of the number and types of all EMALEs identified in each of the four C. burkhardae strains. (D) Nucleotide contributions of EMALEs and Ngaros to Cafeteria genomes. Fractions for each strain are computed based on nucleotides in the assembly (left) and nucleotides in the reads (right) mapping to the different parts of the assembly.

The remainder were partial virophage genomes that were located at contig ends or on short contigs. These cases arise from incomplete assembly rather than from biological truncations, since the assembly algorithm probably terminated due to the presence of multiple identical or highly similar EMALEs within the same host genome – a well-known issue for repetitive sequences (Kolmogorov et al., 2019). With 55 elements, C. burkhardae strain BVI contained nearly twice as many EMALEs as any of the other strains, where we found 27–29 elements per genome (Figure 1C, Supplementary file 1). Compared to the total assembly length, EMALEs composed an estimated 0.7–1.8% of each host assembly (Figure 1D). Contributions calculated from assemblies deviated only slightly (0.01–0.3%) from read-based calculations. Therefore, the assemblies seem to provide a good representation of the actual contribution of EMALEs to the overall host genomes.

EMALEs are genetically diverse

From here on, we focus our analysis on the 33 complete EMALE genomes, which were 5.5–21.5 kb long with a median length of 19.8 kb, and TIRs that varied in length from 0.2 to 2.3 kb with a median of 0.9 kb (Supplementary file 1, Figure 3—figure supplement 1). Their GC-contents ranged from 29.7% to 52.7%, excluding retrotransposon insertions where present. To classify EMALEs we used an all-versus-all DNA dot plot approach (Figure 2). It revealed two main blocks: The first block contained EMALEs with GC-contents of 29.7–38.5% (median 35.3%), whereas EMALEs in the second block had GC-contents ranging from 47.2% to 52.7% (median 49.3%). The C. burkhardae EMALEs can thus be roughly separated into low-GC and mid-GC groups.

Figure 2 with 2 supplements see all
Classification of endogenous virophages based on DNA dot plot analysis.

The self-versus-self DNA dot plot of concatenated sequences of 33 complete EMALE genomes and mavirus reveals two main block patterns, corresponding to EMALEs with low (29–38%) GC-content and medium (47–53%) GC-content. Smaller block patterns define EMALE types 1–8. EMALE identifiers indicate the host strain and contig number where the respective element is found. Multiple EMALEs on a single contig are distinguished by terminal letters. Elements printed in bold represent the type species shown in Figure 3. Inset: GC-content distribution of complete and partial EMALEs labeled ‘complete: TRUE/FALSE’. Some partial EMALEs were too short for type assignment and are thus inconclusive. Retrotransposon insertions, where present, were removed prior to analysis.

Based on the similarity patterns within each block, we further distinguish eight EMALE types, with low-GC EMALEs comprising types 1–4 and mid-GC EMALEs comprising types 5–8 (Figure 2). Representative genome diagrams for each EMALE type are shown in Figure 3, for a schematic of all 33 complete EMALEs, see Figure 3—figure supplement 1. According to this classification scheme, the reference mavirus strain Spezl falls within type 4 of the low-GC EMALEs (Figures 2 and 3, Figure 3—figure supplement 1). Partial EMALEs were classified based on their sequence similarity to full-length type species (Figure 2—figure supplement 1).

Figure 3 with 6 supplements see all
Genome organization of eight EMALE types found in Cafeteria burkhardae.

Shown are schematic genome diagrams of the EMALE type species 1–8; for all 33 complete EMALEs, see Figure 3—figure supplement 1. The reference mavirus genome with genes MV01-MV20 is included for comparison. Homologous genes are colored identically; genes sharing functional predictions but lacking sequence similarity to the mavirus homolog are hatched. Open reading frames are numbered individually for each element. Ngaro retrotransposon insertion sites are indicated where present. The dotted line between EMALE01 and EMALE02 separates a homologous region (left) from unrelated DNA sequences (right) and thus indicates the location of a probable recombination event.

The codon and amino acid composition of EMALE genes clearly correlated with the overall GC-content of the EMALE genomes (Figure 2—figure supplement 2). For each encoded amino acid, we observed a strong shift toward synonymous codons reflecting the overall GC trend, and across amino acids, we observed a shift from those encoded by high-GC codons to those encoded by low-GC codons in low-GC EMALEs and vice versa. This uniform trend across all amino acids likely indicates that selection and evolutionary processes driving GC-content variation in these viruses act on the nucleotide level, rather than on the encoded proteins.

With few exceptions, EMALEs are predicted to encode 17–21 proteins each. None of the encoding genes was found to contain introns. The virion morphogenesis module in EMALE types 1 and 3–7 consists of the canonical virophage core genes corresponding to MCP, PEN, ATPase, and PRO proteins. Type 2 EMALEs likely encode a different set of capsid genes as discussed below, and the truncated EMALE type 8 lacks recognizable morphogenesis genes. Another highly conserved gene in EMALE types 1 and 3–7 is MV14, which is always found immediately upstream of the ATPase (Figure 3, Figure 3—figure supplement 1) and codes for a protein of unknown function that is part of the mavirus virion (Born et al., 2018). MV14 is present in various metagenomic virophage sequences (Paez-Espino et al., 2019) and, based on synteny and protein localization, likely encodes an important virion component in members of the genus Mavirus.

The replication/integration module consists of the rve-INT gene and at least one additional ORF coding for a primase/helicase and a DNA polymerase. Low-GC EMALEs encode a mavirus-related primase/helicase and protein-primed family B DNA polymerase (pPolB) (Figure 3, Figure 3—figure supplement 1). Mid-GC EMALEs, on the other hand, lack the pPolB gene and feature a longer primase/helicase ORF that may include a DNA polymerase domain similar to the helicase-polymerase fusion genes described in PLVs (Krupovic et al., 2016b).

Other mavirus genes frequently found in EMALEs include MV19 (encoding a putative protease domain), and two genes of unknown function, MV08 and MV12. Interestingly, all mid-GC EMALEs encode a predicted tyrosine recombinase (YR) in addition to the rve-INT and thus possess two predicted enzymes for genome integration. YRs have been found in other virophages and likely catalyze integration into giant virus genomes (Desnues et al., 2012; Yutin et al., 2015). Notable genes unique to one EMALE type include a putative DNA methylase and a ribonucleotide reductase small subunit gene found in EMALE07. The Tlr6F protein encoded by EMALE types 1 + 2 is present in diverse MGEs, including other virophages, PLVs, and large DNA viruses of the phylum Nucleocytoviricota (Koonin and Krupovic, 2017; Stough et al., 2019).

In general, genes were syntenic between EMALEs of the same type, whereas gene order was poorly conserved among EMALEs of different types, with the following exceptions: MCP was always preceded by PEN, and ATPase was always preceded by MV14, whereas the MV14-ATPase-PRO-PEN-MCP morphogenesis gene order as seen in mavirus was present only in EMALE types 4–7. EMALE02 represents an interesting case, as it shares 6–7 kb of its 5’ part (we chose the primase/helicase genes to mark the 5’ end of all EMALEs) with EMALE01, while the remaining 11 kb are not closely related to other EMALEs or virophages (Figure 3—figure supplement 2). Genes encoded in the latter region are mostly ORFans, with the exception of an MV12-like gene and divergent MCP and ATPase genes with remote similarity to PLVs (Bellas and Sommaruga, 2021) and adintoviruses (Starrett et al., 2021). EMALE02 may thus be the result of a recombination event that exchanged the canonical virophage morphogenesis module of EMALE01 with capsid genes of a PLV (Figure 3, dashed line). Overall, these observations support the notion that recombination and non-homologous gene replacement are important factors in virophage genome evolution (Yutin et al., 2013).

Core gene conservation and non-homologous gene replacement in EMALEs

To validate our classification scheme for EMALEs and to place them in a phylogenetic context to other virophages, we used maximum likelihood reconstruction on the core proteins MCP, PEN, ATPase, and PRO, as well as on rve-INT (Figure 4). In the resulting phylogenetic trees, EMALE core proteins formed monophyletic clades with mavirus and related sequences from environmental samples, thus significantly expanding the known diversity of the genus Mavirus. The environmental sequences that clustered with EMALE core proteins include a single amplified genome (SAG) from an uncultured chrysophyte (Castillo et al., 2019), the metagenomic Ace Lake Mavirus (ALM) (Zhou et al., 2013), and four additional metagenomes that were identified in a global survey of virophage sequences (Paez-Espino et al., 2019). The chrysophyte SAG is nearly identical to mavirus strain Spezl and indicates that the host range of mavirus extends beyond bicosoecids. The metagenomic sequences either clustered with one of the EMALE types, or branched separately from them, which suggests the existence of additional sub-groups (e.g. M590M2_1006461).

Figure 4 with 1 supplement see all
Phylogenetic reconstruction of conserved EMALE proteins.

Unrooted maximum likelihood trees were constructed from multiple sequence alignments of the four virophage core proteins major capsid protein (MCP), penton protein (PEN), ATPase, and protease (PRO), as well as of the retroviral integrase. Nodes with bootstrap values of 80% or higher are marked with dots. EMALEs are color-coded by type; cultured virophages are printed in bold. ALM, Ace Lake Mavirus; DSLV, Dishui Lake virophage; OLV, Organic Lake virophage; RVP, rumen virophage; TBE/TBH, Trout Bog Lake epi-/hypolimnion; YSLV, Yellowstone Lake virophage. Metagenomic sequences starting with Ga and M590 are derived from Paez-Espino et al., 2019.

Within the Mavirus clade, EMALEs of a given type were monophyletic for each of the four core proteins, which corroborates their dot plot-based classification. It is worth noting that although EMALEs of types 5 and 6 are largely syntenic (Figure 3, Figure 3—figure supplement 1), they were clearly distinguishable in their phylogenetic signatures (Figure 4). A comparison of clade topologies revealed that even within the conserved morphogenesis module, individual proteins differed with regard to their neighboring clades, and low-GC and mid-GC EMALEs did not cluster separately from each other. These observations could suggest that the morphogenesis modules of different EMALE types diversified simultaneously and that adaptation of GC-content may occur rather quickly.

In contrast, phylogenetic analysis of rve-INT proteins revealed separate clades for low-GC and mid-GC EMALEs (Figure 4). Each of these clades was affiliated with different cellular homologs that included MPEs and retroelements. Notably, the rve-INT genes of low-GC EMALEs were located near the 5’ end of the genomes, whereas in mid-GC EMALEs, they were located near the 3’ end (Figure 3, Figure 3—figure supplement 1). These observations suggest that EMALEs encode two different rve-INT versions, one specific for low-GC EMALEs that co-occurs with the pPolB and a shorter primase/helicase ORF, and one specific for mid-GC EMALEs that co-occurs with a longer primase/helicase ORF. The two integrase versions may have been acquired independently, or one version could have replaced the other during EMALE evolution. Such non-homologous gene replacement appears to have taken place among the primase/helicase genes, too, as previously noted for virophages in general (Yutin et al., 2013). EMALEs encode several different versions of primase/helicase genes with a degree of amino acid divergence that precluded their inclusion in a single multiple sequence alignment. The YR proteins encoded by EMALE types 5–8 formed a monophyletic clade and were part of a larger group of recombinases that included virophages from freshwater metagenomes, as well as microalgae and algal nucleocytoviruses (Figure 4—figure supplement 1).

Cafeteria strains differ in their EMALE composition

The four C. burkhardae strains displayed distinct EMALE signatures: strain BVI had the highest number of virophage elements with 13 complete and 42 partial EMALEs, whereas the other three strains had 6–7 complete and 20–22 partial EMALEs each (Figure 1C, Supplementary file 1). EMALE types 1, 3, 4, 5, and 6 were present in every host strain, EMALE07 was found in all strains except Cflag, and EMALE types 2 and 8 were detected in strains BVI and E4-10 only. We found no evidence for sequence-specific genome integration of EMALEs after inspecting the host DNA sequences that flanked EMALE integration sites, which confirms previous reports of mavirus integration (Fischer and Hackl, 2016). EMALEs were flanked by target site duplications (TSDs) that were predominantly 3–5 bp in length, although some were as short as 1 bp or as long as 9 bp (Supplementary file 1). By comparison, mavirus and MPEs generate 5–6 bp long TSDs upon integration (Fischer and Hackl, 2016; Pritham et al., 2007; Kapitonov and Jurka, 2006).

To assess whether homologous EMALEs were found in identical loci in closely related host genomes, we conducted sequence similarity searches with the flanking regions of each of the 33 fully resolved EMALEs. Whenever these searches returned a homologous full or partial EMALE with at least one matching host flank, we considered the EMALE locus to be conserved in these host strains. We found varying degrees of conservation, with examples shown in Figure 3—figure supplement 3. In 11 cases, an EMALE insertion was conserved in at least two host strains (Supplementary file 1): three EMALE loci were shared by all four strains, four were shared by three strains, and another four were shared by two strains. Based on conserved EMALE loci, strains Cflag and RCC970 were most closely related with nine shared EMALE integrations, which is in line with phylogenetic and average nucleotide identity (ANI) analyses of these strains (Hackl et al., 2020). The four C. burkhardae genomes have ANIs of >99% and thus appear to differ mostly based on their content of EMALEs and other MGEs.

The most parsimonious scenario for the origin of EMALEs that are located in identical loci in different host strains is that they derived from a single integration event. For instance, EMALE03 BVI_101 is orthologous to Cflag_017C and RCC970_016A (Figure 3—figure supplement 3C), which suggests that this element initially colonized the common ancestor of C. burkhardae strains BVI, Cflag, and RCC970. Further cases of redundant EMALEs are Cflag_017B & RCC970_016B (EMALE01) and BVI_029 & RCC970_095 (EMALE06). These elements may thus derive from relatively ancient integration events, whereas 18 of the 33 complete EMALEs represent integrations that were unique to a single host strain (Supplementary file 1). Strain BVI contained 10 of these 18 unique integrations, more than twice as many as any other strain.

The genomic landscape around EMALE integration sites ranged from repeat-free flanking regions to complex host repeats (Figure 3—figure supplement 4). Of the 29 different integration sites represented by the 33 fully resolved EMALEs, 18 were located near repetitive host DNA (within 10 kb from the insertion site). These repeats, in addition to EMALE TIRs, multiple copies of the same EMALE type, and the putative heterozygosity of EMALE insertions, occasionally caused assembly problems, as illustrated in Figure 3—figure supplement 3.

Next, we analyzed whether EMALE insertions interrupted coding sequences of the host. Fifteen integration sites were located within a predicted host gene (13 in exons, 2 in introns), four were found in predicted 3’ untranslated regions, and three were located in intergenic regions (Supplementary file 1). These data show that EMALE insertions may disrupt eukaryotic genes with potential negative consequences for the host. The apparent preference for integration in coding regions could be assembly related, driven by increased accessibility of euchromatin, or linked to host factors that could direct the rve-INT via its CHROMO domain (Gao et al., 2008).

EMALEs are predicted to be functional and mobile

Based on genomic features such as coding potential, ORF integrity, and host distribution, most EMALEs appear to be active MGEs. With the exception of EMALE08 and EMALE02, all endogenous Cafeteria virophages encode the canonical morphogenesis gene module consisting of MCP, PEN, ATPase, PRO, as well as MV14. EMALE02 likely encodes more distantly related capsid genes. Therefore, all EMALE types except EMALE08 should be autonomous for virion formation. In addition, all EMALEs contain at least one predicted enzyme for genome integration, an rve-INT in EMALE types 1–7 and a YR in EMALE types 5–8. EMALEs thus encode the enzymatic repertoire for colonizing new host genomes. Finally, the high variability of EMALE integration loci among otherwise closely related host strains strongly argues for ongoing colonization of natural Cafeteria populations by virophages.

The genomic similarity to mavirus implies that EMALEs may also depend on a giant virus for activation and horizontal transmission. Shared regulatory sequences in virophages and their respective giant viruses suggest that the molecular basis of virophage activation lies in the recognition of virophage gene promoters by giant virus encoded transcription factors (Fischer and Suttle, 2011; Claverie and Abergel, 2009; Legendre et al., 2010). We therefore analyzed the 100 nt upstream regions of EMALE ORFs for conserved sequence motifs using MEME (Bailey et al., 2009). For all type 4 EMALEs, which include mavirus, we recovered the previously described mavirus promotor motif ‘TCTA’, flanked by AT-rich regions. This motif corresponds to the conserved late gene promoter in CroV (Fischer and Suttle, 2011; Fischer et al., 2010), thus possibly indicating that all type 4 EMALEs could be reactivated by CroV or close relatives. EMALEs of other types lacked the ‘TCTA’ motif, but contained putative promoter sequences that may be compatible with different giant viruses (Figure 3—figure supplement 5).

MGEs are prone to various decay processes including pseudogenization, recombination, and truncation. Among the 33 fully resolved EMALEs are three truncated elements: Cflag_215 and RCC970_122 (both EMALE04), and BVI_005 (EMALE08) (Figure 3—figure supplement 1). Interestingly, even these shorter elements are flanked by TIRs, which must have regenerated after the truncation event. Whereas most EMALE ORFs appeared to be intact, as judged by comparison with homologous genes on syntenic elements, several EMALEs contained fragmented ORFs (e.g. ATPase and PEN genes in EMALE04 BVI_055B, Figure 3). To test whether premature stop codons may be the result of assembly artifacts, we amplified selected EMALEs by PCR and analyzed the products using Sanger sequencing. When we compared the Sanger assemblies with the Illumina/PacBio assemblies, we noticed that the latter contained several substitutions and indels. For example, the MCP gene of EMALE01 RCC970_016B was split into three ORFs in the Illumina/PacBio assembly, whereas a single ORF was present in the corresponding Sanger assembly (Figure 3—figure supplement 6). None of the Sanger assemblies confirmed fragmentation of conserved virophage genes, emphasizing the importance of independent sequence validation. However, since it was not possible to re-sequence all potential pseudogene locations, we cannot exclude that some EMALE genes may be fragmented. Overall, EMALE ORFs appeared to be intact and are thus likely to encode functional proteins.

Tyrosine recombinase retrotransposons integrate into EMALEs

Strikingly, about one in four virophage elements was interrupted by a GC-rich sequence with a typical length of ~6 kb (Figure 5A, Supplementary file 1). We identified these insertions as retrotransposons of the Ngaro superfamily, within the DIRS order of retrotransposons (Wicker et al., 2007). Ngaro retrotransposons feature split direct repeats with A1–[ORFs]–B1A2B2 structure (Figure 5D), encode a YR instead of the rve-INT that is typical for retrotransposons, and are found in various eukaryotes from protists to vertebrates (Poulter and Goodwin, 2005). In the four C. burkhardae strains, we annotated 80 Ngaro elements, with 10–25 copies per strain (Supplementary file 1). In addition, we found isolated AB repeats scattered throughout the genome, which could have arisen from recombination of 5’ and 3’ repeats and are reminiscent of solo long terminal repeats of endogenous retroviruses (Friedli and Trono, 2015). Solo AB repeats were also occasionally present in EMALEs (BVI_016, BVI_115, Cflag_024, Cflag_214, Cflag_215). Dot plot analysis of concatenated Cafeteria Ngaro sequences revealed four distinct types that showed no similarity at the nucleotide level, but appeared to share the same coding potential (Figure 5B and D). Based on synteny to previously described Ngaros, ORF1 may encode a Gag-like protein; ORF2 encodes predicted reverse transcriptase (Pfam PF00078) and ribonuclease H domains (Pfam PF17919); and ORF3 encodes a predicted YR (Pfam PF00589) with the conserved His-X-X-Arg motif and catalytic Tyr residue. Ngaro YRs are related to putative transposons of bacteria and eukaryotes (Figure 5—figure supplement 1), but bear no sequence similarity to the EMALE-encoded YRs.

Figure 5 with 3 supplements see all
Ngaro retrotransposons in Cafeteria burkhardae.

(A) Genomic profile of an EMALE-integrated Ngaro element showing a GC-content graph (top), open reading frame (ORF) organization of EMALE and Ngaro (middle), and a schematic overview of the three genomic entities (bottom; host: gray, EMALE: blue, Ngaro: red). (B) Self-versus-self DNA dot plot of 80 concatenated Ngaro sequences. Block patterns define Ngaro types 1–4. Ngaros are numbered according to Supplementary file 1, with red numbers indicating retrotransposons inserted in EMALEs. (C) Distribution of Ngaro integration loci in EMALE and host DNA. Ngaro types 1 and 2a show a clear preference for EMALE loci, in contrast to Ngaro types 2b, 3, and 4 that are mostly found in host loci. (D) Coding potential of C. burkhardae Ngaro retrotransposons, shown for one example per type with their host strain and contig numbers listed. Triangles indicate direct repeats. GAG, group specific antigen; RT, reverse transcriptase; RH, ribonuclease H; YR, tyrosine recombinase.

Interestingly, the four Ngaro types in C. burkhardae differed with regard to their integration site preference. Whereas 91% of type 1 Ngaros were inserted in an EMALE, this was the case for only 14% and 8% of type 3 and type 4 Ngaros, respectively. At first glance, type 2 Ngaros were distributed evenly among viral and eukaryotic DNA (47% inside EMALEs, 53% outside EMALEs). However, several type 2 Ngaros were truncated at their 5’ end, featuring ~2 kb long deletions that covered ORF1. All of these deletion variants, which we designate as type 2b, were located in host DNA, in contrast to 82% of the full-length type 2a Ngaros that were inserted in EMALEs (Figure 5C and D). Similarly, most of the type 1 Ngaros that were inserted into eukaryotic chromatin also lacked ORF1. Hence, it appears that ORF1 determines the integration site specificity of Ngaro retrotransposons. We did not observe any preference of host-genome-integrated Ngaros for regions with a GC-content lower than the host average. Out of 20 Ngaro insertion sites with analyzable EMALE flanking regions (including partial EMALEs), 9 were located in intergenic regions, 7 in EMALE genes (1× primase/helicase, 1× pPolB, 1× MV12, 4× MV19), and 4 in TIRs (Supplementary file 1). Considering that intergenic regions comprise only 5–10% of an EMALE genome, we notice a significant bias (p < 1e–5, Fisher’s exact test) toward Ngaro integration in intergenic EMALE DNA, which may be caused either by purifying selection of deleterious Ngaro insertions or by a higher preference for Ngaro integration into EMALE intergenic regions, e.g. due to their lower GC-content.

A possible consequence of retrotransposon insertion is the loss of biological activity and subsequent pseudogenization. However, we found that Ngaro-containing EMALEs did not contain more fragmented genes than Ngaro-free EMALEs (Figure 5—figure supplement 2). Whereas the biological properties of Ngaro retrotransposons and their influence on host-virus-virophage dynamics remain to be explored, the EMALE-Ngaro interactions appear to be convoluted. For instance, an EMALE03 genome in strain Cflag is interrupted by two adjacent Ngaro1 insertions, while the EMALE itself is located inside an Ngaro4 element (Figure 5—figure supplement 3).

Discussion

Virophages represent a recently discovered family of eukaryotic dsDNA viruses that possess interesting genome integration properties and have potentially far-reaching eco-evolutionary consequences. Our genomic survey of the marine bicosoecid C. burkhardae revealed an unexpected abundance and diversity of endogenous virophages, with dozens of elements in a single host genome. Based on DNA dot plots and phylogenetic analysis, we distinguish eight different types of mavirus-related endogenous virophages. Similar to mavirus, these EMALEs could potentially reactivate and replicate in the presence of a compatible giant virus. Mavirus is proposed to act as an adaptive defense system against CroV in Cafeteria populations (Fischer and Hackl, 2016; Koonin and Krupovic, 2016), and our findings suggest that different types of EMALEs may respond to different giant viruses infecting Cafeteria. The assortment of endogenous virophages in a given host genome may thus reflect the giant virus infection history of that population (Blanc et al., 2015). Some EMALEs are present in orthologous genomic loci in two or more host strains and likely date back to the common ancestor of these strains. However, at least half of the EMALE insertions are specific to a given host strain and may thus have been acquired relatively recently. Combined with the overall integrity of EMALEs and the conservation of integrase and capsid genes, these findings suggest that endogenous virophages in C. burkhardae are active MGEs.

Similar to other MGEs, EMALEs are likely prone to various decay processes including pseudogenization, truncation, and recombination; however, the extent to which these mechanisms affect EMALE stability is currently unclear. Although nonsense mutations appear to be present in some EMALE genes, re-sequencing of selected genomic regions did not confirm pseudogene loci in those cases (Figure 3—figure supplement 6). We found three truncated EMALEs (5–15 kb long), but these are not to be confused with the 105 partial EMALEs that resulted from computational limitations during sequence assembly of these multi-copy elements. Recombination most likely created the EMALE02 variant, with no obvious signs of degradation. Based on our EMALE identification approach (see Materials and methods), we are quite confident that we have not missed more severely degraded EMALEs or more divergent virophages, except perhaps for hypothetical elements lacking any recognizable sequence similarity to annotated virophages and possessing a GC-content similar to that of the host genome. Overall, these observations raise the question of how balance between new virophage integrations and loss of existing EMALEs is achieved to prevent overload of the host genome with mobile DNA. Possible explanations include loss by homologous recombination, for example, between TIRs of the same or similar EMALEs, excision of EMALEs via host- or EMALE-encoded integrases or endonucleases, or mechanisms to limit the uptake of new elements. Given that multiple mavirus genome integrations per host cell can be observed within a period of days (Fischer and Hackl, 2016) and that virophage integration rates under natural infection conditions are likely to fall within similar timescales, we hypothesize that the majority of EMALE loss may occur before pseudogenization and degradation can take place. EMALEs would thus have relatively short residence times in the unicellular Cafeteria compared to germline EVEs of multicellular eukaryotes (Barreat and Katzourakis, 2021).

Horizontal transmission of virophages in natural environments is likely limited by their requirement for two concurrent biological entities, namely a susceptible host cell infected with a permissive giant virus. Similar to temperate bacteriophages, persistence in the proviral state may thus be an essential survival strategy for virophages, as underscored by the abundance of EMALEs in C. burkhardae. Endogenous virophages may also be common in eukaryotes outside the order Bicosoecales, although a search in 2015 for provirophages in 1153 eukaryotic genomes found only one clear case in the chlorarachniophyte alga Bigelowiella natans (Blanc et al., 2015). In contrast to rapidly increasing virophage reports from metagenomic sequence mining in recent years (Paez-Espino et al., 2019), we propose that the discovery of host-integrated virophages is still hampered not only by sampling bias, but also by technical limitations. For instance, AT-rich mavirus DNA was severely underrepresented when we sequenced the genome of C. burkhardae strain E4-10M1 with the standard Illumina MiSeq protocol (Fischer and Hackl, 2016). Additional problems arise during binning and assembly procedures. Although our sequencing and assembly strategy was specifically tailored to endogenous virophages and resolved 33 EMALEs in their host genomic context, dozens of EMALEs were only partially assembled, some may contain assembly errors (Figure 3—figure supplement 3), and others may have been missed altogether. Advances in long-read sequencing technologies and assembly algorithms will likely alleviate such problems.

Surprisingly, we found that EMALEs were frequently interrupted by Ngaro retrotransposons, which revealed an additional level of nested parasitism in this microbial system. Cafeteria genomes contain four distinct Ngaro types with different affinities for EMALEs. Deletion of ORF1 in type 1 and type 2 Ngaros coincides with a decreased occurrence of these retrotransposons in EMALEs (Figure 5). Syntenic ORFs in other Ngaros are predicted to encode a Gag-like structural protein (Poulter and Butler, 2015), and Gag proteins of several retroviruses have been linked to integration site specificity (Lewinski et al., 2006; Tobaly-Tapiero et al., 2008). The putative Gag proteins of Cafeteria Ngaros may thus influence whether retrotransposon insertion occurs in an EMALE or in eukaryotic chromatin. So far, retrotransposons have not been described for giant DNA viruses or virophages; however, pandoravirus genomes contain DNA transposons (Sun et al., 2015), and a class of 7 kb long DNA MGEs called transpovirons interacts with the particles and genomes of Acanthamoeba-infecting mimiviruses and their virophages, apparently without affecting viral replication (Desnues et al., 2012; Jeudy et al., 2020). It remains to be studied whether Ngaro retrotransposons use reactivated virophages or giant viruses as vehicles for horizontal transmission, and what effect retrotransposon insertion in EMALEs has on the fitness of the virophage, the host cell, and their associated giant viruses.

In conclusion, we show that endogenous virophage genomes are abundant and diverse in the marine heterotrophic protist C. burkhardae. These mavirus-like EVEs appear to be active and dynamic MGEs with significant potential to shape the genome evolution of their hosts. We present evidence for recombination and gene exchange within EMALEs, and a previously unknown affiliation between virophages and YR retrotransposons. Our findings imply an important role for EMALEs in the ecology and evolution of bicosoecids and are in line with the hypothesis that endogenous virophages provide adaptive defense against giant viruses.

Materials and methods

C. burkhardae cultures and genome sequencing

Request a detailed protocol

C. burkhardae strains BVI, Cflag, E4-10P, and RCC970-E3 were cultured, and their genomes sequenced and assembled as described previously (Hackl et al., 2020).

Detection and annotation of EMALEs

Request a detailed protocol

To identify endogenous virophages, we initially took a manual approach and searched for clearly visible drops in GC-content along contigs when visualized in Artemis Release 16.0.0 (Rutherford et al., 2000). Based on our experience we suspected that integrated virophages would have a considerably lower GC-content than the host genomes. For those low-GC regions we then computed DNA dot plots with Gepard (Krumsiek et al., 2007) to detect TIRs, predicted ORFs with Artemis, manually curated them into gene calls, and also manually assigned functional annotations based on blastp searches against known virophages and public databases. This resulted in the identification of 33 manually curated, TIR-containing high-confidence (‘complete’) EMALEs.

Partial EMALEs that could not be fully assembled and were thus located on contig termini were identified based on sequence similarity to annotated virophages and EMALEs using the BLAST suite. Their boundaries were determined based on GC-content or sequence similarity to already annotated TIRs. Partial EMALEs were annotated with prodigal-2.6.3 (Hyatt et al., 2010). We pre-trained one gene model for mid-GC and one for high-GC EMALEs, and called genes on those genomes separately. Where present, we masked integrated Ngaro elements with ‘N’s prior to training and annotation because we observed in trial runs that automated gene prediction could not properly resolve the gene structure of the Ngaros and would produce fragmented and spurious gene calls for these regions. We also explored other gene prediction programs including PHANOTATE (McNair et al., 2019) and VGAS (Zhang et al., 2019), but found pre-trained prodigal to produce results most consistent with manually curated expert annotations for EMALE genomes.

Detection of Ngaro retrotransposons

Request a detailed protocol

Retrotransposon insertions were first noticed in GC-content graphs as 6–7 kb long GC-rich sequences that interrupted AT-rich EMALEs. DNA dot plot analysis of these regions showed split direct repeats, and conserved domain and tblastx searches predicted RNaseH, reverse transcriptase, and tyrosine recombinase domains. ORFs were annotated in Artemis based on the predicted conserved domains. These initially curated protein sequences and direct repeats were then used to seed blast searches against all four C. burkhardae genomes to identify additional Ngaro elements, which were inspected individually and annotated manually.

Nucleotide contributions of EMALEs and Ngaros to Cafeteria genomes

Request a detailed protocol

To quantify how much of each host genome is comprised of EMALEs and Ngaros, we applied two complementary strategies: (i) We compared the number of nucleotides annotated as EMALEs and Ngaros in the assemblies to the overall assembly sizes, and (ii) we quantified the number of nucleotides in the PacBio reads that we could assign to either of the three fractions – host, EMALE, and Ngaro. The latter approach is less prone to assembly biases, such as overestimation of contributions for elements only present in one allele, or underestimation of contribution due to collapsed repeated copies, or elements not assembled because of low coverage. We aligned the PacBio reads to the assemblies using minimap2 v2.16 (-x map-pb) (Li, 2018) and computed the coverage of the different genomic regions with samtools v1.9 (Li et al., 2009). We did not consider Illumina data for the quantification due to an observed sequencing bias against low-GC regions that would lead to lower EMALE estimates.

Codon usage analysis

To analyze possible correlations between EMALE GC-content and the codon composition of their genes, we counted codons of all genes of complete EMALEs with a custom Perl script and visualized their distribution relative to their GC-content with a custom R script.

Assignment of EMALE and Ngaro types

Request a detailed protocol

We first performed pairwise whole-genome comparisons of all 33 EMALEs plus the reference mavirus genome. Next, we concatenated the EMALE genomes according to their nearest sequence neighbors. We then plotted the resulting concatemer against itself with Gepard (Krumsiek et al., 2007) using a word length of 10, and analyzed the similarity patterns. The dot plot-based classification was confirmed by phylogenetic analysis of virophage core genes. A similar approach was used for type assignment of Ngaro retrotransposons.

For partial EMALEs, we assigned types in an automated manner based on the highest cumulative blastx bitscores to typespecies EMALE genomes. EMALEs with cumulative bitscores below 100 were classified as ‘inconclusive’. For validation, the results were visualized with a beta version of gggenomes (https://github.com/thackl/gggenomes, copy archived at swh:1:rev:b6b6d8f23fab20a8fd3d9904b37329c914b263a5; Hackl et al., 2021).

Phylogenetic analysis of EMALE and Ngaro proteins

Request a detailed protocol

For EMALE core genes, multiple amino acid sequence alignments were constructed with MAFFT using the E-INS-i iterative refinement method (Nakamura et al., 2018). Alignments were manually inspected and trimmed to eliminate long insertions and regions of low sequence conservation. The best model and parameters for a maximum-likelihood phylogenetic reconstruction were estimated with modeltest-NG v0.1.6 (Darriba et al., 2020) and the tree was computed with IQ-TREE v2.0 (ATPase: LG + I + G4, MCP: LG + R4+ F, Penton: LG + I + G4+ F, Protease: LG + R4+ F, all: -B 1000) (Minh et al., 2020). The trees were visualized with ggtree v1.14.6 (Yu et al., 2017). For comparison, we also performed Bayesian inference analysis with MrBayes v3.1.2 using the following settings: rates = gamma, aamodelpr = mixed, number of generations = 1 million (Ronquist and Huelsenbeck, 2003).

For EMALE tyrosine recombinases, we generated a multiple amino acid sequence alignment with MAFFT v7.310 (--genafpair) (Nakamura et al., 2018) and an HMM profile (Eddy, 2011) from the four sequences present in the type species genomes (EMALE05-EMALE08). We then identified closest relatives with jackhmmer on HmmerWeb v2.40.0 for two iterations (-E 1 --domE 1 --incE 0.001 --incdomE 0.001 --seqdb uniprotrefprot) (Potter et al., 2018), realigned the EMALE sequences with the 30 best hits (mafft --genafpair), and trimmed the alignments with trimAl (-automated1) (Capella-Gutiérrez et al., 2009). A maximum-likelihood tree was computed with FastTree v2.1.10 (Price et al., 2010). The tree was visualized with ggtree v1.14.6 (Yu et al., 2017).

Phylogenies for the Ngaro YRs were generated the same way as for the EMALE YRs, however, only one iteration of jackhmmer was run, and only the top 20 database hits were included in the final tree.

Comparative analysis of integration sites and their genomic context

Request a detailed protocol

For each of the 33 fully resolved EMALEs, we copied up to 10 kb (<10 kb if the EMALE was located within 10 kb of a contig border) of the 5’ and 3’ host sequences immediately flanking the TIR of the EMALE, and conducted blastn analyses on each of the four C. burkhardae genome assemblies with each flanking region separately. In case there was another EMALE located in the 10 kb flanking regions, the second EMALE was omitted from the BLAST search and only host sequence was included. After locating orthologous, and sometimes paralogous, sites in each host strain, we identified TSDs by comparing empty and EMALE-containing alleles using pairwise sequence alignments of homologous sites. To analyze the genomic context of EMALEs for repetitive DNA, flanking regions were analyzed by dot plot analysis.

Prediction of putative EMALE promoter motifs

Request a detailed protocol

We predicted putative promoter motifs in EMALE genomes by running MEME-suite v5.1.1 on all 100 bp upstream regions of all coding sequences. We identified the three highest-scoring motifs for each EMALE type individually. Putative motifs were further validated by analyzing their occurrences across all EMALE whole genomes, and motifs not consistently present in multiple intergenic regions were excluded.

Validation of EMALE assemblies by PCR and Sanger sequencing

Request a detailed protocol

We designed primer pairs for selected EMALEs to generate overlapping, 700–1100 bp long PCR products. Primer sequences for the validated EMALE01 as shown in Figure 3—figure supplement 6 are listed in Supplementary file 1. PCR products were obtained using 2 ng of genomic DNA template from C. burkhardae strain BVI, Cflag, E4-10, or RCC970 in 50 μl total volume containing 5 μl 10× Q5 Reaction Buffer (NEB, Frankfurt am Main, Germany), 0.5 U of Q5 High-Fidelity DNA Polymerase (NEB), 0.2 mM dNTPs, and 0.5 μM of each primer. Cycling conditions on a TGradient thermocycler (Biometra, Jena, Germany) consisted of: 30 s initial denaturation at 98°C; 35 cycles of 10 s denaturation at 98°C, 30 s annealing at 56–60°C (depending on the melting temperature of the respective primers), and 45 s to 1 min extension at 72°C; followed by a final extension time of 2 min at 72°C. To check for correct product length and purity, 5 μl of each reaction were mixed with loading dye and analyzed on a 1% (w/v) agarose gel containing GelRed (VWR, Darmstadt, Germany). The remaining PCR mix was purified using a QIAquick PCR Purification Kit (Qiagen, Hilden, Germany) according to the manufacturer’s instructions. Sanger sequencing of PCR products was performed at Eurofins Scientific using the LightRun Tube service. Reads were trimmed, assembled, and mapped to their respective reference EMALEs in Sequencher software v5.2.2 (Gene Codes Corporation, Ann Arbor, MI). Due to the presence of repetitive regions within and across EMALEs in a single host strain and resulting low-quality reads or failed PCRs, the Sanger assemblies typically did not cover the entire length of an EMALE.

Effect of retrotransposon integration into EMALEs

Request a detailed protocol

To test if the insertion of a YR retrotransposon triggered the degeneration of the targeted EMALE, we compared the lengths of conserved genes in Ngaro-containing EMALEs to those without transposon across all 138 EMALEs. Protein lengths were obtained from the sequence files and plotted with a custom R script.

Data availability

Genome sequences of the four Cafeteria burkhardae strains have been published previously (PMID 31964893) and are deposited in GenBank (VLTL00000000.1, VLTM00000000.1, VLTN00000000.1, VLTO00000000.1). DNA sequences and genome annotations of endogenous virophages and Ngaro retrotransposons, as well as multiple sequence alignments used for phylogenetic reconstruction are available via (https://github.com/thackl/cb-emales copy archived at https://archive.softwareheritage.org/swh:1:rev:4fa8b383d5f8f207bf2c9d08421fa2239dbbba94). Additional data items are deposited under https://doi.org/10.5281/zenodo.4632783.

The following previously published data sets were used
    1. Hackl M
    2. Barenhoff D
    3. Heider F
    (2020) NCBI GenBank
    ID VLTL00000000.1. Cafeteria roenbergensis strain RCC970-E3, whole genome shotgun sequencing project.
    1. Hackl M
    2. Barenhoff D
    3. Heider F
    (2020) NCBI GenBank
    ID VLTM00000000.1. Cafeteria roenbergensis strain Cflag, whole genome shotgun sequencing project.
    1. Hackl M
    2. Barenhoff D
    3. Heider F
    (2020) NCBI GenBank
    ID VLTN00000000.1. Cafeteria roenbergensis strain BVI, whole genome shotgun sequencing project.
    1. Hackl M
    2. Barenhoff D
    3. Heider F
    (2020) NCBI GenBank
    ID VLTO00000000.1. Cafeteria roenbergensis strain E4-10P, whole genome shotgun sequencing project.

References

    1. Fenchel T
    2. Patterson DJ
    (1988)
    Cafeteria roenbergensis nov. Gen., nov. Sp., a heterotrophic microflagellate from marine plankton
    Marine Microbial Food Webs 3:9–19.

Decision letter

  1. Antonis Rokas
    Reviewing Editor; Vanderbilt University, United States
  2. George H Perry
    Senior Editor; Pennsylvania State University, United States
  3. Frank Aylward
    Reviewer; Virginia Tech
  4. Chantal Abergel
    Reviewer; CNRS-AMU, France

Our editorial process produces two outputs: (i) public reviews designed to be posted alongside the preprint for the benefit of readers; (ii) feedback on the manuscript for the authors, including requests for revisions, shown below. We also include an acceptance summary that explains what the editors found interesting or important about the work.

Acceptance summary:

Approaching the search of novel viruses while in an endogenized stage, rather than as free virions, the study by Hackl et al., reveals a large diversity of complete and fragmented virophage genomes – termed EMALEs – scattered throughout the genomes of four strains of the marine protist Cafeteria. Given that the activation of the integrated virophage mavirus during infection by the giant virus, CroV, has been shown to have a protective effect on the Cafeteria population, this study provides a tantalizing window into the traces of virophage-giant virus¬-protist interactions in the marine environment. Given the enormous diversity of virophages and giant viruses that have been found in metagenomes with no known hosts, this study is a step towards deciphering the biology of these viruses.

Decision letter after peer review:

Thank you for submitting your article "Virophages and retrotransposons colonize the genomes of a heterotrophic flagellate" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, and the evaluation has been overseen by a Reviewing Editor and George Perry as the Senior Editor. The following individuals involved in review of your submission have agreed to reveal their identity: Frank Aylward (Reviewer #1); Chantal Abergel (Reviewer #2).

The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission.

Essential revisions:

The reviewers made several suggestions and comments that, if addressed, would significantly strengthen your manuscript and clarify its presentation, so i recommend that you consider them and revise accordingly. In particular, please address the following comment:

1) All reviewers were curious/raised questions about the evolutionary analyses used to identify EMALE and NGARO endogenization and loss and so some clarification from the authors about the specific points raised in their analytical process would be beneficial to the article.

Reviewer #1 (Recommendations for the authors):

Small note on nomenclature. In a few situations the authors use the terms EMALEs and virophage interchangeably, which can be a bit confusing because it is not clear if they are referring to the same thing or if there is some implied difference here (i.e. line 102). This is especially true in the Results section, where it would be clearer to just use EMALE for everything since that term was defined specifically.

For the 7 Ngaro elements in EMALE genes (line 326) it would be useful to have a full breakdown of the predicted function of these genes here. And is this section a discussion of all Ngaro elements or only those found within the 33 high quality EMALEs? That was unclear to me.

It would be useful to note the Pfam domain of the Ngaro elements the first time they are mentioned in the text (or if Pfam is not preferred for some reason some alternative protein family, such as an NCVOG or Interpro ID). This would help to more precisely define these elements for readers. It is also unclear in the Methods how the Ngaro elements were identified- I understand this probably involved manual curation to some extent, but it would be useful to know what features the authors looked for when performing these analyses (best BLASTp hits, Pfam domain hits, etc).

Are homologs of the Ngaro genes found in giant virus genomes? This could be advantageous to giant viruses due to their potential to inactivate virophage. It would be interesting if giant viruses promoted Ngaro gene proliferation. Proteins for a collection of NCLDV and their Pfam annotations are publicly available on the Giant Virus Database (https://faylward.github.io/GVDB/), so it may be easy enough to check what the distribution of these elements across the diversity of NCLDV looks like. If the authors chose to do a quick survey of Ngaro elements in other NCLDV it would increase the breadth and impact of the study, but this is not strictly necessary.

The authors provide an interesting discussion on why a previous study found many fewer endogenous virophage (Blanc et al., PNAS 2015). I reviewed the methods of this study and found them to be quite robust, so it is rather surprising that more endogenous virophage were not found. Hackl et al., rightfully point out that assembly issues may mask the presence of these elements. Another important point is that many fewer reference virophage were available in 2015, so homology based methods were more limited at this time.

The skew towards EMALE integration in intergenic regions is interesting. Given the accumulation of junk DNA in many eukaryotic genomes I would actually predict that many degraded EMALEs would not be removed through purifying selection, and that they would accumulate (though given the high effective population size of these marine hosts it is possible). Is it possible that this bias is due to the difficulty in identifying degraded EMALEs? Perhaps there are more degraded EMALEs in the genome but they are just harder to detect? (see comment below).

Perhaps I missed it but I could not find information on how the EMALEs were initially identified and delineated in the Methods. The results state that "To identify endogenous virophages, we combined sequence similarity searches against known virophage genomes with genomic screening for GC-content anomalies. The two approaches yielded redundant results and virophage elements were clearly discernible from eukaryotic genome regions based on their low (30-50%) GC-content (Figure 1A)". To me this seems quite general and it is not clear what tools were used- for example, I am assuming that sequence similarity searches were done at the amino acid level, and if so was this done with the genes already predicted from the previous Sci Data paper? This is an important detail since gene prediction algorithms often miss NCLDV and virophage genes. In the Scientific Data publication the gene prediction was done with "the BRAKER pipeline which utilizes BLAST Augustus and GeneMark-ES. Augustus and GeneMark-ES gene models were trained with publicly available transcriptomic data of C. roenbergensis E4-10P as extrinsic evidence". If EMALE genes are extensively pseudogenized or not expressed would this pipeline still predict them? More details on EMALE detection and some commentary on the possible limitations of degraded EMALE prediction would be useful.

EVEs necessarily exist on a spectrum from intact to degraded, and the reality is that the most degraded elements may consist of only a single pseudogene and may not be confidently detected using any method. I often think of this, since eukaryotic genomes are full of junk DNA with ambiguous origins- one hypothesis is that much of it derives from endogenous viruses that subsequently degrade beyond all recognition. Some additional discussion of this be useful for interpreting the consequences of retrotransposition, since one would expect that these would actually inactivate quite a few EMALEs and effectively turn them into junk DNA.

Reviewer #2 (Recommendations for the authors):

I was only wondering why the authors did not perform an experiment to assess the reactivation of the type 4 EMALEs upon CroV infection as it would nicely demonstrate they are functional MGEs. This would also address the need for another giant virus for EMALEs presenting different promoters.

To me it is not mandatory but would nicely complement the present findings.

The reading of this manuscript also raised one question. Since some groups only possess rve integrase and some present an additional Tyr recombinase, I was wondering if one could be meant for integration into the host genome while the second could be for integration into the giant virus genome, as it is the case for Sputnik which only have the Tyr recombinase and only integrates into the giant virus genome.

I also have a naïve question concerning the pseudogenized EMALEs without Ngaros retrotransposon. The authors addressed this question on the role of Ngaros in EMALE pseudogenization but focused on complete versus absent and did not look for remnants of possible Ngaro. Could pseudogenization in the ones without Ngaros be the result of ancient integrations that disappeared, just leaving for instance the A, B signatures in these EMALEs?

I noticed that in Figure 1 legend B is difficult to read and legends for C and D panels are missing.

In Figure S7, I was wondering what was corresponding to the gap between c023 and c119 for E4-10which appears to be covered by GC content analysis. Same in panel B for RCC970, between c188 and c258 for which a possible explanation is provided, yet the GC content seems to cover this gap.

Reviewer #3 (Recommendations for the authors):

This is a very thorough analysis and raises many avenues for future work.

I found myself speculating on how the biology of the system is working in terms of the gain-and-loss of EMALEs, with one idea for an analysis that could be done. I am not proposing this as additional work, but I am curious if the authors had already looked into something along this vein. Presumably complete EMALEs are relatively recent acquisitions while fragmented EMALEs have been inactivated by random processes (point mutations followed by larger insertion-deletions) and are therefore more ancient. Have you looked at making phylogenies for the core morphology module genes in the EMALE fragments to see if they represent more basal types?

Below are some points where some additional details would provide clarification.

– Were there other GC-low regions that did not correspond to EMALEs? It is curious that the genome seems so uniformly GC high.

– Do lone ngaro elements seem to be integrated in GC-low regions?

L410 "either of the three fractions" What are the three fractions? Do the authors mean "host", "EMALE" and "NGARO"?

Why were only PacBio reads used to quantify the percentage of EMALE/NGARO and not also the Illumina data? Was it due to GC bias?

L415 "samtools v1.9(47)" Missing space before the citation.

Figure 1

The size of (B) is a bit small so it is a bit hard to decipher with the symbols.

The description of (B) could do with a better description that states the lines represent contigs in decreasing size ordered from left to right and top to bottom.

The legend is missing descriptions for parts (C and D).

Figure 2

A description in the legend of what the terms in the plot "complete" "FALSE" and "TRUE" refer to is needed.

Figure 3.

Perhaps instead of the dashed line separating the region of homology between types 1 and 2, another line at the 5' end would better show it as a block of homology.

What program was used to produce the genomic maps?

https://doi.org/10.7554/eLife.72674.sa1

Author response

Essential revisions:

The reviewers made several suggestions and comments that, if addressed, would significantly strengthen your manuscript and clarify its presentation, so i recommend that you consider them and revise accordingly. In particular, please address the following comment:

1) All reviewers were curious / raised questions about the evolutionary analyses used to identify EMALE and NGARO endogenization and loss and so some clarification from the authors about the specific points raised in their analytical process would be beneficial to the article.

We would like to thank the reviewers and editors for their time and effort to evaluate and improve our manuscript! The comments were very helpful and we have taken great care to address all points accordingly. In particular, we have added two sections to Materials and methods where we describe our procedure of identifying EMALEs and Ngaros in the host genomes (see our responses to reviewer comments 1.1, 1.8). We also extended the discussion towards potential mechanisms of EMALE decay and loss, and modified Figures 1 and Figure 3—figure supplement 3 (previously Figure S7) for more clarity.

Please note we have renamed and re-arranged the supplemental figures to associate each of them with a main figure.

Reviewer #1 (Recommendations for the authors):

Small note on nomenclature. In a few situations the authors use the terms EMALEs and virophage interchangeably, which can be a bit confusing because it is not clear if they are referring to the same thing or if there is some implied difference here (i.e. line 102). This is especially true in the Results section, where it would be clearer to just use EMALE for everything since that term was defined specifically.

Thank you for the remark – we tried to eliminate any ambiguous text passages throughout the manuscript.

For the 7 Ngaro elements in EMALE genes (line 326) it would be useful to have a full breakdown of the predicted function of these genes here. And is this section a discussion of all Ngaro elements or only those found within the 33 high quality EMALEs? That was unclear to me.

This section refers to all EMALEs, including partial ones, which is now mentioned in the text. The EMALE genes interrupted by Ngaros are primase/helicase, pPolB, MV12, and 4x MV19, which in addition to Supplementary File 1 (previously Table S3) is now also included in the text. The apparent accumulation in MV19 (predicted peptidase domain) is probably insignificant, since at least two of the affected EMALEs are redundant (found in homologous host genomic loci).

It would be useful to note the Pfam domain of the Ngaro elements the first time they are mentioned in the text (or if Pfam is not preferred for some reason some alternative protein family, such as an NCVOG or Interpro ID). This would help to more precisely define these elements for readers. It is also unclear in the Methods how the Ngaro elements were identified- I understand this probably involved manual curation to some extent, but it would be useful to know what features the authors looked for when performing these analyses (best BLASTp hits, Pfam domain hits, etc).

Pfam domains are now listed in the text. We added a paragraph to Materials and methods called “Detection of Ngaro retrotransposons” that describes our procedure.

Are homologs of the Ngaro genes found in giant virus genomes? This could be advantageous to giant viruses due to their potential to inactivate virophage. It would be interesting if giant viruses promoted Ngaro gene proliferation. Proteins for a collection of NCLDV and their Pfam annotations are publicly available on the Giant Virus Database (https://faylward.github.io/GVDB/), so it may be easy enough to check what the distribution of these elements across the diversity of NCLDV looks like. If the authors chose to do a quick survey of Ngaro elements in other NCLDV it would increase the breadth and impact of the study, but this is not strictly necessary.

We initially had a similar idea, but dismissed it after blast searches against public databases containing giant viruses did not reveal any matches. We now also did a quick comparison of the five Ngaro-type species (1,2a,2b,3,4) elements against the genomes in GVDB. We found one weak hit in one genome (mmseqs translated search, 600 bp, 1.5e-4), but the surrounding regions did not contain any repeat structures that would indicate the presence of a Ngaro-type transposon. Hence, we do not have any evidence that giant viruses contribute to Ngaro proliferation.

The authors provide an interesting discussion on why a previous study found many fewer endogenous virophage (Blanc et al., PNAS 2015). I reviewed the methods of this study and found them to be quite robust, so it is rather surprising that more endogenous virophage were not found. Hackl et al., rightfully point out that assembly issues may mask the presence of these elements. Another important point is that many fewer reference virophage were available in 2015, so homology based methods were more limited at this time.

Indeed, reports of virophage sequences have increased in recent years, which we now point out in the discussion.

The skew towards EMALE integration in intergenic regions is interesting. Given the accumulation of junk DNA in many eukaryotic genomes I would actually predict that many degraded EMALEs would not be removed through purifying selection, and that they would accumulate (though given the high effective population size of these marine hosts it is possible). Is it possible that this bias is due to the difficulty in identifying degraded EMALEs? Perhaps there are more degraded EMALEs in the genome but they are just harder to detect? (see comment below).

You probably meant “The skew towards Ngaro integration in EMALE intergenic regions”? This is an interesting topic, of which we currently know very little about in Cafeteria and thus can only speculate.

First, the genome size of Cafeteria is rather small (~35 Mbp) compared to many other eukaryotes, hence massive accumulation of ‘junk DNA’ is somewhat unlikely.

Second, our study predicts that virophage integration in Cafeteria genomes is an ongoing process. In order to prevent ‘parasitic DNA overload’, mechanisms must exist to either prevent further integration of virophage genomes, or to remove EMALEs from the host genome, e.g. by recombination. This may affect intact as well as degraded elements. It is possible that recombination events leading to EMALE removal increase with increasing numbers of EMALEs, thus restoring a balance faster than via pseudogenization and subsequent loss.

Third, it is not clear how many degraded EMALEs are actually present. We found a few truncated elements (e.g. EMALE type 08), and scattered signs of pseudogenization. However, the few EMALEs with apparent pseudogenes that we were able to PCR-amplify and re-sequence with the Sanger technique turned out to be fully intact, and the premature stop codons were likely an assembly artifact in these cases (Figure 3—figure supplement 6, previously Figure S10). Whereas the actual number of pseudogenes in EMALEs is probably not zero, it is likely to be lower than it may seem at first glance.

Fourth, we did not find evidence that Ngaro insertion leads to pseudogenization in EMALEs.

An alternative possibility is that Ngaros spread via EMALE propagation and that the observed Ngaro integration skew towards EMALE intergenic regions is the result of a selective process that preserves EMALE activity.

Finally, we do not think that we missed degraded EMALEs in our assemblies, please see also our answers to questions 1.1, 1.8 and 1.9.

Perhaps I missed it but I could not find information on how the EMALEs were initially identified and delineated in the Methods. The results state that "To identify endogenous virophages, we combined sequence similarity searches against known virophage genomes with genomic screening for GC-content anomalies. The two approaches yielded redundant results and virophage elements were clearly discernible from eukaryotic genome regions based on their low (30-50%) GC-content (Figure 1A)". To me this seems quite general and it is not clear what tools were used- for example, I am assuming that sequence similarity searches were done at the amino acid level, and if so was this done with the genes already predicted from the previous Sci Data paper? This is an important detail since gene prediction algorithms often miss NCLDV and virophage genes. In the Scientific Data publication the gene prediction was done with "the BRAKER pipeline which utilizes BLAST Augustus and GeneMark-ES. Augustus and GeneMark-ES gene models were trained with publicly available transcriptomic data of C. roenbergensis E4-10P as extrinsic evidence". If EMALE genes are extensively pseudogenized or not expressed would this pipeline still predict them? More details on EMALE detection and some commentary on the possible limitations of degraded EMALE prediction would be useful.

We apologize; this vital method section detailing these aspects should have been included in the initial manuscript. As suspected by the reviewer, EMALEs are not well annotated by eukaryotic annotation pipelines, but have to be annotated separately for good results. To do so, we initially took a manual approach, identified putative endogenous viruses based on clearly visible drops in GC-content along contigs when visualized in Artemis, supported by the presence of terminal inverted repeats in dotplots of those regions. If terminal repeats were missing, we used the GC-content to manually define EMALE boundaries.

For the 33 EMALEs with TIRs, we predicted open reading frames in Artemis and manually curated those into gene calls. For the partial EMALEs we called genes using prodigal with models pretrained on the manually annotated complete EMALEs.

Across the 4 Cafeteria strains, annotated EMALEs account for 85-93% regions with GC-content lower than 56% (measured as 2000bp sequence windows with a GC-content <56%). We manually inspected all locations with deviant GC-content and found the remaining larger regions to be accounted for either by rRNA gene clusters or simple tandem repeats, possibly related to telomeric or centromeric regions.

We also performed blastx searches with the newly identified EMALEs as well as a collection of known virophages against the Cafeteria assemblies to verify that we did not miss other integrated virophages that might not be detectable by GC-content. We are also reasonably confident that our approach is able to detect pseudogenized virophages and that our set of detected EMALEs is comprehensive. However, we of course, cannot exclude that we may have missed viral elements that do not exhibit any obvious deviation in GC-content and at the same time do not share any sequence homology to other known virophages. We modified the discussion to include these considerations.

EVEs necessarily exist on a spectrum from intact to degraded, and the reality is that the most degraded elements may consist of only a single pseudogene and may not be confidently detected using any method. I often think of this, since eukaryotic genomes are full of junk DNA with ambiguous origins- one hypothesis is that much of it derives from endogenous viruses that subsequently degrade beyond all recognition. Some additional discussion of this be useful for interpreting the consequences of retrotransposition, since one would expect that these would actually inactivate quite a few EMALEs and effectively turn them into junk DNA.

In addition to our previous answers, we think that the situation of Cafeteria EMALEs is quite different from other eukaryotic EVEs described so far. Comparison of the four Cafeteria strains revealed that only few EMALE integration loci are conserved among otherwise near-identical host strains. This suggests a fast turnover time for EMALEs, which is not the case for other EVEs that are often described as ‘genomic fossils’. Thus, it may be possible that many EMALEs are removed long before they can pseudogenize, e.g. by TIR-based recombination, homologous recombination between chromosomes, or excised by EMALE- or host-encoded enzymes (integrases, endonucleases). These processes would still apply even when the EMALE is degraded, hence pseudogenized remnants of EMALEs may be rarer than for other eukaryotic EVEs.

Such thoughts are still highly speculative at this early stage of genomic exploration of Cafeteria populations, but we have added a paragraph to the discussion summarizing these considerations.

Reviewer #2 (Recommendations for the authors):

I was only wondering why the authors did not perform an experiment to assess the reactivation of the type 4 EMALEs upon CroV infection as it would nicely demonstrate they are functional MGEs. This would also address the need for another giant virus for EMALEs presenting different promoters.

To me it is not mandatory but would nicely complement the present findings.

We thank the reviewer for the constructive feedback! We fully agree with this comment and have already invested a considerable amount of work into testing EMALE reactivation in response to CroV infection. However, we strongly prefer to publish these findings separately, as they are quite comprehensive and would clearly exceed the limits of the current manuscript.

The reading of this manuscript also raised one question. Since some groups only possess rve integrase and some present an additional Tyr recombinase, I was wondering if one could be meant for integration into the host genome while the second could be for integration into the giant virus genome, as it is the case for Sputnik which only have the Tyr recombinase and only integrates into the giant virus genome.

Absolutely. This is certainly an interesting notion that we share (see L. 147). Given that virophages need two biological entities (host and giant virus) for successful horizontal propagation, the strategy of physically tracking one of these entities appears to be essential. So far, Sputnik integration has only been shown in mimiviruses, and mavirus integration only in its cellular host. Some EMALEs may have combined both strategies but unfortunately, at this point we lack any data that would shed more light on the biological significance of EMALE-encoded Tyr recombinases.

I also have a naïve question concerning the pseudogenized EMALEs without Ngaros retrotransposon. The authors addressed this question on the role of Ngaros in EMALE pseudogenization but focused on complete versus absent and did not look for remnants of possible Ngaro. Could pseudogenization in the ones without Ngaros be the result of ancient integrations that disappeared, just leaving for instance the A, B signatures in these EMALEs ?

A valid question. We analyzed all EMALE genomes (partial and full) by blast using Ngaro A and B repeats as input and found 5 cases of EMALEs (mostly partial elements) with solo AB repeats. One of them is already marked in Figure 3—figure supplement 1 (previously Figure S2): Cflag_215, which is a truncated EMALE04, and recombination of the Ngaro repeats could have led to the truncation of that element. We now include this information in the text.

However, we found no evidence for increased pseudogenization in these cases.

I noticed that in Figure 1 legend B is difficult to read and legends for C and D panels are missing.

Thank you for noticing! We modified the layout for Figure 1 and added the accidentally deleted legends of C and D. We also rearranged the figure slightly to allow for more space for panel B.

In Figure S7, I was wondering what was corresponding to the gap between c023 and c119 for E4-10which appears to be covered by GC content analysis. Same in panel B for RCC970, between c188 and c258 for which a possible explanation is provided, yet the GC content seems to cover this gap.

You have a keen eye! Contigs c023 and c119 actually overlap, i.e. the ends of those two contigs align to the same region on BVI-028. The appearance of GC/repeat-content covering gaps is a technical artifact from the plotting program. We now manually removed the extra GC-plot regions from Figure 3—figure supplement 3 (previously Figure S7).

Reviewer #3 (Recommendations for the authors):

This is a very thorough analysis and raises many avenues for future work.

I found myself speculating on how the biology of the system is working in terms of the gain-and-loss of EMALEs, with one idea for an analysis that could be done. I am not proposing this as additional work, but I am curious if the authors had already looked into something along this vein. Presumably complete EMALEs are relatively recent acquisitions while fragmented EMALEs have been inactivated by random processes (point mutations followed by larger insertion-deletions) and are therefore more ancient. Have you looked at making phylogenies for the core morphology module genes in the EMALE fragments to see if they represent more basal types?

We would like to thank the reviewer for the suggestions and the positive assessment of our work!

Yes, we included partial and full EMALEs in our phylogenetic analyses and the resulting patterns were clear: EMALE genes clustered always based on their type, but not based on whether they were found in full or partial EMALEs. That being said, partial EMALEs are not necessarily more degraded (i.e. containing more pseudogenes, insertions or deletions) than full-length EMALEs, but were rather not fully assembled by the computer algorithm (see also answer to question 1.1, reviewer 1). The same pattern was observed for the two truncated EMALEs of type 4 (Cflag_215 and BVI_002), as shown in Figure 4. We did not include the partial EMALEs in these trees for the sake of readability.

Below are some points where some additional details would provide clarification.

– Were there other GC-low regions that did not correspond to EMALEs? It is curious that the genome seems so uniformly GC high.

Across the 4 Cafeteria strains, annotated EMALEs account for 85-93% regions with GC-content lower than 56% (measured as 2000bp sequence windows with a GC-content <56%). We manually inspected all locations with deviant GC-content and found the remaining larger regions to be accounted for either by rRNA genes or simple tandem repeats, possibly related to telomeric or centromeric regions. EMALEs, thus, are the main source of low-GC regions in these genomes.

– Do lone ngaro elements seem to be integrated in GC-low regions?

No, the GC-content of regions surrounding Ngaros integrated into the host genome does not differ from the general host GC spectrum. We now state that in the text.

L410 "either of the three fractions" What are the three fractions? Do the authors mean "host", "EMALE" and "NGARO"?

Yes, we now state that explicitly in the text.

Why were only PacBio reads used to quantify the percentage of EMALE/NGARO and not also the Illumina data? Was it due to GC bias?

Exactly, especially the AT-rich EMALEs in the older Illumina runs for E4-10 were only poorly recovered. Later on we used PCR-reduced library preparation protocols which mitigated those biases somewhat, although not entirely. The GC-bias in the Illumina data would misrepresent host-EMALE ratios, low-vs-high-GC EMALE ratios, and due to differences in protocols also ratios between the strains. We added this information to the manuscript.

Figure 1

The size of B) is a bit small so it is a bit hard to decipher with the symbols.

The description of B) could do with a better description that states the lines represent contigs in decreasing size ordered from left to right and top to bottom.

The legend is missing descriptions for parts C) and D).

Thank you for pointing this out! We appended the legend of Figure 1B and added the accidentally deleted legends of C and D. Figure 1B is now a bit larger, since we moved the symbol legend to the lower panel. However, for close examinations it will be inevitable to zoom into the digital version of the figure.

Figure 2

A description in the legend of what the terms in the plot "complete" "FALSE" and "TRUE" refer to is needed.

Done.

Figure 3.

Perhaps instead of the dashed line separating the region of homology between types 1 and 2, another line at the 5' end would better show it as a block of homology.

What program was used to produce the genomic maps?

Thank you for the suggestion! We tried to implement the idea, but were not happy with the result as it made the figure slightly less transparent. Hence, we decided to keep the minimalistic dashed line. We used Adobe Illustrator to create the figure based on imported Artemis screenshots of the annotated EMALEs for scaling.

https://doi.org/10.7554/eLife.72674.sa2

Article and author information

Author details

  1. Thomas Hackl

    Max Planck Institute for Medical Research, Department of Biomolecular Mechanisms, Heidelberg, Germany
    Contribution
    Conceptualization, Data curation, Formal analysis, Investigation, Software, Validation, Visualization, Writing – review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-0022-320X
  2. Sarah Duponchel

    Max Planck Institute for Medical Research, Department of Biomolecular Mechanisms, Heidelberg, Germany
    Contribution
    Formal analysis, Visualization, Writing – review and editing
    Competing interests
    No competing interests declared
  3. Karina Barenhoff

    Max Planck Institute for Medical Research, Department of Biomolecular Mechanisms, Heidelberg, Germany
    Contribution
    Investigation
    Competing interests
    No competing interests declared
  4. Alexa Weinmann

    Max Planck Institute for Medical Research, Department of Biomolecular Mechanisms, Heidelberg, Germany
    Contribution
    Investigation, Validation
    Competing interests
    No competing interests declared
  5. Matthias G Fischer

    Max Planck Institute for Medical Research, Department of Biomolecular Mechanisms, Heidelberg, Germany
    Contribution
    Conceptualization, Formal analysis, Funding acquisition, Investigation, Supervision, Validation, Visualization, Writing – original draft, Writing – review and editing
    For correspondence
    mfischer@mr.mpg.de
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-4014-3626

Funding

Gordon and Betty Moore Foundation (5734)

  • Sarah Duponchel
  • Matthias G Fischer

Max Planck Institute for Medical Research Heidelberg

  • Thomas Hackl
  • Sarah Duponchel
  • Karina Barenhoff
  • Alexa Weinmann
  • Matthias G Fischer

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

This work was funded by the Max Planck Society and the Marine Microbiology Initiative of the Gordon and Betty Moore Foundation (Grant #5734). We thank I Schlichting and C Roome for support, and A Koslová, J Reinstein, C Bellas, and B Müller for discussions.

Senior Editor

  1. George H Perry, Pennsylvania State University, United States

Reviewing Editor

  1. Antonis Rokas, Vanderbilt University, United States

Reviewers

  1. Frank Aylward, Virginia Tech
  2. Chantal Abergel, CNRS-AMU, France

Publication history

  1. Preprint posted: November 30, 2020 (view preprint)
  2. Received: August 1, 2021
  3. Accepted: September 28, 2021
  4. Version of Record published: October 26, 2021 (version 1)

Copyright

© 2021, Hackl et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 1,155
    Page views
  • 76
    Downloads
  • 0
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Download citations (links to download the citations from this article in formats compatible with various reference manager tools)

Open citations (links to open the citations from this article in various online reference manager services)

Further reading

    1. Developmental Biology
    2. Genetics and Genomics
    Juliet R Girard et al.
    Research Article Updated

    Mechanistic studies of Drosophila lymph gland hematopoiesis are limited by the availability of cell-type-specific markers. Using a combination of bulk RNA-Seq of FACS-sorted cells, single-cell RNA-Seq, and genetic dissection, we identify new blood cell subpopulations along a developmental trajectory with multiple paths to mature cell types. This provides functional insights into key developmental processes and signaling pathways. We highlight metabolism as a driver of development, show that graded Pointed expression allows distinct roles in successive developmental steps, and that mature crystal cells specifically express an alternate isoform of Hypoxia-inducible factor (Hif/Sima). Mechanistically, the Musashi-regulated protein Numb facilitates Sima-dependent non-canonical, and inhibits canonical, Notch signaling. Broadly, we find that prior to making a fate choice, a progenitor selects between alternative, biologically relevant, transitory states allowing smooth transitions reflective of combinatorial expressions rather than stepwise binary decisions. Increasingly, this view is gaining support in mammalian hematopoiesis.

    1. Cancer Biology
    2. Genetics and Genomics
    Andrew S McNeal et al.
    Research Article

    Benign melanocytic nevi frequently emerge when an acquired BRAFV600E mutation triggers unchecked proliferation and subsequent arrest in melanocytes. Recent observations have challenged the role of oncogene-induced senescence in melanocytic nevus formation, necessitating investigations into alternative mechanisms for the establishment and maintenance of proliferation arrest in nevi. We compared the transcriptomes of melanocytes from healthy human skin, nevi, and melanomas arising from nevi and identified a set of microRNAs as highly expressed nevus-enriched transcripts. Two of these microRNAs—MIR211-5p and MIR328-3p—induced mitotic failure, genome duplication, and proliferation arrest in human melanocytes through convergent targeting of AURKB. We demonstrate that BRAFV600E induces a similar proliferation arrest in primary human melanocytes that is both reversible and conditional. Specifically, BRAFV600E expression stimulates either arrest or proliferation depending on the differentiation state of the melanocyte. We report genome duplication in human melanocytic nevi, reciprocal expression of AURKB and microRNAs in nevi and melanomas, and rescue of arrested human nevus cells with AURKB expression. Taken together, our data describe an alternative molecular mechanism for melanocytic nevus formation that is congruent with both experimental and clinical observations.