The draft genome of strain cCpun from biting midges confirms insect Cardinium are not a monophyletic group and reveals a novel gene family expansion in a symbiont

Genome sequencing, assembly, and annotation

Culicoides punctatus female midges were collected from Leahurst Campus, University of Liverpool, UK using UV light traps and identified from wing morphology and by cytochrome c oxidase subunit 1 barcoding as in Pilgrim et al. (2017). DNA was extracted from single individuals using the QIAGEN DNAeasy™ Blood & Tissue Kit following the protocol for purification of total DNA from Insect. All samples were tested for Cardinium infection using a PCR assay based on 16S rRNA Cardinium specific primers Car-sp-F 5′-CGGCTTATTAAGTCAGTTGTGAAATCCTAG-3′; Car-sp-R 5′-TCCTTCCTCCCGCTTACACG-3′ (Nakamura et al., 2009). Whole-genome sequencing was carried out by the Centre for Genomic Research, University of Liverpool using the Illumina TruSeq Nano library preparation protocol. Two short-insert (∼550 bp insert size) paired-end libraries were constructed from two pooled DNA samples of three individuals each. The libraries were multiplexed and sequenced using 2/3 of a lane on an Illumina HiSeq 2500 platform, yielding 2 × 125 bp paired reads. Adapter removal and quality trimming of the raw Illumina reads were performed with Cutadapt version1.2.1 (Martin, 2011) and Sickle version 1.2 (Joshi & Fass, 2011).

Identification and filtering of symbiont reads were performed using a similar approach to that used previously (Pilgrim et al., 2017). Briefly, a preliminary assembly of the quality trimmed dataset was performed using SPAdes version 3.7.0 (Nurk et al., 2013) using the following parameters (-k 21,33,55,77, –careful, –cov-cutoff 5). The initial contigs were visualized using taxon-annotated GC-coverage plots (Fig. S1) with Blobtools (Kumar et al., 2013; Laetsch, 2016). Additional tblastx searches (Altschul et al., 1997; Camacho et al., 2009) were conducted against a local genomic database consisting of Cardinium genomes—cBtQ1 and cEper1 endosymbionts of the whitefly B. tabaci and the parasitic wasp E. pergandiella, respectively (Santos-Garcia et al., 2014; Penz et al., 2012), that of Cardinium strain cHgTN10 from H. glycines (Showmaker et al., 2018) and the more distantly related Acanthamoeba endosymbiont A. asiaticus (Schmitz-Esser et al., 2010)—with an e-value cut-off of 1e⁻⁶. Cardinium contigs were extracted and checked for contamination by blastx searches against the non-redundant (nr) protein database. Cardinium-specific reads were subsequently retrieved using Bowtie2 (Langmead & Salzberg, 2012) and samtools (Li et al., 2009) and re-assembled de novo using SPAdes as described above. All contigs larger than 500 bp were checked for potential host or other bacteria contamination using blastx searches against nr database and all contaminant contigs were removed from the final assembly. Subsequently, we evaluated the quality of the assembled contigs using the reference-free assembly validation tool REAPR (Hunt et al., 2013). REAPR uses read pairs mapping information to identify potential assembly errors and assign quality scores on each base of the assembly. The error calls were then used to break the pre-assembled contigs at every potential miss-assembly position using the aggressive option “-a.” Finally, the broken assembly was scaffolded using SSPACE (Boetzer et al., 2011) using the default parameters.

The cCpun draft genome was annotated using Prokka version 1.12 (Seemann, 2014) and the completeness was assessed using BUSCO v3 based on the presence of 148 universal bacterial marker genes (Simão et al., 2015). Clusters of Orthologous Groups (COG) functional categories were assigned using the eggNOG database (Huerta-Cepas et al., 2016) while additional domains were assigned by searches against the Pfam protein database (Finn et al., 2016). The k-mer fraction of the filtered reads were computed with Jellyfish v2.2.3 (Marçais & Kingsford, 2011) and used to determine the repeat fraction of cCpun genome using GenomeScope (Vurture et al., 2017). Finally, comparison of the repeat density (repeats ≥ 200 bp and at least 95% identity) between the Amoebophilaceae genomes was performed using MUMmer-plots (Kurtz et al., 2004).

