Introduction
============

Bacteriophages, the viruses that infect and replicate in bacteria, exhibit genomes that are highly dynamic (reviewed by [@evaa089-B23]). They undergo rampant recombination and constant rearrangements. They are also known to play a prominent role in horizontal gene transfers between bacteria (reviewed by [@evaa089-B71]). As such, they can be a source of innovation for the bacteria that are infected as well as for the eukaryotic organisms that carry these bacteria. Aphids (Hemiptera: Aphididae) can host a wide variety of facultative endosymbiotic bacteria ([@evaa089-B29]), in addition to their obligate nutritional endosymbiont *Buchnera aphidicola*. These bacteria are not necessary for reproduction nor survival but can endow their host with a variety of beneficial effects, ranging from survival after heat stress to defense against pathogenic fungi and parasitoid wasps (reviewed by [@evaa089-B56]). Although bacteriophages are ubiquitous wherever bacteria exist, their discovery is scarce among the endosymbiotic bacteria of insects. One such case is the APSE (for ***A**cyrthosiphon **p**isum* **s**econdary **e**ndosymbiont) phage: a lambdoid phage with an isometric head and a short tail (resembling species in the family *Podoviridae*) which was originally isolated from an *Acyrthosiphon pisum* aphid infected by an unidentified secondary endosymbiont ([@evaa089-B73]). This phage has been found so far associated with the gammaproteobacterial symbiont species *Hamiltonella defensa*: an endosymbiont lineage that infects several hemipteran species. In the case of aphids, it is clear that toxin genes carried by this phage confer a protective phenotype against parasitoid wasps by disabling wasp development ([@evaa089-B57]; [@evaa089-B10]).

