For 16S rRNA amplicon sequencing (detailed below), a total of 45 psyllid samples were analyzed, and included 44 different psyllid species. Two of the psyllid samples were the same psyllid species, Trioza urticae, however, the samples were collected from two different locations and time points (Supplementary Table 1). The 44 psyllid species represent 35 genera from five families within Psylloidea. Psyllid species were collected from 11 international geographic locations for this study and are detailed in Supplementary Table 1. Each psyllid sample consisted of an average of two psyllid adults pooled for DNA extractions, because most psyllid species are extremely small and not enough DNA can be obtained from just one individual per sample given our extraction protocol, which preserves voucher specimens for species identification of each sample (see below). Only the Trichochermes walkeri (P1–28) sample consisted of pooled nymphs (Supplementary Table 1). Psyllid DNAs were extracted following a nondestructive DNA extraction protocol (Percy et al., 2018) in order to correctly identify psyllid species with high confidence using both a molecular and morphological diagnostic approach as described in Percy et al. (2018). Briefly, collections were made by sweep net and aspirating directly from the plant, specimens were placed into 95–100% ethanol in the field and stored at −20°C thereafter. Field collected specimens were sorted and identified by DP and are maintained in DP’s personal collection; individuals selected for DNA extraction were then placed (whole or bisected) in the tissue lysate buffer with proteinase K [as per the Qiagen (Germantown, Maryland) protocol] and incubated for 12–15 h at 56°C. After incubation the specimen was retained and preserved in 70–95% ethanol as a DNA voucher. Genomic DNAs was obtained using a Qiagen DNeasy Blood and Tissue Kit and eluted in 200 μl of the elution buffer.

For library preparation of 16S rRNA libraries, custom 16S rRNA primers were synthesized by Integrated DNA Technologies, Inc. (San Diego, CA, United States) using a dual barcode design with 25 primer pairs based on the primer constructs used by the Earth Microbiome Project (EMP¹) similar to Walters et al. (2016) and Caporaso et al. (2012). Our aim was to amplify the 16S rRNA gene region of endosymbionts from psyllids including putative plant pathogens such as Ca. Liberibacter species. In turn, these barcode primers contained the primer sequences (341F and 805R) that target the 16S rRNA region. The primer sequences (341F and 805R) were designed previously by Morrow et al. (2017), because these primers successfully amplified psyllid endosymbionts in addition to plant pathogens such as Ca. Liberibacter species; the forward primer, 16S-341F = 5′-CCTACGGGNGGCWGCAG-3′, and the reverse primer, 16S-805R = 5′-GACTACHVGGGTATCTAATCC-3′. These latter primers span the V3–V4 region of the 16S rRNA gene resulting in a 464 bp amplicon. This primer set was selected, based on results from Morrow et al. (2017), because it did not preferentially amplify the most abundant obligate endosymbiont of psyllids, Carsonella. Carsonella is present in all psyllids at very high titers (Thao et al., 2000b). Moreover, these primers did not preferentially amplify plant plastid 16S rRNA genes from psyllid DNA (Morrow et al., 2017), which can abundantly amplify from psyllid DNA when targeting other variable regions of the 16S rRNA gene (Nachappa et al., 2011).

Amplicon PCR using the above primer set consisted of 4 μl of template DNA, 0.9 μl of each 10 μM primer, and 20 μl of Phusion master mix (ThermoFisher Scientific) per each reaction. The 16S rRNA amplicons were then purified using PureLink^TM PCR Purification Kit (Invitrogen). Using the cleaned PCR product as template, a second PCR was performed with HPLC-cleaned primers (5′-TATGGTAATTGTCCTACGGGNGGCWGCAG-3′ and 5′-AGTCAGCCAGCCGACTACHVGGGTATCTAATCC-3′) to complete the Illumina sequencing construct as described in McFrederick and Rehan (2016). Second PCR reactions consisted of 1 μl of template DNA, 0.5 μl of each forward and reverse primer, and 23 μl of Phusion master mix (ThermoFisher Scientific). For all PCR conditions, we followed PCR protocols as described in McFrederick and Rehan (2016). A total of 18 μl of PCR product from each reaction was normalized using SequalPrep^TM Normalization Plate Kit (Invitrogen). We pooled 4 μl of each of the normalized samples, and performed additional bead clean up. The quality of the libraries was then assessed using a 2100 Bioanalyzer.