Ortholog identification, comparative, and phylogenetic analyses

The genome sequences of the three available arthropod-associated Cardinium strains Cardinium hertigii cEper1 (Penz et al., 2012), Cardinium hertigii cBtQ1 (Santos-Garcia et al., 2014) and Cardinium cSfur (Zeng et al., 2018), the two nematode-associated endosymbionts cHgTN10 and cPpe (Showmaker et al., 2018; Brown et al., 2018) and the Acanthamoeba endosymbiont A. asiaticus (Schmitz-Esser et al., 2010) were obtained from GenBank and used for comparative analyses (accession numbers GCF_000304455.1, GCF_000689375.1, GCA_003351905.1, GCA_003176915.1, and GCF_000020565.1, respectively). The genomes of Cyclobacterium marinum DSM 745 (GCF_000222485.1) and Marivirga tractuosa DSM 4126 (GCF_000183425.1), two free living Bacteroides species, were used as outgroup for the phylogenetic analyses (based on Santos-Garcia et al., 2014). All GenBank retrieved genomes were re-annotated using Prokka software as described above in order to mitigate the effect of inconsistencies due to alternative annotation practices. Orthologous groups of proteins were identified between cCpun, cEper1, cBtQ1, cSfur, cHgTN10, cPpe, and A. asiaticus using an all-vs-all BLAST search and Markov Cluster (MCL) clustering approach as implemented in OrthoFinder method (Emms & Kelly, 2015). Core, accessory and strain-specific orthogroups between the five genomes were visualized with an UpSet plot using the UpSetR package (Conway et al., 2017).

Phylogenetic reconstruction was performed on a set of 278 single copy core protein sequences shared between the six Cardinium genomes, the genome of A. asiaticus and two free living Bacteroides species (Cyclobacterium marinum and M. tractuosa) that were used as outgroup. To this end, a super-matrix was generated by concatenating the protein alignments of the 278 core proteins and trimmed with trimAl version 1.4 (Capella-Gutiérrez, Silla-Martínez & Gabaldón, 2009) using the “automated” option. The best substitution model (LG+F+R4) was selected using ModelFinder (Kalyaanamoorthy et al., 2017) and phylogenetic inference was performed using the maximum likelihood (ML) criterion as implemented in IQ-TREE v1.6.6 (Nguyen et al., 2015). The robustness of the inferred tree was finally assessed with the ultrafast bootstrap approximation method as implemented in IQ-TREE using 1,000 replicates (Hoang et al., 2018). Alternative phylogenetic hypotheses were tested by constrained tree searches using the approximately unbiased (AU) test (Shimodaira & Goldman, 2002) as implemented in IQ-TREE v1.6.6. Additionally, the distribution of the phylogenetic signal across the concatenated super-matrix was calculated as described in (Shen, Hittinger & Rokas, 2017). Briefly, for each of the 278 core protein alignments the log-likelihood score for the best ML tree topology under concatenation and an alternative conflicting topology was calculated under the same substitution model (LG+F+R4). The difference in the gene-wise log-likelihood scores (ΔGLS) between the two alternative topologies was used as a measure of the phylogenetic signal and to visualize the proportion of core genes supporting each conflicting phylogeny. Finally, an independent phylogenetic analysis was performed on a subset of 46 core ribosomal proteins in IQ-TREE v1.6.6 as described above in order to further test the robustness of our phylogenetic inference. Phylogenetic trees were drawn and annotated online using the EvolView tool (He et al., 2016).

Analyses of the DUF1703 gene family expansion

