As part of a field survey seeking to uncover new nematode-associated endosymbionts, approximately 100–500 g of soil and roots were collected using a soil auger or serrated shovel, collecting from the top 15 cm of soil at a 30 cm to 1 m distance from the base of various plants at different sites (Supplementary Table 1). All handling and processing of samples were performed in compliance with USDA APHIS permits to A.M.V.B. Samples were kept cool (<10°C) until processed to isolate nematodes as described previously (Myers et al., 2021). Briefly, soil and roots were placed in Baermann funnels for 3–5 days, allowing nematodes to be collected. Nematodes were further separated from soil by sucrose flotation following (Jenkins, 1964) to remove remaining soil particles, and finally further washed to remove further contaminants by mixing in 40 mL distilled water, centrifuging, and discarding the supernatant.

To characterize nematode community and Wolbachia-like DNA, total nematode community DNA was isolated from each sample. After a brief examination under an inverted microscope, about 800 to 1,400 nematodes from each sample were exposed to five cycles of freeze-thaw to break cuticles. DNA was isolated using the Qiagen DNeasy Blood & Tissue Kit (Valencia, CA, United States) following the manufacturer’s directions. DNA quantity and quality were assessed on the Nanodrop spectrophotometer. An initial sample was prepared for sequencing and analysis without pre-screening as follows: approximately 0.6–1 μg of DNA was used for genomic library preparation with the QIAseq FX 96 DNA Library Kit (Valencia, CA, United States) following the manufacturer’s directions with fragmentation times and AMPure bead size selection steps optimized for 450–550 bp fragments. Library quality and quantity were assessed on the TapeStation 2200 (Agilent, USDA). Libraries were normalized and pooled before sequencing on Illumina HiSeq, with 150 PE cycles performed at Genewiz, Inc. (NJ). After positive results for a plant-parasitic nematode-associated Wolbachia strain were obtained from one sample (following bioinformatics analysis described below), subsequent samples from the same farm were pre-screened for the presence of Wolbachia prior to sequencing with custom primers designed based on published qPCR assay primers (Mee et al., 2015) that showed a high sensitivity for low-titer infections. However, as these primers had mismatches to 16S rRNA gene sequences from plant-parasitic nematode-associated Wolbachia strains, we modified the primer sequences to increase specificity to our targets using aligned sequences. The resulting new primers were: forward Wol-mee-F 5′-CTC ACA GAA AAA GTC CT-3′ and reverse Wol-mee-R 5′-CGC CTT TAC GCC CAA T-3′, with thermal cycle conditions: 95°C for 2 min, 35 cycles of 95°C for 30 s, 59°C for 30 s, 72°C for 1 min, and one cycle of 72°C for 10 min.

To recover nematode and symbiont genomes, reads were de novo assembled (see details in Supplementary Material), then Wolbachia-like contigs were annotated. First, reads for each sample were filtered and trimmed using Trimmomatic version 0.38 (Bolger et al., 2014), and overlaps in paired reads were detected and merged in Pear version 9.11 (Zhang et al., 2014). Filtered paired reads and merged reads were de novo assembled with metaSPAdes version 3.13 (Bankevich et al., 2012; Nurk et al., 2017) using low kmers (25, 33, and 45). Assembly quality was assessed using Quast version 5.0.1 (Gurevich et al., 2013). Assemblies were screened for Wolbachia-like 16S rRNA using a two-step analysis with blastn in Blast + version 2.10.1 (Camacho et al., 2009) (−evalue 10) first to a custom database of Wolbachia 16S rRNA sequences, and then a second blastn to the complete NCBI nt database. Any samples with top blastn hits to PPN Wolbachia from this read-based blast were also considered “positive.” For PPN Wolbachia-positive samples, full genomes were extracted using similar two-step blastn searches, first to Wolbachia genome databases, then to the full nt database.

Based on the high sequence similarity among samples, and their origin from the same farm, we combined these samples for further analysis, to improve coverage and assembly quality. We used an iterative map-assemble approach using bwa version 1.17 (Li and Durbin, 2009) and a subtractive mapping approach as follows. First, we identified non-Wolbachia contigs in the initial assembles using both blastn results and a%GC filter using prinseq-lite.pl in the BRAbB software (Brankovics et al., 2016), then we used bwa mem to map each sample’s reads to this non-Wolbachia data specific to our samples, then we used samtools version 1.9 (Li et al., 2009) and custom scripts to extract unmapped (i.e., Wolbachia-enriched) reads. These enriched reads were concatenated for all samples and assembled in metaSPAdes with a kmers (25, 45, 65, and 99). The resulting assemblies were combined with Wolbachia-like contigs from the individual samples with the new strain. The process of bwa-based subtractive mapping was repeated from the original reads with this new, improved database. The new Wolbachia-enriched reads were again concatenated and de novo assembled again with metaSPAdes. This approach was repeated three times and stopped once the sum of the length of the resulting new strain Wolbachia-like contigs ceased to increase between cycles and inspecting genome contamination and completeness metrics using CheckM v1.0.18 (Parks et al., 2015) at intermediate steps.