All 16S rRNA amplicons from psyllid samples were multiplexed in a single lane for dual index sequencing on the Illumina MiSeq PE300 at the DNA Technologies and Expression Analysis Core Laboratory at the University of California, Davis. The mothur v.1.41.3 pipeline was used to process raw read data following the protocol described in the MiSeq Standard Operating Procedure to process high-quality contigs from overlapping paired-end reads (Schloss et al., 2009; Kozich et al., 2013). To trim off adaptors and control against sequencing errors, low-quality reads were trimmed and removed using the following parameters: the maximum ambiguous base = 0, the maximum sequence length = 440 bp. Chimera vsearch was then used to remove chimeric sequences. To assign taxonomy, the SILVA database (release 132) was used. Operational taxonomic units (OTUs) were clustered at the level of ≥ 97% nucleotide sequence similarity (Stackebrandt and Goebel, 1994).

Morrow et al. (2017) suspected a high amount of OTU splitting when using the primers 341F and 805R (see above) because several OTUs were highly similar in SILVA taxonomy affiliation within a single psyllid sample. These related OTUs potentially can be real and represent different co-occurring strains within a sample, however, they can also represent amplification and/or sequencing error (during library preparation and/or Illumina sequencing), or different 16S rRNA copies from a single bacterium. To help control further against amplification and sequencing errors, in addition to the read quality control filtering steps mentioned above, OTUs with low read abundance per sample were removed from the dataset during the OTU picking step if the total counts of contigs for a particular OTU were equal to and/or less than five. Moreover, to control against environmental or aerosolized contaminates during library preparation five control samples that consisted of water instead of the DNA template were also included during the construction of the 16S rRNA libraries to detect if library contamination occurred. After taxonomy assignments using SILVA (above) we used BLASTn searches against the NCBI non redundant nucleotide database to further annotate 16S rRNA contigs. Furthermore, the aim of this study was to identify common and abundant endosymbionts and plant pathogens within psyllid bodies, therefore, OTUs of mitochondria, chloroplasts, and Carsonella were removed prior to analysis.

The 16S rRNA gene sequences of known free-living bacterial relatives or endosymbionts belonging to three host association categories such as obligate, long-term symbionts that co-occur with obligate endosymbionts, and facultative endosymbionts were obtained from NCBI GenBank to compare AT/GC content with the same 16S rRNA region of sequences in this study (Supplementary Table 2). All NCBI sequences were aligned with the 16S sequences here to obtain the same 16S rRNA gene region using NCBI BLASTn. The percent GC nucleotide composition was calculated for each of these aligned 16S rRNA sequence regions using BBMap (Bushnell, 2014) and Genomics % G–C Content Calculator (Science Buddies 2021). The box and whisker plot were generated based on GC content data using the package “ggplot” (Wickham, 2016) in R v.4.0.2 integrated in Rstudio v.1.3.959 (R Core Team, 2021).