Genome analysis revealed an expansion of the DUF1703 gene family. To analyze this expansion further, a protein sequence alignment of the DUF1703 gene family from Cardinium together with selected Open Reading Frames (ORFs) with sequence similarity retrieved as best BLAST hits form NCBI’s nr database was performed using MAFFT v7 and default parameters (Katoh & Standley, 2013). Ambiguously aligned positions were subsequently removed using trimAl version 1.4 and the “automated” option. A ML phylogenetic analyses was performed with IQ-TREE version 1.6.6 and the phylogenetic tree were constructed and annotated as described above. Additionally, a neighbor-net phylogenetic network was inferred from the translated nucleotide alignment of the cCpun DUF1703 paralogs using SplitsTree version 4.12.6 (Huson & Bryant, 2006; Bryant & Moulton, 2004) and default parameters. A pairwise identity and similarity matrix of the cCpun DUF1703 amino acid sequence paralogs were constructed using the Needleman–Wunsch global alignment method and the BLOSUM62 substitution matrix as implemented in EMBOSS package (Rice, Longden & Bleasby, 2000). Putative signal peptides were predicted on the SignalP 4.1 Server (Petersen et al., 2011) using the sensitive D-cutoff settings. Detection of putative recombination events was performed using the RDP4 software package (Martin et al., 2015). Recombination Detection Program (RDP) implements several methods for detecting recombination signals including MaxChi (Smith, 1992), GENECONV (Padidam, Sawyer & Fauquet, 1999), BottScan (Salminen et al., 1995), Chimera (Posada & Crandall, 2001), and RDP (Martin & Rybicki, 2000). Global parameters were as follow: P-value cutoff was set to 0.001 using a Bonferroni correction and significance was evaluated from a permutation test based on 1,000 permutations. Detected signals were considered significant only when they were confirmed by multiple methods. Inference of recombination signals can be particularly misleading when diverse sequences are analyzed. To avoid such misalignment artefacts, the 25 complete DUF1703 paralogs were grouped into three groups on the bases of nucleotide sequences similarity (>65%) and the analyses was repeated for each group separately. Finally, the results were also confirmed with PhiPack implementing the pairwise homoplasy index (PHI) algorithm (Bruen, Philippe & Bryant, 2006). Residue composition and conservation within the core nuclease PD-(D/E)XK site of the DUF1703 homologs were illustrated with sequence logos using the Skylign tool (Wheeler, Clements & Finn, 2014).

Nucleotide sequence accession numbers

The raw reads and the cCpun draft genome assembly have been submitted to the DDBJ/EMBL/GenBank database under the BioProject accession number PRJNA487198 (WGS project QWJI00000000).

General features of cCpun draft genomes

The final assembly of the cCpun draft genome consists of 57 scaffolds larger than 500 bp (N50 = 41.6 kb, largest scaffold = 116 kb) comprising a total size of 1,137,634 bp (52 scaffolds ≥ 1,000 bp) with an average GC content of ∼33% and an average depth of coverage 90× (Table 1; Fig. S2). Overall, the cCpun genome shares many characteristics with those of the previously sequenced Cardinium strains cEper1, cBtQ1, cSfur, cHgTN10, and cPpe including similar genome size of around one Mb and comparable GC content (33.7–38%) (Table 1). No plasmids were inferred based on the presence of scaffolds with atypically higher read coverage compared with the average coverage of the complete assembly, presenting a contrast to the previously sequenced arthropod-associated Cardinium (cEper1 and cBtQ1) (Table 1; Fig. S2). Nevertheless, we were able to detect several regions with sequence similarity to elements of the two plasmids found in cEper1 and cBtQ1. Matching regions were mainly transposases, suggesting that these might be remnants of ancestral plasmid invasion/s. Although absence of plasmids has also been reported previously for A. asiaticus, the sister species of Cardinium clade (Schmitz-Esser et al., 2010), the presence of low-copy-number plasmids in cCpun cannot be ruled out.

A total of 917 protein coding genes were identified with an average length of 993 bp corresponding to a coding density of around 80% (Table 1; Table S1). cCpun harbors a single set of rRNA genes with the 16S separated from 5S and 23S and encode a complete set of 37 tRNA genes. The identification of 117 out of the 148 BUSCO marker genes (BUSCO score = C: 79% (S: 79%, D: 0%), F: 2.7%, M: 18.2%, n: 148) (Fig. S3) was comparable to that observed for the previously sequenced and complete cEper1 cSfur and cHgTN10 genomes, which suggests that cCpun is a near complete genome. Overall, the redundancy in cCpun as assessed through MUMmer-plots is lower than both A. asiaticus and cBtQ1 previously described as highly repetitive (Santos-Garcia et al., 2014) (Fig. S4). K-mer frequency analysis of the Illumina reads estimated the repetitive fraction of cCpun genome to be circa 13%.