The resulting contigs were assessed by several quality controls to reduce the likelihood of spurious bioinformatic contamination with non-Wolbachia data. To filter out contigs with potential short horizontally transferred Wolbachia-like DNA regions (HGTs), long contigs (>5,000 bp) were removed if coverage was >2 times the average coverage of the longest contigs, using coverage analysis in the pileup.sh in BBMap version 38.9 (Bushnell, 2014). To filter out possible HGTs, contigs were imported into Geneious Prime version 2020.0.4 (Biomatters, Ltd.) and inspected, with contigs >1,000 bp removed if GC content was below 24% or above 42%. Quality was assessed by annotating contigs using Prokka version 1.14.6 (Seemann, 2014) which uses Prodigal for ab initio gene prediction, HMMER3 for protein family profiles, BLAST+ for comparative annotation, Barrnap¹ for rRNAs, Aragorn (Laslett and Canback, 2004) for tRNAs. The resulting genes were then analyzed with blastn to the nt database and with DIAMOND blastx to the full nr database to keep only contigs with the highest similarity to PPN Wolbachia. Finally, to check for regions of possible bifurcating misassembly due to mutational differences in the field sampled specimens, contigs were aligned and checked for blocks of near identity and synteny using Geneious Prime plugins ProgressiveMauve v1.1.1 (Darling et al., 2010) and LASTZ alignment v7.0.2 (Biomatters, Ltd.).

The final assembly quality was analyzed in Quast version 5.0.1 (Gurevich et al., 2013) and re-annotated in Prokka. Completeness and contamination were assessed in CheckM and by evaluating the genome presence of housekeeping genes and tRNAs. Pathways were assessed by gene presence-absence comparisons with other Wolbachia and outgroups Anaplasmataceae using Roary version 3.13 (Page et al., 2015), and using ModelSEED version 2.6.1 (Seaver et al., 2021) which assesses metabolic models in the KEGG and MetaCyc databases. ModelSEED was run with “complete” supplemented in silico media.

To assess possible nematode hosts for the Wolbachia from our sampled nematode communities, we compared 21 sampled nematode communities (Supplementary Table 1) that were extracted and sequenced as described above, as part of a broader study. All partial cytochrome oxidase I (COI) sequences with top blastn match to nematodes were extracted from our scaffolds using blastn and custom scripts. The relative abundance of nematodes in each sampled community was then estimated based on kmer coverage with the equation C = (C_K.R)/(R − K + 1), where C is total coverage, C_K is kmer coverage, K is the length of kmers, and R is read length and normalization to the total number of reads sequenced for each sample. To this matrix of COI hits, we added a column with the normalized coverage of Wolbachia 16S rRNA for the Wolbachia-positive samples. We then calculated and plotted Spearman rank correlation (rho) values and p-values using the R package Hmisc version 4.5-0 ‘rcorr’ which calculates a matrix of Spearman’s rho rank correlation coefficients for all pairs of columns for non-missing elements, using midranks for ties, with method “spearman,” order “hclust,” hclust.method “average” and plotting using ‘corrplot’. Multiple testing correction was performed using the Benjamini and Hochberg (1995) method in R with p.adjust “BH.”

To survey publicly available sequencing projects for potential PPN Wolbachia-like sequences, we developed scripts that take input queries such as geographic locality, sample type, sequencing strategy, and access the NCBI Sequence Read Archive (SRA) data, then analyze downloaded datasets to screen for Wolbachia. Keywords include “rhizosphere,” “soil,” “grassland,” “forest,” “agriculture,” “nematode,” and strategies/platforms such as “MiSeq,” “HiSeq,” “Illumina,” “Amplicon,” or “wgs,” and combinations of these terms. Briefly, these scripts use the E-utilities public API from NCBI to obtain and analyze serially tabular sets of raw read data for SRA projects via concatenated SRR, ERR, and DRR run data, subjecting them to raw read 16S rRNA gene two-step blastn, as described above, using for the custom blastn all PPN Wolbachia, including the new Wolbachia 16S rRNA gene from the present study. For SRA experiments with positive results, reads were further assessed after adapter and quality trimming and filtering using Trimmomatic and merging using Pear as described above. The resulting Wolbachia-like sequences were then analyzed using phylogenetic analysis approaches described below.

To understand the evolutionary relationships among candidate PPN hosts, candidate PPN-type Wolbachia strains from SRA data mining, and the new Wolbachia strain compared to other strains, various phylogenetic analyses were performed. In most cases, the general approach involved first blastn searches of our sequences of interest to the nt database at NCBI to find ∼100–500 most closely related sequences. We then download and aligned these using either Mafft version 1.0.4 (Katoh and Standley, 2013) or Clustal Omega 1.2.3 (Sievers et al., 2011), and trimming and removing duplicates or highly similar sequences within the Geneious Prime version 2020.0.4 (Biomatters, Ltd.) suite, prior to phylogenetic analysis using both maximum likelihood (ML) phylogeny reconstruction was performed in RAxML version 4 (Stamatakis, 2014), assessing bootstrap support from 500 replicates and Bayesian phylogeny estimation with MrBayes version 2.2.4 (Huelsenbeck and Ronquist, 2001; Ronquist et al., 2012) with final phylogenetic trees visualized in FigTree version 1.4.4² with labels and color added in Adobe Illustrator. Specific phylogenetic analysis approaches are as follows.

For inferring relationships among candidate PPN hosts predicted from community correlation analyses (described above), partial nematode cytochrome oxidase 1 (COI) genes that were significantly correlated with Wolbachia-positive samples were aligned and analyzed with the ML GTR Gamma nucleotide model, with rate heterogeneity alpha estimated, and with rapid bootstrapping and search for the best-scoring ML tree (-f a -x 1) and Bayesian analysis with the GTR + G model with 4 categories, and Markov chain Monte Carlo settings of chain length 1,100,000, 4 heated chains, heated chain temp 0.2, subsampling frequency 200, Burn-in length 100,000, with random seed 31,569, and priors with unconstrained branch lengths GammaDir (1,0.1,1,1), checking for convergence with minESS > 200.