Psyllid samples where Ca. Liberibacter was identified, using 16S rRNA sequencing (above), were further examined using a metagenomic approach. The goal of metagenomic sequencing was to generate more sequence data for Ca. Liberibacter loci, to further determine the evolutionary relationships of newly identified Ca. Liberibacter taxa. Psyllid samples that were positive for Ca. Liberibacter, based on 16S rRNA sequencing, include Macrohomotoma gladiata (P1–9), Russelliana capsici (P1–16), and Trichochermes walkeri (P1–28). A preserved herbarium specimen of Solanum umbelliferum (HS9) collected in 2011 that was infected with the Ca. L. psyllaurous haplotype G, based on our previous study Mauck et al. (2019), was also sequenced. Prior to library preparation, each sample was subjected to whole genome amplification using the Repli-G whole genome amplification kit (Qiagen) according to the manufacturer’s instructions. Illumina library preparation and sequencing were conducted by Vincent J. Coates Genomics Sequencing Laboratory at the University of California, Berkeley. To increase the number of reads per sample, each Ca. Liberibacter-infected psyllid/plant library was sequenced twice as paired-end, 150 bp reads on two lanes of Illumina NovaSeq SP. The quality of raw reads for each sample run was verified with FASTQC v.0.11.3 (Andrews, 2010), and low-quality reads and adapter sequences were removed using Trimmomatic v.0.36 (Bolger et al., 2014). For each independent sample run, the trimmed sequences were assembled using MEGAHIT v.1.2.8 (Li et al., 2015), and the resulting contigs were categorized into either unassigned taxonomy or Ca. Liberibacter based on best BLASTn hits against nucleotide databases. To assess the overall quality of Ca. Liberibacter metagenomes, read coverage statistics of putative Ca. Liberibacter associated contigs were calculated using BBMap (Bushnell, 2014). To improve the binning of Ca. Liberibacter specific contigs, contigs from assembly outputs of each independent sample run were analyzed following the Anvi’o metagenomic workflow (Eren et al., 2015). Putative Ca. Liberibacter contigs were selected from both Anvi’o and best BLAST hit results from NCBI BLASTn with the following constraints, E-value = 1e-4 and sequence similarity >70%. Resulting putative Ca. Liberibacter contigs from independent sample runs were merged for each sample using MEGAHIT v.1.2.8. (Li et al., 2015). The genome assemblies were annotated using Prokka v1.14.5 (Seemann, 2014), and BLASTp searches of output files were conducted against non-redundant NCBI databases.

To determine the evolutionary relationships of putative novel Ca. Liberibacter species, orthologous proteins were obtained from metagenomes generated in this study and those available on NCBI Genbank (below) using Orthofinder (Emms and Kelly, 2015). For phylogenetic analyses, we included a total of four draft Ca. Liberibacter genomes from this study and seven complete Ca. Liberibacter reference genomes from NCBI of Ca. L. africanus (NZ_CP004021.1), Ca. L. americanus (NC_022793.1), Ca. L. asiaticus (NC_012985.3), Ca. L. europaeus (GCA_003045065.1), Ca. L. psyllaurous (NC_014774.1), and L. crescens (NC_019907.1) in phylogenetic analyses. Due to low read coverage of Ca. Liberibacter in one of the Ca. Liberibacter positive psyllid samples (M. gladiata, P1–9) a low number of shared orthologs were identified. In consequence, phylogenetic analyses were conducted only with samples obtained from the psyllids R. capsici (P1–16) and T. walkeri (P1–28), in addition to the herbarium sample, S. umbelliferum (HS9), which shared 102 orthologs with the six fully sequenced Ca. Liberibacter reference genomes (above; Supplementary Table 3).

For the ortholog phylogenetic analysis, 102 orthologous sequences from three draft genomes along with the reference sequences (see above) were aligned using MAFFT v7.407 (Katoh and Standley, 2013). Trimming was performed using trimAL v1.2 (Capella-Gutiérrez et al., 2009) with two parameters, gap threshold “gt” and minimum coverage: “cons” in the trimmed alignment. In this case the gt and cons were considered as 0.9 and 20, respectively. All the alignments were visualized using the graphical user interface, and the PHYLIP file was generated from Mesquite v3.51 (Maddison and Maddison, 2018) and further used for phylogenetic tree construction. RAxML version 8.2.12 (Stamatakis, 2014) was used to obtain phylogenetic inferences. The MPI version was executed to run the program parallelly over many connected machines on the cluster. The “PROTGAMMAAUTO” model along with option (-f a) was used to perform rapid bootstrap analysis and to find the best scoring Maximum Likelihood (ML) tree. 100 searches from the parsimony start tree were performed. The Bipartitions tree generated from RAxML was visualized using FigTree v1.4.4 (Rambaut, 2018). The generated tree was then exported to Adobe Illustrator v 25.1 (Adobe Systems Incorporated, United States) for further editing.