Phylogenomic analyses place cCpun as an outgroup of both other insect and nematode Cardinium strains

Recently, a new family named Amoebophilaceae was proposed to include the Cardinium clades as well as the amoeba-associated A. asiaticus (Santos-Garcia et al., 2014). Currently, at least four major phylogenetic clades of Cardinium related bacteria have been described (Nakamura et al., 2009; Edlund et al., 2012) with possible evidence for additional clades (Chang et al., 2010). However, the phylogenetic (evolutionary) relationships between these clades are not clear. Previous phylogenetic studies based on partial 16S rRNA and gyrB sequences failed to provide a consistent phylogenetic placement for the arthropod and the nematode Cardinium clades (Morag et al., 2012; Nakamura et al., 2009).

We established the relationship of this group across a concatenated set of 278 single copy core protein coding genes as well as a subset of 46 ribosomal protein genes shared between the seven Amoebophilaceae genomes. The results of both analyses clearly support the position of the midge Cardinium clade (C) as a sister group to both the other arthropod and nematode Cardinium clades (clades A and B) and confirm that the arthropod-associated Cardinium do not form a monophyletic group (Fig. 1A). Constrained tree tests for two alternative topologies (a) nematode Cardinium as sister group of all other arthropod Cardinium and (b) cCpun and nematode Cardinium as a monophyletic group resulted in significantly worse trees (AU test, p < 0.01). This inference was further supported by analysis of single protein phylogenies (Figs. 1B and 1C). A total of 157 out of the 278 single copy core genes (56%) support the monophyletic grouping of the B group Cardinium strains (cHgTN10, cPpe) with the A group (cEper1, cBtQ1 and cSfur) in exclusion of cCpun (p < 0.001, Fisher’s exact test). In contrast, only 68 genes (24%) support the monophyletic grouping of cCpun with the A group strains while a small subset of genes (n = 52; 19%) supports the monophyletic grouping of cCpun with cHgTN10 and cPpe.

Genome content comparisons estimate both a core Cardinium genome, genes associated with an insect-symbiont lifestyle, and cCpun specific genes and gene families

The OrthoFinder clustering algorithm identified a total of 2,015 ortholog protein clusters across the seven Amoebophilaceae genomes (A. asiaticus, cHgTN10, cPpe, cCpun, cEper1, cSfur, and cBtQ1). The seven genomes share a core of 415 ortholog clusters of which 278 consist of single-copy genes (Fig. 2). The cCpun genome codes for a substantial number of unique proteins (Fig. 2; Table S2). Specifically, among the 812 ortholog clusters predicted for cCpun, 190 clusters—including 204 protein coding genes—were assigned as strain-specific (Fig. 2). Of these genes, 40 were predicted to code for proteins of less than 70 amino acids and likely represent either annotation artefacts or pseudogenised gene fragments.

The majority of cCpun specific proteins, 138 (∼67%), had neither significant matches (blastp, e-value ≤ 10⁻¹⁰) in the NCBI-nr database, nor predicted functional domain. These were assigned as hypothetical proteins. Amongst the remaining 66 predicted cCpun-specific protein clusters, those with ankyrin-repeat domains were particularly well represented in the strain specific set (Table S2). The abundance, diversity and presumably eukaryotic origin ANK repeat containing proteins has long led them to be considered likely to be involved in symbiotic interactions, and this has been demonstrated in a few cases (Siozios et al., 2013; Nguyen, Liu & Thomas, 2014; Voth, 2011; Pan et al., 2008). A total of 46 ANK repeat proteins were present in the cCpun genome, which represents the largest expansion of this gene family in Cardinium, comparable to the expansion of this family in A. asiaticus (54 ANK proteins) (Schmitz-Esser et al., 2010). In total, 20 out of the 46 ankyrin repeat-containing proteins identified in cCpun were not found in the other Cardinium strains, suggesting potential host-specific functions. Among the remaining strain-specific protein clusters, 13 were assigned as putative mobile elements (transposases), three putative transporters, a DNA repair protein RecN, two putative GNAT-family acetyltransferases and a homologue of the hemolysin transporter protein ShlB (Table S2). Finally, a folylpolyglutamate synthase (FolC) homologue involved in the tetrahydrofolylpolyglutamate biosynthesis pathway and a putative riboflavin biosynthesis protein RibBA were also detected. Absence of the complete pathway for the de novo biosynthesis of folate in cCpun suggest that FolC probably participates in the folate salvage pathway (folate to polyglutamate) as suggested also by the presence of a dihydrofolate reductase homologue (De Crécy-Lagard et al., 2007).