For inferring relationships among candidate PPN-type Wolbachia strains from SRA data mining, an initial set of thousands of trimmed and merged reads with top blastn similarity to PPN-type Wolbachia were aligned to reference 16S rRNA genes, then identical sequences were removed. Preliminary phylogenetic analyses in FastTree v2.1.11 (Price et al., 2010) with the GTR model were used to generate a preliminary tree to remove large numbers of sequences from the alignment that displayed exceptionally long branches, similar to the distance separating Ehrlichia/Anaplasma and Wolbachia clades, on the basis that these sequences may reflect either non-Wolbachia alphaproteobacterial, or possibly degrading 16S rRNA pseudogenes or degrading horizontally transferred 16S rRNA fragments. The resulting subset of SRA sequences was analyzed in three separate alignments for the sub-region of the 16S rRNA gene in which they occur, and then together in a larger alignment to include all sequences together. RaxML and MrBayes phylogenies were reconstructed as described above.

For inferring relationships among the new Wolbachia strain and other Wolbachia strains, we focused on three sets of analyses: two which included strain wRad which was the only other strain to date from a PPN, and one which included a larger number of genes even though wRad could not be included. For strain wRad, only 3 gene regions have been sequenced. The first analysis was able to include additional early-branching Wolbachia isolates for which only the 16S rRNA gene sequence was available, and this analysis included additional outgroups Candidatus Neowolbachia serbia and Candidatus Mesenet longicola (GenBank Accessions MH618374.1 and BNGT01000041.1, respectively) and related bacteria from Harpalus pennsylvanicus, as well as additional Rickettsiales and outgroup alphaproteobacterial. The second analysis included more gene regions available for strain wRad, including 16S rRNA, partial CTP synthase and ftsZ, and partial groES and groEL (GenBank Accessions EU833482.1, EU833483.1 EU833484.1, respectively) and corresponding gene regions for other Wolbachia and outgroups with gene regions with blastn matches to these from the NCBI nt and genome databases. Once aligned, the sequence block was concatenated for strains with all regions covered, and the final alignment was stripped of sites with gaps or ambiguities. We produced additional reduced-length alignments from this full gap-stripped alignment to test for potential long-branch attraction artifacts. These reduced alignments were generated by successively stripping sites with the highest evolutionary rates identified using TIGER v2.0 (Cummins and McInerney, 2011). Phylogenetic analyses on resulting blocks were performed as described above.

For a more robust phylogenomic analysis of PPN-type Wolbachia with their outgroups, we prepared an alignment of a larger number of core genes from full Wolbachia and outgroup genomes downloaded from NCBI. We used Roary to generate a codon-based alignment (core_gene_alignment.aln) block for analysis. As for the previous analyses, to control for potential alignment artifacts, we removed all positions with ambiguities or gaps and performed TIGER analysis to create additional shorter alignments with high evolutionary rate sites progressively removed, to potentially reduce the effects of long-branch attraction. ML and Bayesian analyses were performed as described above for each alternative alignment. Finally, using the initial codon-based alignment, we performed in silico translation, removed gaps and ambiguities, and performed ML analysis of amino acid sequences was performed with RaxML using the PROTGAMMAGTR substitution model with empirical base frequencies and 500 bootstrap replicates.

To evaluate gene repertoire differences among Wolbachia strains and outgroups and analyze functional gene ontology (GO) enrichment, we first performed gene annotation on our new Wolbachia strain assembly as well as 200 other draft Anaplasmataceae genomes from NCBI using Prokka with identical parameters for all assemblies. We then created ortholog sets using Roary v3.13.0 (Page et al., 2015) which performed blastp on gff files from Prokka, with parameters -e for codon-aware alignment in PRANK (Löytynoja, 2014) and -i 60 to allow for detection of distant orthologs with outgroups. Initial Roary results suggested some draft genomes were potentially incomplete (too few genes) or contaminated (too many genes), so this genome set was reduced to the best genomes, comprising 93 draft genomes from Wolbachia strains plus 9 outgroup draft genomes. Roary gene clusters without clear gene annotations, listed as “group_#,” were re-examined and if possible assigned gene names by cross-referencing Prokka gene calls to gene and GO annotations for Wolbachia strains downloaded from UniProtKB database and other databases (MetaCyc, KEGG). These databases were also used to create a master GO annotation file for Anaplasmataceae for downstream enrichment analyses. Pangenome and core genome comparisons based on Roary gene_presence_absence.csv outputs were depicted with the online Venn drawing tool http://bioinformatics.psb.ugent.be/webtools/Venn/. Functional GO enrichment was assessed using topGO version 2.4.0 (Alexa and Rahnenfuhrer, 2020) which assesses GO-term graph topology ad creates test statistics using the algorithm ‘weight01.fisher’ which returns corrected p-values not affected by multiple testing. TopGO was implemented in R using the script aip_topgo_usage.consider_universe.R³ for multiple gene subsets (depicted in Venn diagrams) using ‘diff’ and ‘universe’ gene sets.