For the 16S rRNA and endosymbiont phylogenetic analyses, similar to above, Ca. Liberibacter and endosymbiont sequences from 16S rRNA gene amplicon sequencing and NCBI GenBank were aligned using MAFFT v7.407 (Katoh and Standley, 2013) and trimmed manually. Mesquite v3.51 (Maddison and Maddison, 2018) was then used to generate the PHYLIP file. To generate the ML tree, the file was submitted to the CIPRES web server to run the RAxML-HPC BlackBox with “GTRGAMMA” model. For the Ca. Liberibacter, Enterobacteriaceae, or Wolbachia phylogeny the outgroup was set as L. crescens, Escherichia coli strain U 5/41, and Wolbachia endosymbiont of Cimex lectularius, respectively, and other settings were kept as default. Visualization and further editing of the bipartitions tree was performed using FigTree v1.4.4 (Rambaut, 2018) and Adobe Illustrator v 25.1 (Adobe Systems Incorporated, United States). The Psylloidea cladogram of 44 psyllid species, that was modified from Burckhardt et al. (2021), was redrawn using the package “ape” (Paradis et al., 2004) in R v.4.0.2 integrated in Rstudio v.1.3.959 (R Core Team, 2021).

If phylogenetic analyses and BLAST identify a potentially new Ca. Liberibacter species DNA--DNA hybridization (DDH) and the average nucleotide identity (ANI) values will be computed pairwise with the new Ca. Liberibacter species and the reference Ca. Liberibacter genomes used above for phylogenetic analyses. DDH values were calculated using the Genome-to-Genome Distance Calculator 3.0 (GGDC²)(Meier-Kolthoff et al., 2013). Briefly, according to Meier-Kolthoff et al. (2013) genome distances were calculated with BLAST+ local alignments and by calculating genomic distance as identities/HSP length (Formula 2). Distance values were then converted to estimated DDH values and their associated confidence intervals using a generalized linear model. The ANI of common genes and the percentage of conserved DNA among the Ca. Liberibacter genomes were calculated using the ANI calculator³ (Rodriguez and Konstantinidis, 2016).

We obtained additional psyllid specimens from species that were identified as Ca. Liberibacter positive from our 16S rRNA amplicon sequencing to determine how prevalent Ca. Liberibacter infection is in these species. We were able to obtain additional psyllid specimens from two of the three Ca. Liberibacter positive species, which includes Trichochermes walkeri (N = 9 samples) and Macrohomotoma gladiata (N = 2 samples) (Supplementary Table 4). We also obtained three additional psyllid species to screen that feed on the same host plant genus (Rhamnus sp.) as one of our Ca. Liberibacter infected species here, T. walkeri, to determine if these psyllid species were infected with Ca. Liberibacter as well; these psyllid species include Trioza rhamni (N = 4 samples), Cacopsylla alaterni (N = 6 samples), and Cacopsylla rhamnicola (N = 7 samples) (Supplementary Table 4). Depending on specimen availability, for each psyllid sample a single adult individual and/or 2–3 nymphs were pooled for each DNA extraction (Supplementary Table 4).