Candidate proteins related to the adaptation of Cardinium to arthropod hosts (as opposed to Amoeba and nematode) were identified as being in the four arthropod-associated Cardinium strains (cCpun, cSfur, cEper1, and cBtQ1), and not Amoebophilus and the nematode-associated Cardinium strains (cHgTN10 and cPpe). The four strains from whitefly, wasp, planthopper and midge uniquely share 11 ortholog protein clusters (Fig. 2). Among them we observed the virulence-associated E family protein previously detected in the plasmids harbored by cEper1 and cBtQ1 (Penz et al., 2012; Santos-Garcia et al., 2014) and a nicotinamide mononucleotide transporter.

cCpun possesses both afp-like and type IX secretion systems

Intracellular microbes utilize a variety of specialized protein secretion systems in order to invade and interact with their eukaryote host (Tseng, Tyler & Setubal, 2009; Dale & Moran, 2006). A common characteristic of the Amoebophilaceae genomes is that all encode for a putative afp-like protein secretion system presumably involved in host-microbe interactions (Penz, Horn & Schmitz-Esser, 2010; Penz et al., 2012; Hurst et al., 2007). This system was also observed in the cCpun genome (Fig. 3) (Penz, Horn & Schmitz-Esser, 2010; Penz et al., 2012; Santos-Garcia et al., 2014). The organization of the AFP-like genes clusters is conserved between the four Amoebophilaceae genomes and suggests operon-like structures (Fig. 3).

We additionally identified seven components of the type IX secretion system (T9SS) in cCpun, a system related to gliding motility and pathogenicity in several members of the phylum Bacteroidetes (McBride & Zhu, 2013; McBride & Nakane, 2015). cCpun is the third Cardinium strain reported to retain components of the T9SS system (Santos-Garcia et al., 2014; Zeng et al., 2018). Four of these protein clusters with homology to the core components of the T9SS (GldK, GldL, GldM, GldN) are shared between cCpun, A. asiaticus, cBtQ1, and cSfur while an additional three proteins with homology to the lipoproteins GldD, GldJ, and GldH are uniquely shared between cCpun and A. asiaticus with exception the GldJ which was also found in the cSfur genome in two identical copies (Table S3). More recently, core components of the T9SS secretion system were found on the plasmid of Cardinium cBtQ1 (Santos-Garcia et al., 2014).

Originally described in Flavobacterium johnsoniae, the T9SS is unique among the phylum Bacteroidetes having important role in secretion of proteins involved both in gliding motility and pathogenicity (McBride & Nakane, 2015; Sato et al., 2010). The presence of the Gld homologs in cCpun as well as A. asiaticus supports an ancestral origin of the T9SS machinery which was subsequently lost from cEper1 and the nematode clade (cHgTN10 and cPpe). The functional role of the T9SS components in Cardinium is unknown. The gene set identified as present in the clade is small compared to that known for active T9SSs (which may have more than 18 components). The low number of genes identified may either reflect co-option of other (unidentified) genes into the secretion process, or a function outside of secretion. However, it is tempting to speculate that the T9SS machinery in Amoebophilaceae has progressively been replaced by the AFP-like protein secretion system. This hypothesis is supported by the complete absence of Gld homologs in both cEper1 and the nematode strains, which suggests that the T9SS is dispensable and likely undergoing gradual loss due to genome reduction processes (Toft & Andersson, 2010).

The cCpun genome contains an expansion of the DUF1703 gene family