To assess genes and gene regions that may be important functionally, through signatures of increased or decreased purifying or directional selection, we applied an analysis approach that assessed pairwise dN/dS followed by GO enrichment tests on various high and low dN/dS gene sets. These analyses involved first generating new pairwise codon-based nucleotide alignments of orthologs generated in Roary, for each pair of Wolbachia strains. We then generated sliding window blocks using the function split.java in KaKs Calculator 2.0 (Zhang et al., 2006) to create 1,200 bp length windows with overlaps of 600 bp, preserving codon positions. We then used KaKs Calculator to assess dN/dS on all windows, specifying genetic code 11. This program assesses Ka/Ks (or dN/dS) while controlling for multiple substitutions per site and using a maximum likelihood framework for model selection and AICc for model averaging (Zhang et al., 2006). In addition to multiple substitution corrections performed by this software, we also set a maximum Ks cutoff for our output data of 2 for Wolbachia pairs in the major supergroups, and 2.5 for Wolbachia in the PPN supergroup. Annotated genes were matched and counted as within windows if they occurred across at least 300 bp of a given window. Output dN/dS values were grouped into subsets for topGO functional enrichment analysis, partitioning genes into sets with values into top 10%, top 25%, the bottom 10%, and bottom 25% of each pairwise comparison.

Screening for PPN-type Wolbachia was initially performed on 16 samples (Supplementary Table 1) with the initial screening of raw reads showing 10 out of 16 samples with reads mapping to the Wolbachia 16S rRNA gene. However, blastn analysis of these reads showed that most of these read hits were more similar to non-PPN Wolbachia strains than to PPN strains, with only one sample (P3-11) showing high similarity to wPpe from the plant-parasitic nematode P. penetrans and wRad Wolbachia from R. similis. This sample was from a fruit tree farm in Los Fresnos, Texas (26.1585 N 97.3844 W) consisting of pooled soils from the following fruit trees: mango (Mangifera indica; Sapindales: Anacardiaceae), guava (Psidium guajava; Myrtales: Myrtaceae), pomello (Citrus maxima; Sapindales: Rutaceae), sugar-apple (Annona squamosa; Magnoliales: Annonaceae), and sapodilla (Manilkara zapota; Ericales: Sapotaceae). After the discovery of this initial positive result, separate samples were collected from 14 individual fruit trees at the same farm, from which the following five fruit tree samples were found to be PCR-positive and were processed for community shotgun sequencing: sugar-apple, avocado (Persea americana; Laurales: Lauraceae), plantain (Musa × paradisiaca AAB; Zingiberales: Musaceae), guava, and “Rosetta fruit” (Malus pumila “Niedzwetzkyana”; Rosales: Rosaceae) (Supplementary Table 1).

Based on initial alignments of contigs from six Wolbachia-positive samples which originated from the same farm, inter-sample divergence was found to be low, therefore, to increase coverage, final genome assembly was performed for the six samples with reads pooled together. The final assembly (NCBI GenBank accession JAIXMJ000000000) of the new Wolbachia strain, hereafter, designated wTex (named for its location in southern Texas), consisted of 192 scaffolds with a total length of 1,013,022 bp, maximum scaffold length of 57,862 bp, N50 of 10,082, with 33.49% GC, coding density 0.809, with 989 predicted genes, and a full set of rRNA and tRNA genes (3 and 38, respectively), and average coverage of 15.96X (Supplementary Table 2). CheckM-based genome completeness was 84.37% based on 368 markers from 63 Rickettsiales genomes, with a contamination score of 0.64 CheckM. Prior to this final assembly, during iterative steps to improve the assembly CheckM completeness scores decreased as contamination scores decreased. Comparative genome features across similar Wolbachia strains are shown in Table 1. Among the predicted genes, 398 (40.2%) had no known function. GC content (Figure 1) was intermediate between that of clades C and D Wolbachia strains from filarial nematodes and clades A and B from arthropods, whereas assembly length was shorter than that of the majority of clades A and B Wolbachia strains and was more similar to that of clades C and D Wolbachia.

Correlation analyses (Supplementary Figure 1) showed Wolbachia wTex was positively associated with nematode hits to the partial COI gene of Heteroderidae sp. CD2526 (GenBank accession MK033155.1) (rho 0.659044735, BH-corrected p-value 0.02339631) and Ptycholaimellus sp. M1 (GenBank accession KX951909.1) (rho 0.713184097, BH-corrected p-value 0.006335523), however, hits matching the latter were very short and these scaffolds were at very low coverage. Scaffolds identified as matches to the COI of Heteroderidae sp. CD2526 were generally longer > than 1,000 bp with longer blast matches and were therefore extracted for further analysis. This produced 9 scaffolds with similarities to this COI hit. Similar sequences from GenBank were downloaded and aligned and analyzed by ML and Bayesian methods, which showed 7 sequences clustered with support in a clade together with Helicotylenchus spp. while 2 sequences clustered in a sister-clade with Rotylenchus spp. (Supplementary Figures 2, 3).

Phylogenetic analyses for the 16S rRNA gene alone (Supplementary Figure 4) which included the broadest array of outgroups and additional early-branching Wolbachia in this study and an alignment block length of 1,573 nucleotide positions, resulted in a strongly supported clade for PPN Wolbachia strains wTex, wPpe, and wRad, in a basal position in the Wolbachia phylogeny. The closest early-branding sister taxa to this PPN Wolbachia clade included a well-supported cluster from various aphids including Pentalonia nigronervosa wPni-like strains (accessions NZ_JACVWV01000005.1, KJ786949.1, KJ786950.1) along with strains from the conifer aphid Cinara cedri (AY620430.1) and the trunk-feeding aphid Stomaphis sinisalicis (KF751211.1), and strains wBta from the whitefly Bemisia tabaci (KF454771.1) and wBry from the spider mite Bryobia spec. V (EU499316.1). The largest sequence difference in the 16S rRNA region among PPN-type Wolbachia strains was 4.211%, whereas the average difference between the PPN-clade and wPni-like strains was 3.934%. A separate early-branching sister clade consisted of closely related isolates of a Wolbachia from the fungal-feeding mold mite, Tyrophagus putrescentiae (Supplementary Figure 4).