DNA was extracted from whole psyllids using the Quick-DNA^TM Microprep Kit (Zymo Research). DNA quality and yield was assessed by gel electrophoresis and the SpectraMax QuickDrop Micro-Volume Spectrophotometer (Molecular Devices). Based on Ca. Liberibacter phylogenetic analyses here of 16S rRNA amplicon samples (see above), T. walkeri and M. gladiata were infected with Ca. Liberibacter taxa that were closely related to Ca. L. asiaticus and Ca. L. psyllaurous, respectively, in consequence, modified universal Ca. Liberibacter primers, LG774F/LG1436R (Morris et al., 2017), were used for screening additional psyllid specimens. This primer pair is expected to amplify a 684 bp fragment from the 16S rRNA gene region of Ca. L. asiaticus, Ca. L. psyllaurous, and Ca. L. europaeus. PCR was performed in 25 μl reactions consisting of 2 μl of DNA template, 2.5 μl of 10X Thermopol buffer (NEB), 0.5 μl of dNTP mix (2 mM each), 0.5 μl of each 10 μM primer, and 0.125 μl of Taq DNA polymerase (NEB). The Bio-Rad C1000 Touch^TM Thermal Cycler was used with the following protocol: 95°C for 30 s, followed by 30 cycles of 95°C for 30 s, 58°C for 1 min, and 72°C for 1 min, then 72°C for 5 mins. Only T. walkeri samples (N = 9 samples) were positive for Ca. Liberibacter based on PCR screening with universal Ca. Liberibacter primers (above). In turn, amplicons from T. walkeri were purified with QIAquick Gel Extraction Kit (Qiagen) and cloned into the pGEM T-easy vector (Promega) following the manufacturer’s instruction. Sanger sequencing was conducted by the Genomics Core facility at the University of California, Riverside.

Amplicon sequencing of the bacterial 16S rRNA gene region from 45 psyllid samples, which includes 44 psyllid species (Supplementary Table 1), generated a total of 7,516,747 contigs. Following quality control and criteria filtration, 198,560 high-quality contigs (305,128 reads) were retained. The total number of contig counts per psyllid sample ranged from 17 to 4,531 with an average of 798. These contigs were clustered into 91 OTUs, and each psyllid sample possessed an average of 2–3 OTUs, (min = 1, max = 7, and SD = 1; Supplementary Table 5). Non-target OTUs, such as mitochondria, chloroplasts, and Carsonella, only represented 7% of the total reads for high-quality contigs. Seven OTUs were present in two out of the five control samples. Specifically, three OTUs were unique to one control sample (P2–96) and four OTUs were unique to the second control sample (P3–96) (Supplementary Table 6). None of the OTUs in the control samples were present in any of the psyllid samples. The contig counts for control samples were relatively low with an average of 21 contigs per OTU (min = 7, max = 80, and SD = 27) (Supplementary Table 6). In turn, sequences present in control samples may have been airborne contaminants that were present at low abundance during the library preparation process.

From 16S rRNA amplicon sequencing, OTUs that belong to four bacterial phyla (Proteobacteria, Bacteroidetes, Firmicutes, and Fusobacteria) were identified from 45 psyllid microbiomes (Figure 1). The majority of OTUs, among all psyllid samples examined, reside within Proteobacteria (83 OTUs out of 92 OTUs) (Supplementary Table 5), and all psyllid samples possess Proteobacteria OTUs (Supplementary Table 5). Within Proteobacteria, OTUs reside from two classes, Alphaproteobacteria (13 OTUs) and Gammaproteobacteria (70 OTUs), and were found in 24 and 36 psyllid samples, respectively (Supplementary Table 5).

The majority of OTUs (77%) from psyllid samples belong to the class Gammaproteobacteria, and primarily reside within the bacterial Class Enterobacteriaceae (68% of OTUs) (Figure 1). These Enterobacteriaceae OTUs were unique to their psyllid species, except for OTUs 8, 33, 72, 76, and 162, which occur in more than one psyllid sample (Supplementary Table 5 and Supplementary Figure 1). Approximately 55% of Enterobacteriaceae OTUs were not identified to the genus level, however, ∼8, 1, 1, and 2% were identified as Arsenophous, Sodalis, Pantoea, and Enteric bacteria, respectively (Figure 1 and Supplementary Table 5).