Expansion and contraction of gene families in microbial genomes constitute a major source of both genetic and functional novelty, contributing to their adaptation to changing environments (Bratlie et al., 2010). Despite a tendency for evolution to eliminate redundancy and streamline genomes, endosymbiotic bacteria and intracellular pathogens often contain multi-gene families. Interestingly, the majority of the expanded gene families in these host-associated microbes encode putative effector proteins enriched in eukaryotic domains including ANK, LRR, and TPR repeats, F-box and U-box domains (Domman et al., 2014; Wu et al., 2004; Siozios et al., 2013; Schmitz-Esser et al., 2010).

Inspection of the cCpun genome revealed the presence of an expansion of hypothetical proteins related to the DUF1703 protein family (Knizewski et al., 2007) not previously observed in other Cardinium genomes, or other heritable microbes. A total of 25 gene paralogs coding for hypothetical proteins of this family were identified (Fig. 4). The DUF1703 family contains a group of modular proteins consisting of an N-terminal AAA-ATPase like domain (Pfam ID: PF09820) and a C-terminal PDDEXK_9 nuclease domain (Pfam ID: PF08011). In addition to the 25 paralogs, six genes were found to contain only the AAA-ATPase like domain whilst two genes contained only the nuclease domain (Fig. 4B). All partial genes were detected near the borders of the cCpun scaffolds and may be artefactually truncated. Our estimate of gene family size is thus conservative.

The members of the DUF1703 gene family display in cCpun are diverse, as attested by an average amino acid identity of just 39% amongst members (Fig. S5). This extensive divergence of paralogs suggests that the expansion of this gene family is not recent. Moreover, the pairwise comparison suggest at least three main expansion waves (Fig. S5). Phylogenetic analysis indicates that all but one of the Cardinium cCpun DUF1703 carrying protein sequences form a single cluster closely related to those found in Simkania, an intracellular bacterium member of Chlamidiales known to be associated with protozoa (Fig. 4A). Notably, two of the simkania’s paralogs are encoded on its pSn plasmid, suggesting possible roots for horizontal dissemination of the DUF1703 genes. The exception is the gene CCPUN_02500, which forms a distinct group with homologs identified in Cardinium strain cSfur and the only intact DUF1703 carrying homolog in cHgTN10, and which is closely related to homologs found in Rickettsia and metagenomically-recovered sequences belonging to uncultured members of the Bacteroidetes and Gammaproteobacteria (Anantharaman et al., 2016).

Although larger, the expansion of the DUF1703 gene family is not unique to the cCpun genome. Amongst the most recently sequenced Cardinium genomes (cSfur and cPpe) we identified smaller expansions of the DUF1703 family (Figs. 4A and 4B). In contrast, the genomes of cEper1, cBtQ1, and cHgTN10 contain only a single gene homolog whilst no homologs were detected in A. asiaticus or free-living relatives (Fig. 4B). Reconstruction of the phylogenetic relationships between the homologs clearly show that members from the same organism group together suggesting that independent expansions took place after divergence from the common ancestor. Surprisingly, the eight paralogs identified in cPpe genome are more closely related to their Rickettsia counterparts than the rest of the Cardinium homologs, indicating possible independent acquisition. Our results suggest that the DUF1703 genes have probably originated in Cardinium after they diverged from A. asiaticus, presumably by horizontal gene transfer (HGT) with later expansion in the lineage leading to cCpun, cSfur, and cPpe.

Phylogenetic network analyses revealed several reticulation events within the DUF1703 gene family in cCpun indicating frequent recombination among gene family members (Fig. 4C). We further investigated the extent of recombination using different methods implemented in RDP4 software (Martin et al., 2015). Due to the limited sequence similarity between the members of the DUF1703 family we restricted our analyses to group of sequences sharing at least 65–70% nucleotide similarities since misalignment artefacts can confound the identification of true recombination signals. We detected evidence of intragenic recombination in all examined groups with multiple methods (Table S4) suggesting that DUF1703 paralogs in cCpun readily recombine. Despite the extensive recombination, no apparent homogenization between the members of this gene family is observed as suggested by the limited sequence similarity and the absence of monophyletic clustering of cCpun paralogs. Overall, our results point to a HGT scenario for the origin of Cardinium DUF1703 gene family with subsequent expansion in the cCpun genome, and variation produced both by mutation and recombination.