Phylogenetic analyses for the concatenation of 3 gene regions (16S rRNA, ftsZ, and groEL) to better resolve the phylogeny including all three PPN-type strains wTex, wPpe, and wRad resulted in a strongly supported clade for PPN Wolbachia strains (Figure 2). The supported clade was obtained for both ML and Bayesian analyses, which produced similar tree topologies, and for all alternative alignments including the full sequence alignment, the gap-stripped alignment, and alignments with progressive stripping of high evolutionary rate sites resulting in alignment lengths of 2,368 to 4,674 nucleotide positions. Bootstrap support for the PPN Wolbachia clade was 85–100%, with support of 100% for the majority of the alternate alignments. The PPN Wolbachia clade formed a strongly supported earliest branch for the genus Wolbachia (supergroup L in Figure 2) for alignments with full data, gap-stripped data, and data with high evolutionary rate sites stripped that included over 68% of the original alignment length, however, progressive stripping of such high-rate sites placed strain wPni at the root of the tree with or without >50% bootstrap support and reduced support for other basal branches.

Phylogenomic analyses of 100 genome-wide protein-coding genes (orthologs from Roary), resulted in similar results to the 3-region results (above), producing a strongly supported clade for PPN Wolbachia strains wTex and wPpe (Figure 3). ML bootstrap support and Bayesian posterior probabilities were largely similar for these analyses, showing support for this clade of 100% bootstrap and posterior of 1. These results were consistent for alternative alignments including the full sequence alignment, the gap-stripped alignment, an alignment with 3^rd codon positions removed, a translated alignment, and alignments with progressive stripping of high evolutionary rate sites resulting in alignment lengths of 17,382 to 87,628 nucleotide positions. The exception to this high support was alignments stripped to less than 20% of the original positions, in which there was very little support for any Wolbachia clade. As with the 3-region phylogenetic analysis, the 100-gene phylogenetic analyses produced a strongly supported earliest branch position for the PPN Wolbachia clade (supergroup L in Figure 3) in some, but not all alternative alignments. For example, for alignments with full data, full data with gaps and ambiguities stripped, and full data with 3rd codon positions stripped, there was 100% bootstrap support for the PPN Wolbachia clade forming the earliest branch for Wolbachia. However, for alignments with high evolutionary rate sites stripped such that the alignment was 25–43% of its original length, the earliest branch of genus Wolbachia varied, sometimes placing wPni at the root of the tree just before PPN Wolbachia, and sometimes placing strains wFol and wCfeT as sisters to PPN Wolbachia, together forming the earliest branch. However, one alignment with high-rate sites stripped, with 22% of the original length, produced 91% bootstrap support for PPN Wolbachia at the root of the genus. Phylogenies based on amino acid alignments, however, also placed strain wPni at the root with 82% bootstrap support.

Sequence read archive (SRA) database screening of 3,400 amplicon experiments from soils and rhizospheres produced 81 sequence runs (i.e., SRR/ERR/DRR files) with the highest blastn similarity to the 16S rRNA gene from PPN-type Wolbachia wTex, wPpe, or wRad. Following read-merging, there was 4,535 top sequences read blastn matches to PPN-type Wolbachia among these runs. After removal of identical sequence reads and sequences leading to exceptionally long branches in phylogenies and sequences placed ambiguously between the Ehrlichia/Anaplasma and Wolbachia clades, there were 61 unique sequences unambiguously grouped with the Wolbachia clade. These sequences were from 24 separate SRA runs originating from the United States, France, Germany, Sweden, Switzerland, Japan, India, and Malaysia. Although phylogenetic analyses with these short sequences produced generally low bootstrap support and low posterior probabilities, tree topologies suggest several clades (Figure 4). Group 1 comprised 19 sequences clustered at the root of the Wolbachia tree along with Wolbachia from PPNs with broad geographic origins (France, Germany, Sweden, Switzerland, Japan, Malaysia, and in the United States, from Florida, Michigan, California, Appalachia) from diverse ecosystems. Groups 2 and 3, from Malaysia and India, clustered near the root of the Wolbachia CDF supergroups. Groups 4 through 8, comprising 28 sequences, were all from a beech forest in Germany and sequences clustered mostly with various Wolbachia strains from Collembola. Group 9 was a distinct cluster of 6 sequences from the same beech forest with similarity to Wolbachia from quill mites, while Group 10 formed a cluster of 6 sequences allied with Wolbachia from Curculio sp. (weevils) from Minnesota. Phylogenetic trees for separate sub-regions produced similar results (Supplementary Figures 5–7).

Specific predicted pathways and genes of interest based on known or hypothesized functions were searched for in the wTex assembly. The wTex genome was similar to wPpe and other Wolbachia in having conserved pathways for glycolysis and the tricarboxylic acid cycle and pathways for biosynthesis of nucleotides including the pentose phosphate pathway, and peptidoglycan and fatty acids, but lacking genes for key steps or most steps of other biosynthetic processes including amino acid, vitamin and co-factor, and carbohydrate synthesis, suggesting incomplete pathways. Both PPN-type Wolbachia strains wTex and wPpe had full-length predicted genes for heme synthesis that were also conserved across outgroups and other Wolbachia, including hemA (encoding 5-aminolevulinate synthase EC 2.3.1.37), hemB (encoding delta-aminolevulinic acid dehydratase EC 4.2.1.24), hemC (encoding porphobilinogen deaminase EC 2.5.1.61), hemE (encoding uroporphyrinogen decarboxylase EC 4.1.1.37), ctaB (encoding protoheme IX farnesyltransferase EC 2.5.1.141), hemF (encoding oxygen-dependent coproporphyrinogen-III oxidase), and hemH (encoding ferrochelatase EC 1.3.3.3).