Within the class Alphaproteobacteria, five Wolbachia (Anaplasmataceae) OTUs were present in 19 out of 45 psyllid samples (∼42% of all psyllid samples). Wolbachia OTU1 was present in 19 psyllid samples/species (Supplementary Figure 1 and Supplementary Table 5), and the other four Wolbachia OTUs co-occur with Wolbachia OTU1 in three psyllid samples, Cornegenapsylla allophyli, Cacopsylla pruni, and Trioza vitiensis (Supplementary Table 5). Other alphaproteobacterial OTUs include the genus Sphingomonas and Rickettsia (Figure 1). The genus Spingomonas was detected in four psyllid samples (∼9% of all psyllid samples) with two different OTUs; OTU1799 was detected in Pachypsylla celtidismamma, while OTU72 was present in Livia junci, Triozoida limbata, and Arytaina devia ssp. insularis (Supplementary Table 5 and Supplementary Figure 1). The psyllid sample, C. pruni, was infected with two different Rickettsia OTUs, OTU127 and OTU367, and Heterotrioza chenopodii was detected with Rickettsia OTU183 (Supplementary Table 5).

BLAST results from psyllid microbiome 16S rRNA sequences reveal that the vast majority (80%) of sequences had best BLAST hits to insect endosymbionts on NCBI GenBank and displayed an average of >∼96% nucleotide sequence similarity to known insect-associated endosymbionts (Supplementary Table 5). The remaining sequences were primarily found to be plant (2%) or environmental associated microbiomes (10%) (Supplementary Table 5). Sequences that are related to known plant pathogens were ∼5%, including the three Ca. Liberibacter OTUs and a Xylella related OTU (OTU155) (Supplementary Table 5). Only two OTUs (OTUs 618 and 1,355) (2%) in two samples had best BLAST hits to human associated microbes, and therefore these latter OTUs may be contamination from psyllid collection or extraction (Supplementary Table 5).

The obligate endosymbiont of aphids (Buchnera), in addition to other obligate host-restricted endosymbiont species in insects, possess 16S rRNA gene sequences that are significantly higher in AT richness compared to their free-living relatives (Lambert and Moran, 1998; Nováková et al., 2009). Therefore, we analyzed GC content of the 16S rRNA gene region sequenced for OTUs in this study (Supplementary Table 2) and compared them to free-living relatives and known obligate, long-term, and facultative endosymbionts of insects. These free-living relatives and known insect endosymbionts were chosen because they belong to the same bacterial classes that the majority of OTUs from our study were classified as based on Silva, as average GC content can vary with taxonomy (Supplementary Table 2). Specifically, we analyzed: (1) the percentage of GC content from 16S rRNA sequences and (2) the phylogenetic patterns of 16S rRNA sequences that were in the average GC percent range of known primary and long-term insect endosymbionts of insects (Figure 2).

We found that Carsonella, the obligate endosymbiont of psyllids, displays the lowest average GC content (∼35% GC and SD = 0.9) compared to all sequences analyzed in this analysis for the same 402-nt 16S rRNA gene region (Supplementary Table 2 and Figure 2). In contrast, the obligate endosymbiont of aphids, Buchnera, which also belongs to Gammaproteobacteria, shows an average of ∼47% GC (SD = 0.3) content (Figure 2), which was very similar to Profftella, the Gammaproteobacterial long-term symbiont that is harbored within the psyllid genus Diaphorina (Figure 2 and Supplementary Table 2). Other Gammaproteobacterial long-term endosymbionts identified in the psyllids Ctenarytaina eucalypti and Heteropsylla cubana displayed an average of ∼52% GC content (Figure 2 and Supplementary Table 2). Obligate and long-term symbionts that belong to Alphaproteobacteria and Gammaproteobacteria such as Wolbachia and Arsenophonus, respectively, displayed an average of ∼48% (SD = 0 and 0.9, respectively) GC content for this 16S rRNA gene region (Figure 2 and Supplementary Table 2). Interestingly, facultative Wolbachia endosymbionts display very similar average GC content (∼ 47% GC, SD = 0.4) compared to the average of obligate Wolbachia endosymbionts. In contrast, facultative Arsenophonus symbionts displayed an average of ∼52% GC (SD = 1.6) content (Figure 2 and Supplementary Table 2). The free-living relatives of these latter endosymbionts and sequences from this study display an average of ∼53% GC (SD = 1.4), which is in the range of several facultative endosymbionts analyzed in this study (Figure 2 and Supplementary Table 2). Overall, we found that the average percentage of GC content of obligate and long-term endosymbionts ranged from ∼41 to 49% GC (Figure 2 and Supplementary Table 2).