To gain a better insight into the role of DUF1703 proteins we sought to investigate the distribution and abundance of proteins containing the AAA-ATPase and PDDEXK_9 domains in other prokaryotes and eukaryotes. We searched the Pfam database for protein sequences containing the two domains and exhibited similar architecture with Cardinium homologs. In most cases, DUF1703 containing genes occurred in low copy number per genome. Most species carried fewer than four copies whilst only 9.8% of the species contained 10 copies or more (Fig. 5), ranking cCpun among the species with the largest number of DUF1703 paralogs. Species with higher abundance of DUF1703 paralogs are scattered across the prokaryotic taxonomy suggesting that DUF1703 protein expansion has occurred on multiple occasions within bacteria.

The reason for the expansion of the DUF1703 gene family in cCpun and its putative functional role is yet unknown. It is notable that DUF1703 genes have been also identified in the Rickettsia endosymbiont infecting biting midges (Pilgrim et al., 2017). Mirroring the pattern for midge Cardinium, the midge Rickettsia genome also contains multiple DUF1703 paralogs compared to other Rickettsia species with evidence of intragenic recombination (p < 0.001, PHI test, 1,000 permutations). However, Cardinium cCpun and Rickettsia DUF1703 carrying genes are phylogenetically unrelated (Fig. 4A) suggesting independent evolutionary histories, and independent expansion of this gene family in the two groups of midge symbionts. These data suggest this gene family may have a particular function in symbiosis with midges.

The biological role of the DUF1703 is still unclear. A recent transcriptomic study of the Cardinium strain cEper1 in its host E. suzannae showed that its only DUF1703 gene homolog is moderately transcribed in both sexes (Mann et al., 2017). Notably, a putative signal peptide cleavage site was predicted for 10 out of 25 DUF1703 paralogs in cCpun (Table S5) suggesting that they are potentially secreted, acting against DNA/RNA outside of the symbiont. Surprisingly, no signal peptides were detected in any of the paralogs identified in cSfur and cPpe (data not shown). It is noteworthy that an intact DUF1703 homolog of bacterial origin has been previously reported as component of the Maternal-Effect Dominant Embryonic Arrest (“MEDEA”) factor, a selfish genetic element reported in Tribolium castaneum (Lorenzen et al., 2008). PD-(D/E)XK nucleases constitute a large and functionally diverse superfamily of proteins which includes among others restriction endonucleases, Holliday junction resolvases, transposases, and DNA repair enzymes (Steczkiewicz et al., 2012). Recently, dual PD-(D/E)XK nuclease domains have been identified in a wide range of toxins from diverse intracellular bacteria (Gillespie et al., 2018; Lindsey et al., 2018). More interestingly, some of these domains have been directly linked with the induction of reproductive parasitism in the form of CI in Wolbachia (Beckmann, Ronau & Hochstrasser, 2017). Structural comparison of the PD-(D/E)XK core nuclease site from Cardinium and Rickettsia DUF1703 homologs and that of the CI-like toxins show considerable differences, especially in the sequence between the catalytic residues (Asp, Glu, and Lys) (Fig. S6). In addition, the AAA-ATPase domain associated with the DUF1703 nuclease is not found in the CI-like toxins of Wolbachia and related proteins (Gillespie et al., 2018) which might suggest these proteins have different functions. The biological role of Cardinium DUF1703 proteins remains to be determined.

Putative horizontal gene transfers as a source of genes in the cCpun genome

Horizontal gene transfer has been previously reported as the source of several genes in A. asiaticus, cEper1, and cBtQ1 (Penz et al., 2012; Santos-Garcia et al., 2014; Schmitz-Esser et al., 2010). Many of the HGT genes were found to be shared with members of the Alphaproteobacteria that have an intracellular lifestyle, especially species within the Rickettsiales order, consistent with HGT within the shared environment of the cell.