Riboflavin synthesis and transport genes were notably absent in wTex and wPpe, except for ribB (encoding 3,4-dihydroxy-2-butanone 4-phosphate synthase EC 4.1.99.12) and one outstanding gene annotated as ribN (riboflavin transporter) in wPpe that was not present in any Wolbachia or outgroup strains. Similarly, biotin synthesis genes were absent in wTex and wPpe, with no homologs found matching either outgroup Anaplasmataceae-type biotin genes or the frequently horizontally transferred ‘BOOM’ (biotin synthesis operon of obligate intracellular microbes) operon genes, bioA (encoding adenosylmethionine-8-amino-7-oxononanoate aminotransferase EC 2.6.1.62), bioB (encoding biotin synthase EC 2.8.1.6), bioC (encoding malonyl-[acyl-carrier protein] O-methyltransferase EC 2.1.1.197), bioD (encoding ATP-dependent dethiobiotin synthetase EC 2.6.1.62), bioF (encoding 8-amino-7-oxononanoate synthase 2 EC 2.3.1.47), and bioH (encoding pimeloyl-[acyl-carrier protein] methyl ester esterase EC 3.1.1.85). However, both wTex and wPpe assemblies had a predicted bioY gene (encoding a biotin importing transporter protein), and birA (encoding the biotin ligation protein bifunctional ligase/repressor EC 6.3.4.15), and diverged variants of the biotin utilization genes pccB (encoding propionyl-CoA carboxylase beta chain EC 2.1.3.-) and fabD (encoding malonyl CoA-acyl carrier protein transacylase EC 2.3.1.39). The latter two genes had such low sequence similarity to other Wolbachia strains and outgroups that they were not clustered as homologs to similarly annotated copies in Roary.

Additional predicted genes in wTex that were shared between PPN-type Wolbachia wTex and wPpe, but unique to these strains (i.e., not found in other Wolbachia) included a 1,350 nt gene lysC (encoding lysine-sensitive aspartokinase 3 EC 2.7.2.4, involved in the first step of the lysine biosynthesis via diaminopimelate ‘DAP’ pathway (see Figure 5) which also leads to methionine biosynthesis via de novo pathway, and threonine biosynthesis). This gene had the closest blastn and blastx similarity to genes from the alphaproteobacteria Candidatus Midichloria mitochondrii (hereafter Midichloria) (accession NC_015722.1), or distant genera of Vibrionales (Photobacterium, Vibrio), but no homologs to the standard Wolbachia or Anaplasmataceae variants of the lysC gene. Directly adjacent to this predicted lysC gene was a unique 987 nt PPN-type Wolbachia variant of asd2 (encoding aspartate-semialdehyde dehydrogenase 2 EC 1.2.1.11, which catalyzes the second step in lysine biosynthesis via the DAP pathway), with closest blast matches to genes from Midichloria and next to distantly related bacteria Photobacterium, and Vibrio. Although wPpe had both variants of asd2, wTex had only the Midichloria-like variant. These Midichloria-like asd2 and lysC genes were located adjacent to Wolbachia-like carA (encoding carbamoyl-phosphate synthase small chain EC 6.3.5.5, the first step of pyrimidine and arginine synthesis) and dapA (encoding 4-hydroxy-tetrahydrodipicolinate synthase EC 4.3.3.7, catalyzing the third step in the lysine biosynthesis pathway). The asd2 and lysC gene operon was syntenic, with the same gene order and orientation as that of Midichloria, although the latter had a branched-chain amino acid transaminase (BCAT EC 2.6.1.42) in the place of the dapA in PPN-type Wolbachia, whereas the distantly related versions of asd2 and lysC in other Wolbachia strains and outgroup Anaplasmataceae were not located together in tandem (Supplementary Figure 8). Phylogenetic analyses of all asd2 and lysC variants showed that the PPN-type Wolbachia and Midichloria versions of these genes have phylogenetically diverged from all other Anaplasmataceae/Wolbachia variants with high bootstrap support for a sister-clade to Midichloria for both genes (Supplementary Figures 9, 10). Nucleotide identity between Midichloria-like and PPN-type asd2 and lysC genes was ∼63 and 60%, respectively, while identity between wTex and wPpe homologs were ∼79 and 81%, respectively. Given the finding of these lysine biosynthesis pathway genes (Figure 5), and the observation that the missing final gene required for lysine synthesis lysA (encoding diaminopimelate decarboxylase EC 4.1.1.20) was previously predicted as a bacteria-to-eukaryote HGT in whiteflies (Luan et al., 2015) we searched the complete metagenome assemblies comprising mostly nematode contigs using blastn and blastp to see if we could detect lysA in a contig matching nematodes, but no such match was found. However, coverage was low, leaving some uncertainty about missing genes.