A second signature of long-term endosymbiont evolution is congruent phylogenetic patterns between endosymbiont insect taxa and their insect hosts. For this analysis we only constructed phylogenies of endosymbiont sequences if their 16S rRNA sequence’s percent GC content was between the average range of obligate and long-term endosymbionts (∼41–49%) (above). We found 34 sequences that fell within this range and they either belonged to Enterobacteriaceae or Wolbachia. In turn we constructed two separate phylogenies for these two distinct classes of bacteria (Figure 3 and Supplementary Figure 2). Compared to the psyllid host phylogeny, Wolbachia sequences appear to be recently diverged, largely sharing the same OTU (OTU1) and are horizontally transmitted across all psyllid families examined (Supplementary Figure 2). For the Enterobacteriaceae phylogeny four distinct bacterial lineages were supported with bootstrap values of 87% and higher. However, only two clades (hereafter named clade 1 and 2) appear to have congruent phylogenies with the psyllid host phylogeny at the family level for Triozidae and Psyllidae (Figure 3). More sequence information for both bacterial and psyllid phylogenies are required, however, to provide robust bootstrap support at a finer resolution.

Three Ca. Liberibacter OTUs (OTU526, OTU841, and OTU443) were identified in three different psyllid species, Macrohomotoma gladiata (P1–9), Russelliana capsici (P1–16), and Trichochermes walkeri (P1–28). In M. gladiata, OTU526 has the best BLAST hit to Ca. L. asiaticus with 98% nucleotide sequence similarity (Supplementary Table 5). The OTU841 in R. capsici has the best BLAST hit to Ca. L. americanus with 96.8% nucleotide sequence similarity (Supplementary Table 5). The OTU443 from the psyllid sample T. walkeri, displayed the best BLAST hit to “unclassified Ca. Liberibacter” with 98.3% nucleotide sequence similarity (Supplementary Table 5). To further determine the identity and evolutionary relationships of these latter Ca. Liberibacter taxa psyllid DNA samples that were positive for Ca. Liberibacter were subject to metagenomic sequencing for downstream phylogenetic analyses.

To further analyze Ca. Liberibacter taxa identified from 16S rRNA amplicon sequencing, we performed metagenomic sequencing on the three Ca. Liberibacter positive psyllid samples (M. gladiata, R. capsici, and T. walkeri) and one herbarium specimen infected with Ca. L. psyllaurous haplotype G from a previous study (Mauck et al., 2019). An average of 66,126,618 high quality reads were generated for each independent psyllid/plant sample run (SD = 5,675,882; Supplementary Table 7). Metagenomic assemblies for the four samples yielded an average of ∼0.9 Mb (SD = 0.3 Mb) of Ca. Liberibacter specific sequences per sample with approximately 322–1,637 contigs per sample (Table 1). The Ca. Liberibacter assembly from M. gladiata had very low sequence coverage and a low N50 value (Table 1) and therefore was not included in downstream phylogenetic analyses. The remaining three Ca. Liberibacter assemblies had variable levels of read coverage per contig because DNA samples were of low concentration and obtained from Repli-G whole genome amplification (Supplementary Table 8), however, a total of 102 shared orthologous Ca. Liberibacter genes (Supplementary Table 3) were still obtained and used for subsequent phylogenetic analyses.