In accordance with previous observations of symbiont genomes, our results indicate that HGT has likely shaped the accessory genomes of cCpun (Table 2). The majority of the accessory genes of cCpun for which homologs could be assigned in the database are more similar to corresponding genes of bacterial species outside Bacteroidetes, with a bias to genes within the Proteobacteria having closest sequence similarity (Table 2; Fig. S7). For cCpun-specific genes, closest sequence matches lay within bacterial species known to be associated with other arthropods including Rickettsia and Wolbachia, as well as the amoeba-associated bacteria Candidatus Paracaedibacter acanthamoebae and Candidatus Jidaibacter acanthamoeba (Table 2). Four of these genes clustered with gene sequences from torix group Rickettsia, which are also found in midges. Three of these genes encode putative transposases, and one is a hypothetical protein that in other Rickettsia is located on a plasmid hypothesized to be important in determining the host-symbiont interaction (Gillespie et al., 2015).

Among the putatively horizontally exchanged gene set were ORFs encoding a carbonic anydrase (CA), an amino acid permease, and a putative chromosome-partitioning protein. Finally, two cCpun-specific genes encoding hypothetical proteins had their closest homologs within Aedes mosquitoes (Table 2; Fig. S7). Notably, the two proteins also have partial similarities with a large ankyrin repeat containing protein (Aasi_1610) previously identified in A. asiaticus (Schmitz-Esser et al., 2010). Although both cCpun proteins had their ten top hits assigned to Aedes sequences, the partial similarities to A. asiaticus suggest that they might be fragments of an Aasi_1610 distant homolog. Note, the number of these genes derived from HGT may be even higher since the majority of the accessory genes did not have any significant matches on the GenBank database, and many of these likely represent HGT events from as yet uncharacterized genomes.

The presence of CAs gene is interesting. Among the Amoebophilaceae, CA homologs were detected only in cCpun, cSfur, and cHgTN10 and not in other Cardinium strains nor A. asiaticus, Notably, the three Cardinium homologs do not form a monophyletic group, with cHgTN10 and cSfur homologs being clustered together and more closely associated with a putative CA previously identified in the Rickettsia endosymbiont previously found in biting midges (Pilgrim et al., 2017) (Fig. S8). Our results suggest that the Cardinium CA homologs have independent evolutionary histories and probably originated from independent horizontal transfer events into the three genomes.

The function of these CAs is not clear. CAs are ancient and ubiquitous multi-class zinc-containing metalloenzymes that catalyze the interconversion of CO₂ to bicarbonate (Smith & Ferry, 2000; Smith et al., 1999) and are involved in a variety of biochemical processes including respiration and pH homoeostasis (Gai et al., 2014). Studies have shown that CAs are essential for microbial growth in free living bacteria under ambient air with low levels of CO₂ (Mitsuhashi et al., 2003; Merlin et al., 2003; Kusian, Sültemeyer & Bowien, 2002). However, whilst CAs are common in many bacterial groups, they are less commonly observed in the genomes of obligate intracellular bacteria (Ueda, Nishida & Beppu, 2012). Studies suggest that intracellular pathogens may rely on CAs for virulence and survival within the host cell (Valdivia & Falkow, 1997), possibly through regulating the phagosome pH during the infection (Nishimori et al., 2014). An intriguing hypothesis is whether CAs might actually play a role in the survival of Cardinium outside of the host in comparable way to the role of CAs in free living bacteria, and thus facilitating its horizontal transmission. Interestingly, the plant-mediated horizontal transmission of Cardinium bacteria between phloem sap-feeding insects has been previously reported, supporting such a scenario (Gonella et al., 2015).

cCpun lacks a biotin or other B-vitamin biosynthetic pathways, indicating it is unlikely to act as a source of these vitamins to its haematophagous host. Indeed, putative homologs of the complete biotin transport system (BioY: CCPUN_01590, BioM: CCPUN_08370, and BioN: CCPUN_08380) were detected, suggesting that cCpun may depend on external provision of biotin from the host. The presence of a complete biotin transporter gene set contrasts with other Cardinium genomes, which lack these transporters, but may carry complete operons for the synthesis of biotin, lipoeta and pyridoxal 5’-phosphate (vitamin B6) (Penz et al., 2012). Exception is the recently sequenced strain cSfur which encode for both a biotin transport system and a complete operon for biotin synthesis (Zeng et al., 2018).