Another predicted lysine metabolism gene, shared only among PPN-type Wolbachia (wTex_00187) and wPni, was most similar (67% amino acid identity) to a saccharopine dehydrogenase (SDH) gene (EC 1.5.1.9) in the conifer aphid (Cinara cedri), involved in lysine catabolism (Figure 5). There was also synteny conservation in wTex and wPpe for a 1,308 nt predicted gene (wTex_00248) which had no blastn or blastx matches in the nt/nr databases, which was adjacent to hemA, which catalyzes the first step of heme synthesis. Another unique gene in wTex and wPpe was a 1,494 nt predicted gene tlcA (encoding ADP, ATP carrier protein 1) with its closest match being 70% nucleotide similarity to homologs in Rickettsia. There was a 729 nt predicted gene with partial blastn similarity to a squalene/phytoene synthase/isoprenoid synthase gene (potentially involved in carotenoid synthesis) from the conifer aphid C. cedri and a Wolbachia strain from the gall mite Fragariocoptes setiger. The latter was adjacent to dxr (encoding 1-deoxy-D-xylulose 5-phosphate reductoisomerase EC 1.1.1.267). Two other shared genes predicted in wTex and wPpe were arrayed in synteny: a 663 nt gene annotated as deoC (a gene involved in carbohydrate degradation EC 4.1.2.4) with high blastx similarity to BON (bacterial OsmY and nodulation) domain-containing protein genes in unrelated bacteria such as Flavobacteria, Holosporaceae, Cardinium, and Rickettsia; and an adjacent 246 nt predicted gene with DUF2188 domain and no matches in either blastn or blastx to the nt/nr databases. The first of these two genes had partial length homology (69–80% amino acid identity) to genes from two Wolbachia strains, wCfeT from cat fleas and a strain from the gall mite F. setiger. There were also predicted def (encoding peptide deformylase EC 3.5.1.88) and pterin-4-alpha-carbinolamine dehydratase genes (PCBD1 EC 4.2.1.96) involved in phenylalanine metabolism to tyrosine through tetrahydrobiopterin with distinct homologs in PPN-type Wolbachia that more closely matched non-Wolbachia, including eukaryotes (Culicoides midges) and distantly related bacteria (Francisella), or Rickettsia and an ameba endosymbiont (Candidatus Nucleicultrix amoebiphilia), respectively.

We searched for cytoplasmic incompatibility factor genes, cifA and cifB. From the 92 Wolbachia genomes analyzed in Roary, there were 152 cifA and cifB-like genes identified, yet these shared no homology with any predicted genes in wTex or wPpe or clade C or D Wolbachia. Next, we searched for WO prophage-like or plasmid-associated genes. Strain wTex showed no homology to the more than 2,000 phage or prophage-type genes detected in these 92 Wolbachia genomes analyzed. Nor did wTex have any detected homologs to the 84 plasmid-type genes or in this set of genomes. Other enriched features of PPN-type Wolbachia are described in the section below.

Analysis of gene repertoire overlap between PPN-type Wolbachia strains wTex and wPpe showed 501 predicted shared genes between these strains with an additional 400 and 437 genes only found in wTex and wPpe, respectively. Compared to wPpe, strain wTex was enriched for GO processes DNA-transposition, thiamine biosynthesis, and thiamine diphosphate biosynthesis (Supplementary Table 3), whereas the wPpe strain was enriched for GO processes mismatch repair and one-carbon metabolism (Supplementary Table 4). The thiamine enrichment in wTex arises from three genes, thiE (encoding thiamine-phosphate synthase EC 2.5.1.3), thiM (encoding hydroxyethylthiazole kinase EC 2.7.1.50), and thiD (encoding hydroxymethylpyrimidine/phosphomethylpyrimidine kinase EC 2.7.1.49 EC 2.7.4.7) which occurred as duplication in two scaffolds. Of these genes, only thiM was present in wPpe, whereas the variants of these genes in wTex had close matches only to one Wolbachia strain (wCfeT) and otherwise were most similar to these genes in the spirochete Brachyspira.

Comparison of the shared predicted genes between these two strains (i.e., core genes) compared to the four most closely related strains with complete genome assemblies (namely, wPni from the banana aphid Pentalonia nigronervosa, wFol from the springtail Folsomia candida, wCfeT from the cat flea Ctenocephalides felis, wChem from the bedbug Cimex hemipterus), not counting the 290 universally shared genes among these strains, showed the core PPN-type Wolbachia gene repertoire was most similar to that of wPni, followed by wFol, then wCfeT, and wChem (Figure 6A). This mirrored their phylogenetic places described above, with the four strains having 132, 104, 86, and 66 shared genes, respectively, not counting the 290 universally shared genes between all these strains. GO enrichment for the set of shared core PPN-type Wolbachia genes not shared with these other strains or shared with at most two other strains (Figure 6A) showed enrichment for biological processes diaminopimelate, pseudouridine, and chorismite biosynthesis, gluconeogenesis, and lysine biosynthesis via diaminopimelate, while GO enriched metabolic processes were metallo-aminopeptidase activity, lyase activity, monooxygenase activity, and ATP binding (Supplementary Table 5).