A phylogenetic analysis was first conducted with Ca. Liberibacter 16S rRNA gene sequences that were derived from both 16S rRNA amplicon sequencing and metagenomic sequencing. The Ca. Liberibacter species from T. walkeri (P1–28) formed a clade together within other Ca. L. psyllaurous taxa with 92% bootstrap support (Figure 4). The Ca. Liberibacter species from R. capsici (P1–16) forms a distinct lineage apart from other Ca. Liberibacter species (Figure 4). This latter Ca. Liberibacter taxa shares a more recent common ancestor with Ca. L. americanus compared to other Ca. Liberibacter species (Figure 4). Only the 16S rRNA sequence of Ca. Liberibacter from M. gladiata (P1–9) was obtained through 16S rRNA amplicon sequencing and this sequence appears to be more closely related Ca. L. asiaticus (Figure 4 and Supplementary Table 5). The sequence divergence between the 16S PCR-derived sequence and the metagenomic-derived 16S sequence for both P1–16 and P1–28 was present, however, bootstrap support was very high (100%) and nucleotide differences were very small (as the scale on the figure indicates 0.006 nucleotide substitutions per site). This sequence divergence between PCR and metagenomic 16S sequences may be attributed to sequencing error or nucleotide substitutions due to multiple 16S rRNA operons per Liberibacter genome and/or co-occurring strains in pooled psyllid samples.

Since the latter 16S rRNA phylogenetic analysis was constrained in only using a 402 nt sequence region, which was amplified by 16S amplicon sequencing, a more robust phylogenetic analysis was subsequently conducted. The second phylogenetic analysis consists of a 27,922 amino acid alignment of 102 orthologous proteins. From this analysis the Ca. Liberibacter species from R. capsici (P1–16) forms a distinct clade apart from the other Ca. Liberibacter species based on the high amount of amino acid changes per site (Figure 5). The Ca. Liberibacter taxa isolated from T. walkeri (P1–28) and S. umbelliferum significantly clustered within Ca. L. psyllaurous (also known as Ca. L. solanacearum) with 100% bootstrap support (Figure 5).

To further determine if the Ca. Liberibacter species from R. capsici (P1–16) is a distinct Ca. Liberibacter species compared to known Ca. Liberibacter reference genomes we conducted both DDH and ANI analyses. Pairwise DDH estimates for all reference Liberibacter species genomes and Ca. L. capsica ranged from 18.9 to 31.1% (Supplementary Table 9). These values are well below the 79% cutoff commonly applied to members of the same subspecies and the 70% cutoff applied to members of the same species. Pairwise ANI estimates for all reference Ca. Liberibacter species genomes and Ca. L. capsica ranged from 69.72 to 84.49% (Supplementary Table 10). These values are well below the 95% cutoff commonly applied to members of the same species for ANI analyses (Goris et al., 2007). Collectively, these results provide further support that Ca. L. capsica (P1–16) is a distinct Ca. Liberibacter species.

The purpose of the Ca. Liberibacter screening was to (1) determine if the Ca. Liberibacter infection is prevalent in the psyllid species that were identified with Ca. Liberibacter OTUs from 16S rRNA gene amplicon sequencing, and (2) if the Ca. Liberibacter infection occurs in psyllid species that share the same plant host genus as T. walkeri (Rhamnus); samples screened were based on sample availability, see “Materials and Methods” and Supplementary Table 4. Therefore, all DNA extracts from these latter psyllid specimens were screened with the primer set LG774F/LG1463R, which universally amplifies ∼0.7 kb of Ca. Liberibacter 16S rRNA fragments (see “Materials and Methods”). Only T. walkeri samples (N = 9, Supplementary Table 4) produced a positive PCR band using this primer set, and therefore were cloned and Sanger sequenced one PCR clone from each sample for BLAST analysis.

Based on NCBI BLAST results, the best BLAST hit (E-value = 0) from the nine T. walkeri samples was the taxon Ca. L. psyllaurous with >97% sequence identity (Supplementary Table 11). Although all the T. walkeri samples matched Ca. L. psyllaurous with high sequence identity, there was small nucleotide differences between sample sequences, which may be attributed to the multiple 16S rRNA operon copies within a single Ca. Liberibacter genome and/or strain differences within or among psyllid individuals (Supplementary Table 11). These screening results, in terms of identity (Ca. L. psyllaurous) and sequence variation, are consistent with 16S rRNA amplicon sequencing and metagenomic sequencing from the T. walkeri (P1–28) sample, and in turn, this may indicate that Ca. Liberibacter psyllaurous is highly prevalent in some T. walkeri populations.