Microscopic analysis and genome sequencing of this phage, first termed APSE-1 (for bacteriophage 1 from APSE), revealed a circularly permuted and terminally redundant double-stranded DNA molecule of 36,524 base pairs (bp) ([@evaa089-B73]). APSE-1 encodes for a protein product showing low sequence similarity with the shiga-like toxin B subunit of several bacteriophages. Further whole-genome sequencing of a second APSE genome from *A. pisum* (termed APSE-2) as well as specific regions of other APSE phages revealed variation across the so-called "types" ([@evaa089-B51]; [@evaa089-B21]). These "types" are characterized by carrying a unique set of genes including (putative) toxins from three protein families: shiga-like toxin, cytolethal distending toxin (CdtB), and YD-repeat toxin. So far, seven types of APSE have been described ([table 1](#evaa089-T1){ref-type="table"}). They are associated with *Hamiltonella* strains from six aphid species (from two subfamilies) and the whitefly *Bemisia tabaci*.

###### 

Previously Characterized APSE Toxins

  Insect Host              *Hamiltonella* Strain   APSE Type                                   Toxin         Reference
  ------------------------ ----------------------- ------------------------------------------- ------------- --------------------------------
  *Acyrthosiphon pisum*    NA                      APSE-1                                      Shigga-like   [@evaa089-B73]
  *Acyrthosiphon pisum*    5AT                     APSE-2                                      CdtB          [@evaa089-B51]; [@evaa089-B21]
  *Acyrthosiphon pisum*    NY26                    APSE-2                                      CdtB          [@evaa089-B10]
  *Acyrthosiphon pisum*    82B                     APSE-2                                      CdtB          [@evaa089-B51]; [@evaa089-B21]
  *Acyrthosiphon pisum*    ZA17                    APSE-2[^a^](#tblfn2){ref-type="table-fn"}   CdtB          [@evaa089-B46]
  *Acyrthosiphon pisum*    WA4                     APSE-2                                      CdtB          [@evaa089-B46]
  *Acyrthosiphon pisum*    A1A, A2F, A2H           APSE-3                                      YD-repeat     [@evaa089-B21]
  *Acyrthosiphon pisum*    AS3                     APSE-3                                      YD-repeat     [@evaa089-B57]
  *Acyrthosiphon pisum*    AS5                     APSE-3                                      YD-repeat     [@evaa089-B57]
  *Acyrthosiphon pisum*    R7                      APSE-3                                      YD-repeat     [@evaa089-B57]
  *Aphis fabae*            H76                     APSE-3                                      YD-repeat     [@evaa089-B22]
  *Aphis craccivora*       5ATac                   APSE-4[^b^](#tblfn3){ref-type="table-fn"}   Shigga-like   [@evaa089-B21]
  *Uroleucon rudbeckiae*   NA                      APSE-5[^b^](#tblfn3){ref-type="table-fn"}   Shigga-like   [@evaa089-B21]
  *Chaitophorus* sp.       N4                      APSE-6                                      CdtB          [@evaa089-B21]
  *Aphis fabae*            H402                    APSE-6                                      CdtB          [@evaa089-B22]
  *Bemisia tabaci*         NA                      APSE-7                                      CdtB          [@evaa089-B21]

Note.---Summary of the previous works characterizing APSE phage toxins. NA, not available.

Later referred to as APSE-8 ([@evaa089-B10]).

Variants of APSE-1, see [@evaa089-B21].

Differences in the level of protection have been observed in different associations of aphids species/genotypes and *Hamiltonella*--APSE combinations. For example, experimental work has shown that *Hamiltonella* endosymbionts can confer varying levels of resistance against different parasitoid wasp species, and even no defense at all ([@evaa089-B4]; [@evaa089-B40]; [@evaa089-B32]; [@evaa089-B41]). It has been shown that *Hamiltonella* carrying APSE-3 confer a strong defensive phenotype, whereas the ones carrying APSE-2 generally confer moderate defense against parasitoids ([@evaa089-B57]; [@evaa089-B46]; [@evaa089-B10]). This protection is indeed largely dependent on the APSE "type" ([@evaa089-B21]; [@evaa089-B46]), and experiments have conclusively demonstrated that APSE-3, carried by *Hamiltonella defensa* strain AS3, can confer its protective phenotype when transplanted to the naturally occurring APSE-free and nonprotective strain A2C ([@evaa089-B10]). The beneficial fitness effect of APSE-bearing *Hamiltonella* has been shown to be conditional to the presence of APSE and the environmental pressure of parasitoid infection. In one study, the authors found that, although apparently no cost to infection with *Hamiltonella* could be detected in aphids not exposed to parasitism, a significant decline in the frequency of *Hamiltonella* was observed ([@evaa089-B55]). In another study, the loss of APSE led to an increase in intra-aphid *Hamiltonella* abundance and was associated to "severe" deleterious effects on aphid fitness ([@evaa089-B74]).

In the current work, we sought to explore and expand our knowledge on the diversity of APSE phages across aphid species. For this purpose, we assembled and annotated full genomes for 16 APSE phages that infect *Hamiltonella* endosymbionts across ten insect species, including several aphids (from five subfamilies) and the whitefly *B. tabaci*. We performed phylogenetic analyses on *Hamiltonella* and its associated APSE phages in order to understand the diversity and evolutionary trajectory of this defensive phage. Through comparative genomics, we investigated genome rearrangement and more specifically look at the evolution of the toxin cassette. We found evidence suggesting that recombination takes place among APSE types in aphids and that most variation in gene content is observed in two main regions of their genomes. Finally, analysis of the toxin cassettes revealed that mobile elements might be involved in at least some of the variation observed in this symbiotically relevant genomic region.

Materials and Methods
=====================

Aphid Collection, DNA Extraction, and Sequencing
------------------------------------------------

Complete genome sequences of APSE phages were retrieved from the NCBI database for APSE-1 (also known as *Hamiltonella virus APSE1*) and *Hamiltonella* strain 5AT. We then used the genome data available for five *Hamiltonella* strains from the aphid *A. pisum* and two from the whitefly *B. tabaci* to extract integrated APSE phages. Additionally, we gathered sequencing data from eight different aphid species (from four genera belonging to four subfamilies) collected between the years 2008 and 2016: some of which were previously known to host *Hamiltonella* ([@evaa089-B48]). Specimens were kept in 70% ethanol at 6 °C. For whole-genome sequencing, we prepared DNA samples enriched with bacteria following a slightly modified version of the protocol by [@evaa089-B14] as described in [@evaa089-B34]. For this filtration protocol, ∼15 aphids from one colony were pooled together. Extracted DNA was used to prepare custom paired-end libraries in Genoscope as in [@evaa089-B44]. These libraries were sequenced using either 151- or 251-bp paired-end reads chemistry on a HiSeq2500 Illumina sequencer. For full details on specimen collection, species identification, and accession numbers for the samples from which APSE phages were extracted, see [supplementary table S1](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online.

*Hamiltonella* and APSE Genome Assembly and Annotation
------------------------------------------------------

First, we scanned the aforementioned *Hamiltonella* assembled genomes using the sequence of the phage attachment site (also known as attP) and the last 62 bp with a 90% identity threshold in UGENE v1.29.0 ([@evaa089-B54]). For noncircularized genomes, we first scanned the contigs using BlastN ([@evaa089-B2]) for these sequences in order to identify the contig or scaffold where the putative APSE phage resided in. Illumina sequences from the eight newly sampled aphid species were first right-tail clipped (requiring a minimum quality threshold of 20 and a minimum length of 75 bp) using FASTX-Toolkit v0.0.14 (<http://hannonlab.cshl.edu/fastx_toolkit/>, last accessed March 4, 2020). Additionally, PRINSEQ v0.20.4 ([@evaa089-B68]) was used to remove reads containing undefined nucleotides as well as those left without a pair after the filtering and clipping process. The resulting reads were assembled using SPAdes v3.11.1 ([@evaa089-B6]) with the \--only-assembler option and k-mer sizes of 33, 55, 77, 99, and 127. From the resulting contigs, those that were shorter than 200 bp were dropped. The remaining contigs were binned using results from a BlastX ([@evaa089-B2]) search (best hit per contig) against a database consisting of the Pea aphid's proteome and a selection of aphid's symbiotic bacteria proteomes ([supplementary table S2](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online). When no genome was available for a specific bacterial lineage, closely related bacteria were used. The assigned contigs were manually screened using the BlastX web server (searching against the nr database) to ensure correct assignment. This binning process confirmed the presence of *Buchnera*, *Hamiltonella*, and its corresponding APSE phage. APSE was confidently assigned to *Hamiltonella* by checking for chimeric APSE--*Hamiltonella* chromosome contig ends in the corresponding insertion sites and/or by checking that no other facultative symbionts were sequenced in the sample. The resulting contigs were then used as reference for read mapping and individual genome assembly using SPAdes, as described above. *Hamiltonella* draft assemblies included APSE-assigned reads. APSE individual assemblies were visually screened for inconsistencies using Tablet ([@evaa089-B49]).

The resulting APSE genomes underwent a draft annotation using Prokka 1.14.4 ([@evaa089-B69]). This was followed by noncoding RNA prediction using infernal v1.1.2 ([@evaa089-B52]) (against the Rfam v14.1 database, [@evaa089-B36]; [@evaa089-B37]) and tRNAscan-SE v2.0.5 ([@evaa089-B13]). We then performed manual curation of the annotations on UGENE ([@evaa089-B54]) through online BlastX searches of the intergenic regions, open reading frame (ORF) finding feature in UGENE, and through DELTA-BLAST ([@evaa089-B7]) searches of the predicted ORFs against NCBI's nr database and the InterProScan v5 web server ([@evaa089-B33]; [@evaa089-B50]). SignalP v5.0 ([@evaa089-B1]) and Phobius v1.01 ([@evaa089-B35]) were used to predict signal peptides. ORFs were considered to be putative functional proteins (and thus not pseudogenes) if seemingly essential domains for the function were found or if the ORFs displayed truncations but retained identifiable domains. Short ORFs (≤300) that were conserved but consistently pseudogenized in all but one APSE genome were annotated as miscellaneous features. The origin of replication was determined using originX ([@evaa089-B75]). Once the first APSE genome was curated, we used these proteins as a "genus"-specific database in Prokka to annotate the subsequent genomes and iteratively added novel proteins to the database. Codon usage was calculated separately for closed high-quality *Hamiltonella*--APSE pairs using the Cousin web server ([@evaa089-B9]).

Finally, phage-like regions were searched in *Arsenophonus* spp. genomes using the web server for PHASTER ([@evaa089-B77]; [@evaa089-B3]). The retrieved phage-like regions were then manually inspected in UGENE, and identity of genes versus those of APSE was assessed using the online web server of BLAST. The EMBOSS suite program primersearch ([@evaa089-B63]) was used to predict possible amplification targets in the *Arsenophonus* genomes using the primers reported in [@evaa089-B26].

Phylogenetic and Recombination Analyses
---------------------------------------

In order to reconstruct the phylogeny of *Hamiltonella* endosymbionts, we used seven gene sequences (*accD*, *dnaA*, *gyrB*, *hrpA*, *murE*, *ptsI*, and *recJ*) following [@evaa089-B21]. Genes sequences were gathered from the [Supplementary data](#sup1){ref-type="supplementary-material"} for [@evaa089-B45] ([https://doi.org/10.5281/zenodo.2566355](10.5281/zenodo.2566355), last accessed March 4, 2020) then used to identify orthologs in newly sequenced *Hamiltonella* strains using BlastX. Due to the low coverage for *Hamiltonella* draft genomes for two samples (3702 and 3692), we were unable to recover sequences for all seven genes used for phylogenetic inference, and thus these specimens were excluded from the *Hamiltonella* phylogenetic analysis. For each gene, sequence alignments were conducted using MUSCLE v.3.8.31 ([@evaa089-B27]) and visually checked using SeaView v4 ([@evaa089-B28]). We then removed divergent and ambiguously aligned blocks using Gblocks v0.91b ([@evaa089-B70]). A partitioned scheme was selected using PartitionFinder v2.1.1 ([@evaa089-B39]) with one data block defined for each codon position (first, second, and third) in each gene. The best model of evolution of each partition was selected among those available in MrBayes, using the Bayesian Information Criterion metric under a *greedy* algorithm. Finally, we concatenated the resulting alignments and ran a Bayesian phylogenetic inference using the GTR + I + G model in MrBayes v3.2.7 ([@evaa089-B64]) running two independent analyses with four chains each for 1,500,000 generations and checked for convergence.

For performing both phylogenetic inferences and analyzing the genetic differences across APSE, we first ran an orthologous protein clustering analysis using OrthoMCL v2.0.9 ([@evaa089-B42]; [@evaa089-B15]) on the full predicted proteome for the phage genomes ([supplementary table S3](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online). We then extracted the single-copy core proteins for phylogenetic reconstruction (30 protein groups). Sequences were aligned using MAFFT v7.450 (maxiterate 1,000 localpair) ([@evaa089-B38]). Divergent and ambiguously aligned blocks were removed using Gblocks. Substitution model was selected using ModelTest-NG v0.1.6 ([@evaa089-B19]) using the Akaike information criterion criterium. Bayesian phylogenetic reconstruction was conducted with MrBayes as described above by running chains for 300,000 generations and checking for convergence. Given prior evidence for intragenic recombination in APSE proteins ([@evaa089-B20]), we constructed a second data set removing all genes where recombination has putatively occurred. Recombination was tested for running PhiPack ([@evaa089-B11]) on the aligned proteins and results of these analyses, for the putative nonrecombining proteins based on $\Phi_{\omega}$, can be found in [supplementary table S4](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online. Bayesian phylogenetic reconstruction was ran as described above.

To infer the relationships of the 14 toxin-cassette- and lytic-region-proteins (lysozyme and holins) in APSE genomes, we searched for similar sequences using the online BlastP web server versus the NCBI's nr database and selected the top 50 nonredundant hits. These were filtered for a minimum of 70% query coverage and a 35% sequence identity. We aligned each of them using MAFFT and manually removed divergent and ambiguously aligned blocks. For each gene, model selection was done with ModelTest-NG. Bayesian phylogenetic inference was run for 300,000 generations in MrBayes as described above.

All resulting trees were visualized and exported with FigTree v1.4.1 (<http://tree.bio.ed.ac.uk/software/figtree/>, last accessed March 4, 2020). All files used for phylogenetic analyses as well as PHASTER phage search results and comparative genomics files used in this study are available in [https://doi.org/10.5281/zenodo.3764739](10.5281/zenodo.3699313) (last accessed March 5, 2020).

Results and Discussion
======================

APSE Phage Genomes
------------------

We successfully de novo assembled genomes for eight complete APSE phages from eight different aphid species and extracted an additional eight from both closed and draft *Hamiltonella* endosymbionts from the pea aphid *A. pisum* and the whitefly *B. tabaci* ([table 2](#evaa089-T2){ref-type="table"}). Additionally, we recovered draft genomes for all eight *Hamiltonella* strains hosting the newly sequenced APSE phages. All recovered phage genomes have a genome size between 33,476 and 41,343 bp. They have a very conserved G + C content ranging between 42.81% and 45.04% and code for an average of 38 coding sequences (CDSs) with 3 to no pseudogenes ([table 2](#evaa089-T2){ref-type="table"}). Most APSE phages preserve an intact tRNA-Lys-(UC)UU, which has been lost in 5D and ZA17, both of which still keep an internally degraded tRNA pseudogene. An analysis of 17 bacterial hosts and 37 phage genomes found that the presence of tRNAs in these phages tended to correspond to codons that were highly used by the phage's protein coding genes but were rare in the host's ([@evaa089-B5]). However, we found virtually no difference in the lysin codon usage between pairs of *Hamiltonella*--APSE genomes ([supplementary table S5](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online). Nonetheless, we found that the *Hamiltonella* strain 5AT only codes for a tRNA-Lys-CUU, following the pseudogenization of the conserved tRNA-Lys-UUU, which recognizes the most commonly used codon for lysine. Conversely, its associated APSE preserves its own tRNA-Lys-UUU. This suggests that the presence of a tRNA-Lys-UUU in the APSE phage could drive the loss of this tRNA in the bacterial host, which could in turn provide a selective advantage to carrying this phage.

###### 

Genome Assembly and Annotation Statistics

  Insect Host                   *Hamiltonella*                              APSE                                       
  ----------------------------- ------------------------------------------- ------- -------- --------- ------- ------- ---
  *Bemisia tabaci*              MEAM1[^a^](#tblfn5){ref-type="table-fn"}    NA      38,988   NA        43.36   39(2)   4
  *Bemisia tabaci*              MED-Q1[^a^](#tblfn5){ref-type="table-fn"}   NA      38,949   NA        43.34   39(2)   4
  *Cinara confinis*             2801[^b^](#tblfn6){ref-type="table-fn"}     178×    38,949   1,638×    43.94   38(2)   3
  *Cinara* sp. 3046             3046[^b^](#tblfn6){ref-type="table-fn"}     5×      37,531   38×       43.96   38(2)   3
  *Cinara pinimaritimae*        2836[^b^](#tblfn6){ref-type="table-fn"}     27×     35,191   597×      43.84   38(0)   3
  *Cinara cuneomaculata*        2628[^b^](#tblfn6){ref-type="table-fn"}     777×    36,094   6,498×    43.88   37(1)   3
  *Drepanosiphum platanoidis*   3702[^b^](#tblfn6){ref-type="table-fn"}     3×      36,931   111×      43.86   38(2)   3
  *Eriosoma grossulariae*       3692[^b^](#tblfn6){ref-type="table-fn"}     3×      38,420   38×       42.90   42(3)   3
  *Acyrthosiphon pisum*         5AT[^a^](#tblfn5){ref-type="table-fn"}      28.7×   39,867   NA        42.91   41(2)   3
  *Acyrthosiphon pisum*         NY26[^a^](#tblfn5){ref-type="table-fn"}     643×    39,887   NA        42.86   41(2)   3
  *Acyrthosiphon pisum*         5D[^a^](#tblfn5){ref-type="table-fn"}       NA      39,146   NA        42.81   41(2)   3
  *Acyrthosiphon pisum*         ZA17[^a^](#tblfn5){ref-type="table-fn"}     567×    39,145   NA        42.81   41(2)   3
  *Acyrthosiphon pisum*         APSE1[^c^](#tblfn7){ref-type="table-fn"}    NA      36,524   NA        43.89   37(2)   3
  *Acyrthosiphon pisum*         MI47[^a^](#tblfn5){ref-type="table-fn"}     NA      36,522   NA        43.88   37(2)   3
  *Periphyllus testudinaceus*   2671[^b^](#tblfn6){ref-type="table-fn"}     8×      41,343   57×       44.37   38(0)   3
  *Acyrthosiphon pisum*         AS3[^a^](#tblfn5){ref-type="table-fn"}      535×    38,992   NA        45.04   36(0)   3
  *Cinara watanabei*            3293[^b^](#tblfn6){ref-type="table-fn"}     194×    33,476   25,300×   44.46   34(2)   3

Note.---Assembly and genome statistics for complete APSE genomes and their *Hamiltonella* hosts. NA, not available.

Extracted from publicly available *Hamiltonella* genomes.

Sequenced in this study.

Phage isolated from an unidentified secondary symbiont of *Acyrthosiphon pisum* ([@evaa089-B73]).

All APSE have two noncoding RNAs (ncRNAs) that code for small bacterial RNAs (sRNA): a c4 antisense family sRNA and an isrK Hfq binding family sRNA. In the case of the former, a c4 sRNA has been identified in P1 and P7 phages from *Escherichia coli* as a regulator of the *ant* gene, an antirepressor, through the binding to complementary regions in the upstream region of this gene ([@evaa089-B18]). Through the analysis of the putative so-called a′ and b′ antisense regions of the c4 sRNA, we found that they consistently matched the 3′-end region of the distant DNA polymerase gene and the upstream noncoding region of the adjacent BRO family, N-terminal domain-containing protein. Three BRO family proteins (BRO-A, BRO-C, and BRO-D) from the *Bombyx mori nucleopolyhedrovirus* (BmNPV) have been suggested to influence host DNA replication and/or transcription, following the finding that they are able to bind DNA from the host with a stronger affinity for single-stranded than for double-stranded DNA ([@evaa089-B76]). Regarding IsrK, its expression pattern in *Salmonella typhimurium* suggested that this sRNA might be involved in the regulation of virulence ([@evaa089-B59]). Lastly, both APSE phages harbored by *Hamiltonella* infecting *B. tabaci* MEAM1 and MED-Q1 biotypes code for an additional sRNA-Xcc1 family sRNA upstream of a protein of unknown function (ORF6N domain-containing protein) also unique to these APSE. In *Xanthomonas campestris* pv. *vesicatoria*, *sRNA-Xcc1* is under the positive control of two important virulence regulators, suggesting that they might play a role in pathogenesis ([@evaa089-B16]). All this evidence points toward an important regulatory function of these sRNAs in APSE and highlight them as targets to study their downstream regulatory effects in both APSE and *Hamiltonella* gene expression.

APSE Gene Content Variation and Toxin Cassettes
-----------------------------------------------

APSE genomes were highly conserved in both gene content, order, and sequence identity ([supplementary fig. S1](#sup1){ref-type="supplementary-material"}, [table S3](#sup1){ref-type="supplementary-material"}, and file S3, [Supplementary Material](#sup1){ref-type="supplementary-material"} online). The most variable genes, in terms of sequence identity, across APSE were those of a BRO family, N-terminal domain-containing protein (orthogroup APSEcp_004), the DNA polymerase (APSEcp_008), the lyzozyme (APSElp_003), the small subunit of the terminase (APSEcp_020), and the tail proteins (needle, spike, and fiber assembly; APSEcp_028 and APSEcp_033-34). Regarding the DNA polymerase gene, we observed differential loss of the intein domain (a protein intron). The lyzozyme gene was by far the one that showed the most sequence divergence. A phylogenetic analysis suggested two well supported lyzozyme phylogroups (so-called "P13" and "F," [supplementary fig. S2](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online), consistent with previous findings ([@evaa089-B21]). From our analyses, "P13" is restricted to APSE-7 and APSE-1 (including the APSE-1 variants APSE-4 and APSE-5).

One marked difference in CDS content between APSE harbored by *Hamiltonella* symbionts of aphids and those harbored by the whitefly *B. tabaci* was in the protein flanked by the BRO family, N-terminal domain-containing protein and the putative P-loop NTPase. Although all APSE from aphid-infecting *Hamiltonella* code for a putative phage regulatory protein (Rha family), those infecting the whitefly code for an AntA/AntB antirepressor family protein. To our knowledge, the function of both of these proteins is not well understood, and thus the significance of this difference remains to be studied.

Altogether, we identified two main regions of genomic variation in APSE phages and propose an APSE classification based on the composition of these two variable sequence stretches based on [@evaa089-B21], and keeping the main type numbering (X), with a slight revision (X.Yz) based on our findings ([fig. 1](#evaa089-F1){ref-type="fig"}): The first variable region (X) is the toxin cassette, with small variations within the types (Y), and the second one is found around the vicinity of the DNA polymerase gene (z). In regards to the toxin cassettes, we found five different "types" or "gene sets" represented in the newly assembled APSE genomes. We found seven phages that could be classified within APSE-7. These were associated with *Hamiltonella* from two whiteflies and five species of aphids (from two subfamilies). We found that although the CdtB family protein was conserved across APSE-7, its companion hypothetical protein was not. Both proteins in the cassette showed identifiable signal peptides, suggesting their export. The CdtB family protein is also present in APSE-2; however, the structure of the cassette is unlike that of APSE-7. It codes for an additional two putative toxin genes, and in the case of APSE-2.2 subtypes, two flanking hypothetical proteins, all with signal peptides. In the case of APSE-2.1a, we found that the hypothetical proteins flanking CdtB were missing. Nonetheless, we found a different hypothetical protein (containing a signal peptide) flanked by insertion sequence pseudogenes. This suggests that, at least in this case, the new member of the toxin cassette was likely mobilized by the action of transposable elements. As in a previous study by [@evaa089-B21], we found that all APSE-2 and -7 types had both holin family genes (lambda and HP1) followed by a different lysozyme gene phylotype, so-called "F" and "P13," respectively. The APSE-1 type was represented by two nearly identical APSE genomes, both associated with *Hamiltonella* strains from *A. pisum*. Both presented the same four hypothetical proteins, including two previously identified as putative subunits of a shiga-like toxin ([@evaa089-B73]; [@evaa089-B21]). APSE-3 types were found to be associated with *Hamiltonella* hosted by two distantly related aphid species, *A. pisum* and *Periphyllus testudinaceus*. They showed a very conserved cassette, encoding for a hypothetical protein, with no identifiable signal peptide, and an RHS-repeat putative toxin (previously referred to as "YD-repeat-containing" by [@evaa089-B21]) with a signal peptide. Finally, the cassette found in 3293 (*Hamiltonella* from *Cinara watanabei*) represents a completely novel type, which we classified as APSE-8. This toxin cassette encodes for a single putative virulence-associated protein with no recognizable signal peptide. Previous studies have classified the APSE-2 variant found in *Hamiltonella* strain ZA17 as APSE-8 ([@evaa089-B10]; [@evaa089-B24]; [@evaa089-B17]; [@evaa089-B60]). However, we found that the only difference in its toxin cassette was the pseudogenization of the first putative toxin gene within this region when compared with APSE-2 5AT and NY26. Therefore, we have considered it to be a subtype of APSE-2, in agreement with [@evaa089-B46].

![Variable regions in APSE genomes. Genome synteny plots of the two variable regions across APSE phages. On the bottom left, a color key for the features shown in the plot. On the right, proposed classification of APSE in types and subtypes based on the variable regions composition.](evaa089f1){#evaa089-F1}

Regarding the second variable region, the main differences laid in the presence/absence of the KilA-N domain-containing protein and its companion hypothetical protein, a putative transcriptional regulator upstream of the DNA polymerase gene. How this changes impact the phage fitness and phenotype remains to be investigated.

Although most genes in the cassette had only very distant or no affiliations to genes in the databases, both the RHS-repeat and CdtB toxins had many close hits in the NCBI's nr database (see Materials and Methods). In both cases, they formed a well-supported monophyletic group ([fig. 2](#evaa089-F2){ref-type="fig"}). As aforementioned, the RHS-repeat toxin is, to the best of our knowledge, restricted to APSE-3 types. On the contrary, the CdtB protein is actually present in several APSE phages, namely APSE-7 and APSE-2. The CdtB proteins encoded in the two different APSE types do not form distinct subclades, rather the proteins of APSE-2 are nested within CdtB from APSE-7 ([fig. 2*B*](#evaa089-F2){ref-type="fig"}). Due to previous evidence for intragenic recombination ([@evaa089-B20]), we ran a recombination test using PhiPack ([@evaa089-B11]) and found no significant evidence for recombination ($\Phi_{\omega}$*P* value = 7.02e-02). However, NSS- and Max $\chi^{2}$-tests did found significant recombination among CdtB genes (*P* value = 4.00e-03 and *P* value = 2.00e-03, respectively), and thus, the relationships among APSE-encoded CdtB proteins should be interpreted with caution.

![APSE toxin gene phylogenies. Bayesian phylogeny of (*A*) RHS-repeat toxins and (*B*) CdtB proteins. Numbers at nodes represent Bayesian posterior probabilities. "\*" = 1. Phylogenies were midpoint rooted.](evaa089f2){#evaa089-F2}

In both cases, these proteins were closely related to CdtB encoded by *Arsenophonus* spp. hosted by two whitefly species: *Aleurodicus floccissimus* and *B. tabaci*. Additionally, another four proteins in the toxin cassettes had high similarity and a close phylogenetic relation to proteins in *Arsenophonus* spp. (orthogroups APSEtp_002, APSEtp_007, APSEtp_013, and APSEtp_014). Previous analyses have identified APSE-like DNA polymerase genes in *Arsenophonus* symbionts of whiteflies, aphids, louse flies, bat flies, and psyllids ([@evaa089-B30]; [@evaa089-B26]). The diagnostic polymerase chain reactions used in these studies targeted the so-called P45 (DNA polymerase) and P3 (virulence-associated protein E family protein). Analysis of three available high-quality *Arsenophonus* genomes (*Arsenophonus nasoniae* strain FIN \[CP038613.1CP038621.1\], *Arsenophonus* endosymbiont of *Nilaparvata lugens* strain Hangzhou \[JRLH01000001.1\], and *Arsenophonus* endosymbiont of *Aleurodicus floccissimus* \[OUND01000000.1\]), revealed that these reported hits actually corresponded to prophage regions in their genomes (PHASTER annotations available in [https://doi.org/10.5281/zenodo.3764739](10.5281/zenodo.3764739), last accessed March 5, 2020). On closer inspection, we found that only one region in the chromosome of *Arsenophonus nasoniae* strain FIN contained both of these proteins. This region, as did others, preserved high sequence identity to one or a few APSE proteins. However, the regions were not in fact APSE phages: that is, they consistently showed a different gene order, insertion of proteins not characteristic of fully sequenced APSE phages, and/or pseudogenized proteins that are otherwise conserved across APSE analyzed in this study (including the DNA polymerase). Although there are currently no full genome sequences for the strains analyzed by [@evaa089-B30] or [@evaa089-B26], this genomic scan, along with the presence of putative toxins in *Arsenophonus* spp. closely related to those present in APSE, strongly suggests that *Arsenophonus* bacteria have indeed been historically infected by close relatives of the *Hamiltonella*--associated APSE phages, and that these phages have left an imprint in the bacterial genomes. The retention of RHS-repeat- and CdtB-APSE-like toxin genes in at least two different strains of *Arsenophonus* suggests that these proteins could in fact endow the *Arsenophonus* symbiont with a protective phenotype. This could in turn provide a positive "APSE-free" fitness effect of carrying the endosymbiont when confronted with the environmental pressure of parasitoid infection. So far *Arsenophonus*, while quite widespread and very diversified in insects ([@evaa089-B53]), has not been credited with conferring protection against parasitoids of its insect hosts.

Phylogenetic History of *Hamiltonella* and APSE
-----------------------------------------------

Phylogenetic analyses supported previous observations ([@evaa089-B20]) that closely related *Hamiltonella* strains were harbored by distantly related aphid species (i.e., aphids belonging to distinct subfamilies) and hosted different APSE types ([fig. 3*A*](#evaa089-F3){ref-type="fig"}). A noteworthy case is that of *Hamiltonella* strains ZA17 and 5D: They belong to two different clades but host nearly identical APSE phages. Conversely, the closely related *Hamiltonella* strain 5D and MI47 host rather different APSE types. APSE-7 is notably present in both *Hamiltonella* clade C and in the distantly related lineages of this symbiont hosted by the whitefly *B. tabaci*. This provides strong evidence and supports previous findings, based on single-gene phylogenies ([@evaa089-B20], [@evaa089-B21]), of horizontal transfer of these viral entities across *Hamiltonella* strains and lineages. The horizontal transfer of phages across lineages of endosymbiotic bacteria has been suggested for *Wolbachia*, where horizontal transfers of the WO phage seem to occur across both related and divergent *Wolbachia* ([@evaa089-B47]; [@evaa089-B8]). This horizontal transfer could be facilitated by the coexistence of two divergent *Hamiltonella* lineages within the same host, similarly to what has been observed for *Wolbachia* ([@evaa089-B67]; [@evaa089-B62]). Indeed, multiple symbiont species and strains of the same species can coexist within the same population of an aphid species ([@evaa089-B66]; [@evaa089-B72]; [@evaa089-B31]; [@evaa089-B65]; [@evaa089-B48]), offering an opportunity for coinfection and transfer or recombination of their phages. Through the reconstruction of a phylogenetic network using the concatenated single-copy shared proteins of APSE, we observe evidence for recombination across APSE types ([fig. 3*B*](#evaa089-F3){ref-type="fig"}). This is a feature that has previously been observed for these phages ([@evaa089-B20]) and is a general feature of phage genomes (reviewed by [@evaa089-B23]). In fact, only 18 out of 30 single-copy shared orthologous groups of proteins showed no significant evidence of intragenic recombination based on a $\Phi_{\omega}$ test ([supplementary table S4](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online). From these 18, nine showed significant evidence of recombination as judged by two additional tests implemented in PhiPack: NSS and Max $\chi^{2}$. A concatenated nucleotide alignment of the nine putative nonrecombinant genes revealed significant intergenic recombination ($\Phi_{\omega}$*P* value = 0.00e + 00). These results suggest genome-wide recombination for APSE phages, which is expected to complicate the inference of phylogenetic relationships. Regardless, Bayesian phylogenies were reconstructed for the three sets of single-copy shared proteins and can be found in [supplementary figure S3](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online. Given the rampant recombination of their genomes, it is therefore hard to conclude on the phylogenetic history of these viruses.

![Phylogenetic relations of *Hamiltonella* and APSE phages. (*A*) Dendrogram of phylogenetic relationships among *Hamiltonella* endosymbionts based on single-copy shared proteins. Numbers at nodes represent Bayesian posterior probabilities. "\*" = 1. (*B*) Phylogenetic network of single-copy shared proteins as calculated by SplitsTree.](evaa089f3){#evaa089-F3}

One common feature of all phylogenetic analyses was that the APSE phages infecting *B. tabaci*-associated *Hamiltonella* strains were consistently highly divergent from those of aphid-associated ones. Additionally, the SplitsTree network reconstruction suggests little recombination (relative to that observed among aphid-associated APSE) between these two groups of APSE infecting different hosts (aphids and *B. tabaci*). This is congruent with the *Hamiltonella* symbiont phylogeny which positions those infecting the whitefly *B. tabaci* as a separate sister clade to that made up of all currently available aphid-hosted ones ([fig. 3*A*](#evaa089-F3){ref-type="fig"}). Previous experimental studies of *Hamiltonella* infections in *Sitobion avenae* aphids have shown that strains of this symbiont establish infections more easily when transferred from the same host species ([@evaa089-B43]). Additionally, these infections were found to be more stable. The same study also found that infection success was higher when the transferred symbiont strain was more closely related to the one that was originally present in the host. All this strongly suggests that these two lineages infecting different host groups can follow separate evolutionary histories, possibly driven by both ecological factors and other barriers limiting the horizontal transfer of the two distinct symbiont lineages.

Conclusion
==========

Facultative symbionts in insects can confer conditional fitness benefits to their hosts ([@evaa089-B61]; [@evaa089-B12]; [@evaa089-B25]), and *Hamiltonella* and its APSE-conferred defense against parasitoid wasps is one of the best studied cases ([@evaa089-B58]). The results presented in this study support previous studies reporting a variable toxin cassette across APSE as well as intragenic recombination. Additionally, the whole-genome sequencing, annotation, and phylogenetic analyses of three new APSE types conducted here have revealed important features of these protective entities: 1) their genomes are highly conserved in terms of gene content with most of the variation localized to the toxin cassette; 2) alterations can occur within a cassette, with pseudogenization events affecting the putative toxin genes (i.e., APSE-2.2a) or IS-mediated insertion of novel genes; 3) APSE harbored by whitefly-associated *Hamiltonella* seem to have undergone a separate evolutionary history, in line with their hosts; and 4) *Arsenophonus* spp. might have historically had association with APSE-like phages, with some of these having left an imprint in their genomes. This study provides the first genome-wide overview of the diversity of five APSE phage types hosted by *Hamiltonella* endosymbionts, improving our understanding of these phages. The features highlighted in this study will be of great value in the design and interpretation of future experiments of APSE-mediated *Hamiltonella* parasitoid defense.

Supplementary Material
======================

###### 

Click here for additional data file.

This work was supported by a PhD grant from the INRAe ECOFA department (to J.R.), the Marie-Curie AgreenSkills+ fellowship program cofunded by the EU's Seventh Framework Programme (FP7-609398 to A.M.-M.), the Marie Skłodowska-Curie fellowship program *H2020-MSCA-IF-2016* (MicroPhan: 746189 to E.J.), and the *France Génomique* National Infrastructure, funded as part of the *Investissemnt d'Avenir* program managed by the Agence Nationale pour la Recherche (ANR-10-INBS-09 to C.C.). This publication has been written with the support of the AgreenSkills+ fellowship program which has received funding from the EU's Seventh Framework Programme (Grant Agreement No. FP7-609398; AgreenSkills+ contract). We are grateful to the genotoul bioinformatics platform Toulouse Midi-Pyrenees (Bioinfo Genotoul) for providing help and/or computing and/or storage resources. The authors are grateful to the CBGP-HPC computational platform. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

 
=

**Data deposition:** The annotated genome sequences from this project have been deposited at the European Nucleotide Archive under the accessions PRJEB15504, PRJEB31183, PRJEB37342, PRJEB37357, and PRJEB37362--PRJEB37365. The low-coverage *Hamiltonella* drafts are available in [supplementary file S4](#sup1){ref-type="supplementary-material"}, [Supplementary Material](#sup1){ref-type="supplementary-material"} online.

[^1]: These authors contributed equally to this work.