Comparison of genome repertoires of PPN-type Wolbachia strains wTex and wPpe to that of the pangenomes of remaining major supergroup clusters (supergroups C, D, and F denoted “CDFpan” and supergroups A and B denoted “ABpan”) (Figure 6B) showed more core PPN Wolbachia genes shared with CDFpan than with ABpan (683 vs. 649, respectively). Specifically considering genes shared with only CDFpan or ABpan, there was a similar pattern (i.e., 39 for CDFpan vs. 12 for ABpan). The AB pangenome was larger than the CDF pangenome (7,848 vs. 4,892 genes, respectively), therefore alternatively, the differences could be calculated as proportions of each pangenome cluster shared with PPN-type Wolbachia. The latter comparison showed PPN-type Wolbachia shared 0.1396 of their pangenome with CDFpan whereas 0.0817 of their pangenome with ABpan. Conversely, there were more accessory genes (not shared between PPN-type Wolbachia strains) that were shared with ABpan than with CDFpan (22 vs. 11 for wTex and 29 vs. 25 for wPpe), although this difference was not found when controlling for the approximately doubled number of accessory (non-shared) genes in ABpan compared with CDFpan (6,062 vs. 2,994 genes, respectively).

The GO enrichment analyses for PPN-type Wolbachia strains wTex and wPpe compared to other Wolbachia showed several significantly enriched functions (Figure 7 and Supplementary Tables 6–8). GO enrichment was compared for three sets of overlap (Figure 7), first for PPN-type Wolbachia compared to the pangenome of all Wolbachia, including genes shared between pangenomes (Figure 7A and Supplementary Table 6), then for the core shared genes from PPN-type Wolbachia compared to the pangenome of all Wolbachia (Figure 7B and Supplementary Table 7), then for the total pangenomes of PPN-type Wolbachia compared to the pangenome of all Wolbachia (Figure 7C and Supplementary Table 8). Among PPN-type Wolbachia, there was significant GO term enrichment (Figure 7A) for various cellular biosynthetic processes, including several nutrient pathways: protoporphyrinogen IX biosynthesis (part of heme synthesis), thiamine and thiamine diphosphate biosynthesis, lysine biosynthesis via diaminopimelate and diaminopimelate biosynthesis, fatty acid biosynthesis, cellular amino acid biosynthesis. Among GO enrichment for core functions shared by wTex and wPpe (Figure 7B) were many GO terms including nutrient pathways such as protoporphyrinogen IX biosynthesis, lysine biosynthesis via diaminopimelate and diaminopimelate biosynthesis, and fatty acid biosynthesis. For pangenomes of wTex and wPpe excluding shared genes with other Wolbachia (Figure 7C), enriched GO terms included thiamine and thiamine diphosphate biosynthesis, cellular amino acid biosynthesis, and cobalamin biosynthesis.

Analysis of dN/dS to investigate signatures of purifying or positive selection within PPN-type Wolbachia (wTex and wPpe) produced a range of dN/dS values from 0.0119336 to 0.334742, with a mean dN/dS of 0.08011301, considering values with Ks < 2. For genes within the dN/dS alignment block that fell within the top or bottom 10% or 25%, topGO enrichment analyses were performed against the pangenome of PPN-type Wolbachia (Supplementary Figure 11 and Supplementary Tables 9–12). Enriched GO categories for the highest values of dN/dS indicative of lower-than-average purifying selection included transcription antitermination, protein folding and transport, and various nutrient metabolism processes including arginine metabolism, N2-acetyl-L-ornithine:2-oxoglutarate 5-aminotransferase activity (a part of arginine metabolism), diaminopimelate biosynthesis (a lysine precursor), lysine biosynthesis, pyridoxal phosphate (vitamin B6) binding (Supplementary Tables 9, 11). Notable GO categories with the lowest dN/dS values (indicative of highest purifying selection) included translation and key metabolic functions associated with an energy (e.g., ATP binding and TCA) and metal ion, iron-sulfur cluster, and heme-binding (Supplementary Tables 10, 12).

The analysis of enriched GO categories for gene sets with highest and lowest dN/dS values associated with transitions between early-branching clades of Wolbachia was analyzed (Figure 8 and Supplementary Tables 13–20). Specifically, universally conserved lowest dN/dS (bottom 10% values) shared among early-branching clades of Wolbachia (strains wTex, wPpe, wPni, wFol, and wCfeT – branch topology and branch lengths extracted in Figure 8 from Figure 3) showed 7 GO categories including DNA topoisomerase type II activity and respirasome. In contrast, conserved lowest dN/dS gene sets showed diverse differences in GO terms unique to each branch. PPN-type Wolbachia strains wTex and wPpe showed unique (not shared) enrichment for heme binding, phosphorelay signal transduction system, ribosome, and zinc ion binding. Uniquely enriched GO terms from the lowest dN/dS gene set in the branch from PPN-type Wolbachia to wPni included 12 terms including protoporphyrinogen IX biosynthetic process. Uniquely low dN/dS gene set GO enrichment in the branch from wPni to wFol included four-way junction helicase and translation elongation factor activity, while the respective enriched terms from wFol to wCfeT were 4 iron, 4 sulfur cluster binding (Figure 8). Various GO categories were uniquely enriched for the highest dN/dS gene sets (top 10%) in these branches, such as for PPN-type Wolbachia, in terms of protein folding and transport, chaperone binding, and 3-dehydroquinate synthase activity (part of aromatic amino acid synthesis via the shikimate pathway) (Figure 8). For the branch from PPN-type Wolbachia to wPni, high dN/dS enrichment included glucosamine-1-phosphate N-acetyltransferase activity (involved in amino sugar metabolism) and many cellular membranes and cell shape/division terms. For the branch from wPni to wFol, high dN/dS enrichment included 4-hydroxy-tetrahydrodipicolinate synthase activity (involved in lysine biosynthesis) and various ribosome and rRNA processing functions, as well as shikimate metabolism (Figure 8). For the branch from wFol to wCfeT, high dN/dS enrichment included 7 terms (Figure 8), including thiamine pyrophosphate binding (involved in thiamine transport).