We obtained all publicly available genomes of the genus Liberibacter from the National Center for Biotechnology Information (NCBI) and Pathosystems Resource Integration Center (PATRIC) database on March 24th, 2022. We obtained 36 genomes of ‘Ca. Liberibacter asiaticus,’ 2 genomes of ‘Ca. Liberibacter americanus,’ 2 genomes of ‘Ca. Liberibacter europaeus,’ 2 genomes of ‘Ca. Liberibacter africanus,’ and 8 genomes of ‘Ca. Liberibacter solanacearum’ (Supplementary Table S1). Additionally, we included 2 genomes of Liberibacter crescens (Table 1). Genome completeness was assessed using CheckM v 1.2.0 (Parks et al., 2015), based on a set of conserved single copy genes (Supplementary Table S2). Three genomes of ‘Ca. Liberibacter asiaticus’ were excluded (genomes Tabriz3, SGCA_1, and SGpsy) because their genome sizes were considerably lower than the accepted 1.2 Mb genome size of Las (0.769 Mb and 0.23 Mb) or their genome completeness was less than 80% based on CheckM, likely reflecting incomplete genomes.

All genomes retrieved from public databases were annotated using Prokka v 1.14.6 in order to obtain uniform annotations (Seemann, 2014). Following annotation, a genus-wide pan-genome was constructed using Roary v. 3.13.0 with Rhizobia grahamii as the outgroup (Table 2). Due to the high divergence amongst Liberibacter species, we decreased the percentage sequence identity to 50% to encompass all possible orthologs following previously established methods (Chase et al., 2018; Rodriguez and Martiny, 2020). Using the pan genome reference fasta file output by Roary, we extracted a representative nucleotide sequence of each core and accessory gene using cdbfasta.¹ Theses gene sequences were translated using the transeq command from EMBL-EBI (Madeira et al., 2019). Functional categorization of the core genes and accessory genes was performed using the program eggNOG-mapper v. 2 (Huerta-Cepas et al., 2017, 2019). KEGG annotation was performed using the GhostKoala online server (Kanehisa et al., 2016). Additionally, we measured the pairwise average nucleotide identity (ANI) between all genomes used in this study using OrthoANI v. 1.2 (Lee et al., 2016).

The core gene alignment was used to build a phylogenetic tree representing the entire genus. The core gene alignment produced by Roary was polished using the gBlocks program (Castresana, 2000). The polished core gene alignment was inputted to RAxML v. 8.2.12 to build a phylogenetic tree (Stamatakis, 2014). In order to further evaluate and curate the isolates, we created a distance matrix from the RAxML phylogenetic tree in newick format using the Tree and reticulogram REConstruction (T-REX) server (Boc et al., 2012). Visualization of the phylogenetic tree built from core gene sequences was performed using the Interactive Tree of Life (iTOL; Letunic and Bork, 2021). Additionally we explored the phylogenetic relatedness of the accessions based on accessory gene content. We constructed a distance matrix based on the presence-absence patterns of accessory genes and then used that matrix to build a neighbor-joining tree. To build the tree, we calculated the Euclidean distances from the presence-absence matrix using the dist function in R (R Core Team, 2021). The distance matrix based on accessory gene content was then used as input to build the tree using the ape package in R (Paradis et al., 2004). Using the Analysis of Similarities (ANOSIM) in the R vegan package, we tested the association between the accessory gene content and species identity (Oksanen et al., 2020).

To investigate the rearrangement history of Liberibacter and to generate a multiple whole genome alignment, we used the progressiveMauve system available through Mauve (Darling et al., 2004). Before performing multiple genome alignment with Mauve, we used the Mauve Contig Mover (MCM) on assemblies with >1 contig using the Las isolate Psy62 or Lso isolate CLsoZC1 as the reference genome, depending on the analysis (Darling et al., 2007). ProgressiveMauve was run with default settings and the Locally Collinear Blocks (LCBs) that were found across all genomes were extracted and concatenated into a multiple genome alignment. We performed progressiveMauve analysis with 10 different sets of genomes: a set of only single contig assemblies (n = 16), a set of genomes comprising two representative sequences per Liberibacter species (n = 12), and 8 different species-specific analyses (Table 3). For the interspecific analyses, Psy62 was used as the reference genome for MCM. However, for the Lso species-specific analyses, we used CLsoZC1 as the reference genome for MCM. For the species-specific analysis of Las genomes, we used subsets of Las genomes with n = 10 or n = 8 because Mauve was too computationally intensive to run the analyses on all 36 Las genomes. We therefore investigated three sampling subsets. The first subset was composed of all 10 single-contig Las assemblies (chromosome level); the second was composed of a random subsample of 8 Las genomes from the entire set of 36 genomes; and the third set included a random subsample of 8 Las genomes from the set of 26 Las genomes composed of >1 contig. We selected 8 Las genomes for the second and third analysis sets in order to match the number of available Lso genomes (n = 8).

We utilized the codeml set of models from PAML v.4.9 to calculate the ratio of nonsynonymous to synonymous substitution rates, known as dN/dS or ω, on the sets of single-copy core and single-copy accessory genes (Yang, 1997, 2007). For all models tested, we only considered genes with a minimum of four sequences (i.e., four different accessions carry the gene), as this is the minimum number recommended for codeml analysis.² The codeml analysis in PAML requires the input of an unrooted, maximum-likelihood phylogenetic tree for each sequence analyzed, which we constructed using RAxML v. 8.2.12 using the gene sequence alignment as input for tree construction (Stamatakis, 2014).

We performed the codeml analysis under various models and compared the model results using likelihood ratio tests (Yang, 2007). Using the various models, we could test the null hypothesis that ω = 1.0 against the alternative of positive selection (ω > 1) using two different methods. First, in the global test approach, we calculated a single estimate of ω for each gene across its entire phylogeny. To test for the validity of the estimated value of ω, we performed a comparison of two models using the likelihood ratio test: one that estimated a single ω for each gene from the data (Model = 0 with Fix_omega = 0 in the codeml control file) and another model that calculated the likelihood if ω = 1.0 (Model = 0 with Fix_omega = 1 and Omega = 1 in the codeml control file). Using this set of models, we could identify genes with evidence of positive selection when the initial ω estimate was >1.0 and when the likelihood of the two models differed significantly based on p < 0.01 after FDR p-value correction. The second analysis consisted of the Sites models approach which tests for genes with variable selection pressure across sites, or codons, of the sequence. For each gene, we first compared Sites models M0 and M3 using the likelihood ratios to test for the presence of heterogeneity in evolutionary rates across codons. Any gene with a significant result for the M0–M3 comparison was then compared using Sites models M1a and M2a to identify specific codons with evidence of positive selection (ω > 1.0). For all summary statistics of ω, we excluded estimates of ω that were greater than 20 as potentially unreliable due either to low dS or poorly resolved optimization.

All statistical analyses were performed in R v 4.1.2 (R Core Team, 2021). In order to compare the abundance of genes involved in KEGG functional pathways and ANI across species, we first assessed the normality of the data using the Shapiro–Wilk test and the variance of the data using Levene’s test available through R. For all comparisons of KEGG functional pathways, we performed ANOVA analyses using the R function aov and identified which groups were significantly different from each other by computing a post hoc Tukey HSD test with the R function TukeyHSD. For ANI comparisons, the Shapiro–Wilk test for normality was significant, so we proceeded with a Kruskal-Wallis analysis using the kruskal.test function in R. The Kruskal-Wallis test revealed significant differences between the groups, so we followed with the Pairwise Wilcoxon Rank Sum Test using the R function pairwise.wilcox.test with Benjamini and Hochberg p-value correction for multiple testing.

To investigate the genus wide diversity and evolution of Liberibacter species, we obtained 52 Liberibacter whole genome sequences from public databases, representing one culturable and five nonculturable Liberibacter species. We also retrieved the genome sequence from Rhizobium grahamii as an outgroup sequence. Our complete sample of 53 whole genome sequences encompassed six Liberibacter species, with one species represented by 36 accessions (Table 1; Supplementary Table S1). We constructed the genus-wide pan-genome and categorized each gene as either a core or accessory gene in order to construct a core gene alignment to build a phylogeny. We defined a core gene as a gene that was present in a majority (95%) of the samples being tested (Tettelin et al., 2008). Across the entire genus, we identified 439 core and 8,487 accessory genes. Among accessory genes, 5,807 were singletons, meaning they were only identified in a single accession. We excluded singletons from further analysis and also excluded genes with an average number of sequences per isolate greater than one (due to splitting errors). Our final dataset consisted of 436 single-copy core genes and 2,639 single-copy accessory genes for downstream analysis (Table 2).

With the genus wide pan-genome information, we constructed a core gene alignment to use as input to build a maximum likelihood phylogeny (Figure 1). As expected, the sequences clustered by species, with L. crescens basal among Liberibacter species (Wang N. et al., 2017). Following the split of L. crescens from the rest of the Liberibacter species, Lam and Leu sequences formed a clade that diverged from Lso, Laf, and Las. Overall, the phylogeny’s topology was highly supported with 100% bootstrap support across all nodes leading to species groups. Across all nodes in the tree, including within species clades, the median bootstrap support was 61.0% and the mean was 58.7% (Supplementary Figure S1). The lowest bootstrap supports were primarily found at nodes near the tips of the tree and only within species groups, such as in the Las clade, suggesting lower evolutionary divergence among these samples and their core gene sequences. The species tree based on core genes recapitulates the previous analysis of Thapa et al. (2020).

We performed functional analyses on the core and accessory gene sets using the COG annotation scheme (Figure 2). After excluding the category of S - ‘function unknown’ (which was the second largest category), the largest COG category among core genes was J - ‘translation, ribosomal structure, and biogenesis’ (96 genes), followed by F - ‘nucleotide transport and metabolism’ (32 genes) and H - ‘coenzyme transport and metabolism’ (31 genes). The accessory genes differed, because the largest categories (excluding ‘function unknown’ genes, which was the largest category) were L - ‘replication, recombination, and repair’ (233 genes), M - ‘cell wall/membrane/envelope biogenesis’ (166 genes), and E - ‘amino acid transport and metabolism’ (135 genes).

To investigate Liberibacter intra- and interspecies divergence, we calculated the genome-wide average nucleotide identity (ANI) using OrthoANI between all possible pairs of accessions included in this study, resulting in 1,378 pairwise ANI calculations (Figure 3). ANI is a commonly used metric to demarcate species and genus boundaries in bacteria (Jain et al., 2018). OrthoANI calculates high resolution ANI across the entire genome by fragmenting the genome sequences into 1,020 bp long fragments that are used for analysis (Lee et al., 2016). Across all accessions (including the outgroup species), the mean ANI value was 87.13% with a minimum value of 62.85% and a maximum of 99.99%. As expected, comparisons performed within species had a significantly higher ANI compared to the ANI calculated from across species comparisons (p < 2.2 × 10⁻¹⁶, Wilcox Test). Specifically, average ANI for between species comparisons was 76.2% (excluding calculations against the outgroup) while average ANI within species was 99.8%. As expected, ANI both within and between Liberibacter species were significantly different to the ANI calculated against the outgroup species (mean ANI against outgroup = 66%, p < 2.2 × 10⁻¹⁶, Kruskal–Wallis test; Figure 3).

Next, we investigated ANI within the Liberibacter species. We hypothesized that Lso isolates may have greater nucleotide diversity, and therefore lower mean ANI, compared to other Liberibacter species, like Las, which in turn may underlie the broader host range of Lso. We focused on the within-Lso and within-Las ANI comparisons as these two species had the greatest number of available genomes with 8 and 36 genomes, respectively. The mean within-Las ANI was 99.89% while the mean within-Lso ANI was 98.17% (Supplementary Figure S2). ANI was significantly higher for within-Las comparisons than within-Lso comparisons, indicating that Lso has greater nucleotide diversity compared to Las genomes (p < 0.0001, Wilcox Test). All other within-species calculations for Laf, Lam, Leu, and L. crescens resulted in ANI values of >99%, but only two genomes were available for each of these species. Overall, however, these analyses suggest that Lso has greater within-species nucleotide diversity compared to other Liberibacter species.

We also investigated genomic differences among species. There were significant differences in total genome size between species (p = 0.001, Kruskal–Wallis test); for example Las isolates tended to have significantly smaller genomes compared to L. crescens and Leu isolates (p < 0.05). We also identified species’ differences by counting the core and accessory genes for each species individually (Table 2; Supplementary Figure S3). Intraspecies pan-genome analysis returned an average of 985 core genes per species with a minimum of 786 core genes in Las and a maximum of 1,291 core genes in L. crescens. A total of 536 genes were shared among the pathogenic species (Supplementary Figure S3A). The accessory gene count was more variable among the species, with an average of 291 accessory genes per species but a minimum and maximum of 44 (in Leu) and 768 (in Las). We recognize that the number of accessory genes is positively correlated with the number of genomes analyzed up to a saturation point (Bosi et al., 2015); hence, we are undoubtedly underestimating the number of accessory genes in species with only two genomes available. Among species with only two genomes available, Laf had 229 accessory genes, Lam had 65, and Leu had 44 resulting in an average of 112 accessory genes in these species.

Using the COG functional annotation scheme, we explored the species-specific core and accessory gene functions together (Figure 4). As expected (Merfa et al., 2019), L. crescens tended to have more genes with annotations in particular COG functional categories such as E - ‘amino acid transport and metabolism,’ H - ‘coenzyme transport and metabolism,’ O - ‘post-translational modification, protein turnover, and chaperones,’ P - ‘inorganic ion transport and metabolism,’ Q - ‘secondary metabolites biosynthesis, transport, and catabolism,’ and T - ‘signal transduction mechanisms,’ which may contribute towards L. crescens ability to be cultured or may be due to an inherently larger chromosome. Excluding L. crescens, the distribution amongst the COG categories was fairly consistent among the unculturable ‘Ca. Liberibacter’ species. COG categories J, L, and M revealed interesting patterns amongst Las, Laf, and Lso, as all three species tended to have higher counts in these categories compared to Lam and Leu. These categories were ‘translation, ribosomal structure, and biogenesis,’ ‘replication, recombination, and repair,’ and ‘cell wall/membrane/envelope biogenesis,’ respectively.

We compared and contrasted the core and accessory gene content between Las and Lso isolates, as these two species had >2 genomes available, allowing for greater confidence in the core and accessory gene representation in these species. Las and Lso had a total of 786 and 822 core and 768 and 547 accessory genes, respectively. In regard to core genes, the majority (75%) of core genes were the same between the two species as they shared 689 core genes in common. Each species had a set of unique core genes, however, with Las having 97 and Lso having 133 unique core genes. The accessory genes proved to differentiate the two species further, as Las and Lso only shared 354 accessory genes which represented just 37% of the total genes. The Lso species harbored a set of 414 unique accessory genes that were not a part of the Las accessory gene set, while Las only had a set of 193 unique accessory genes. Overall, this comparison suggests that Lso not only harbors more genetic diversity within core genes, as measured by ANI, but also more variation in accessory gene content.

Since we found differences in gene content among the Liberibacter species, we compared and contrasted the abundance of genes involved in different processes essential to the biology of the cell: nucleotide, amino acid, energy, and carbohydrate metabolism using KEGG functional category definitions (Kanehisa et al., 2016). Liberibacter bacteria are largely restricted to the phloem of plants, which is rich in sugars and amino acids, so we hypothesized that the different species of Liberibacter would exhibit differences in the number of genes involved in the metabolism of these compounds due to differences in pathogenic ability. For each genome analyzed, we counted the number of genes in each category and investigated the species-specific differences (Figure 5). We found a significant difference in the number of amino acid, energy, and nucleotide metabolism genes between the different species (p < 0.001, ANOVA). However, the number of genes involved in carbohydrate metabolism was not significantly different between the species (p = 0.13), likely reflecting the importance of carbohydrate metabolism and constraint against losing this function. Isolates from the species L. crescens had significantly higher abundance of amino acid and energy biosynthesis related genes compared to all other ‘Ca. Liberibacter’ species. In contrast, the Las isolates had a significantly higher number of nucleotide biosynthesis genes compared to Lam, Leu, and Lso isolates.

We also investigated large scale genome rearrangements and changes in genomic synteny across the Liberibacter genus, as this provides another way to investigate evolutionary relationships within the genus. Previous work has demontrated extensive genome rearrangement in Liberibacter genomes, but only in a handful of genomes (Thapa et al., 2020; Tan et al., 2021). We have used progressiveMauve to preform extended whole genome alignment on two sets of Liberibacter genomes: a set of 12 genomes representing the two highest quality genomes from each of the six species and a second set representing the 16 available chromosome-level genome assemblies representing five species of Liberibacter: Laf, Lam, Las, Lso, and L. crescens.

The alignment of 16 assemblies revealed a total of 136 synteny blocks, or locally collinear blocks (LCBs), across all 16 genomes. An LCB refers to a conserved stretch of the genome with no evidence of rearrangement history. The 136 LCBs spanned a sequence length of 1,167,307 bp and had an average pairwise identity of 60.5% across isolates (Table 3). We compared the phylogenetic relationships in the tree built from the core gene alignment to a phylogeny built from pairwise distances inferred from the Mauve whole genome alignment (Supplementary Figure S4). Considering just the tips of the tree, 12/16 of the isolates maintained the same position in the tree; all four isolates with different positions were Las isolates. The congruence between the two trees indicates that the core gene alignment and whole genome alignment recapitulate similar phylogenetic histories, but with evidence for phylogenetically informative rearrangements in Las.

To investigate the genus wide rearrangement history, we also performed progressiveMauve alignments with the set of 12 most-contiguous genomes representing all six species and within individual species (Table 3). We performed the alignment several times: once all together with all 12 genomes and again once for each species individually to compare the amount of LCBs, which reflects the amount of rearrangement (Darling et al., 2007). The alignment with all 12 of the representative genomes identified156 LCBs spanning a total of 1,148,204 bp, with a pairwise sequence identity of 44.6%. As previously described, we compared the phylogeny built from this whole genome alignment to the phylogeny built from the core gene sequences (Supplementary Figure S5), and we found no differences between the two trees.

Additionally, we performed progressiveMauve alignments for each species individually in order to investigate any species-specific rearrangements and to compare LCB counts within and between species. In contrast to genus-wide analysis, species-specific mauve analysis retrieved very few LCBs – i.e., an average of 11.25 within species. Lso isolates had the most synteny blocks with 46 LCBs identified, followed by the Lam group with 22 LCBs (Table 3). Because the Las group had the most genomes available (n = 36) and Mauve is computationally intensive, we performed the progressiveMauve analysis with three different sets of Las genomes: single contig assemblies (n = 10), a random set of eight genomes from the entire set of 36 Las genomes, and a random set of eight genomes composed of >1 contig from a set of 26 Las assemblies. The analysis revealed 1, 3, or 11 LCBs depending on the set of genomes analyzed (Table 3). Overall, these analyses indicate that the Liberibacter genus has an extensive rearrangement history across species, with 156 synteny blocks shared across the entire genus. As expected, there were fewer LCBs within species, reflecting fewer obvious genome rearrangements within species, but the Lso group was notable for its relative high number of LCBs.

Genome rearrangement, reduction, and recombination is prevalent in the Liberibacter genus, and it has been hypothesized to be important towards pathogenicity and virulence evolution (Tan et al., 2021). To better understand these issues, we first investigated the evolutionary patterns of the accessory genes by comparing the core gene phylogeny against a phylogeny based on the accessory gene composition (Figure 6). The accessory gene phylogeny recapitulates a similar clustering pattern as the core gene phylogeny, which has accessions clearly grouped by species. However, the correlation between the two distance matrices used to build each tree was non-significant based on the Mantel test (p = 0.086, R = 0.158, Mantel Test). This discordance was prominent in the overall species clustering as the placement of the Laf accessions was positioned at a deeper portion of the accessory gene tree compared to the core gene tree. The core gene tree positioned the Laf accessions more closely to Las in a single clade that diverges from the Lso isolates, while the accessory gene tree placed the Laf accessions farther from the Las group and closer to the Leu and Lam clade. The overall pattern of accessory gene content suggests that isolates of the same species have more similar accessory gene content and also that gene exchange between species are not occurring at a rate to alter phylogenetic signal, with the except of the placement of Laf accessions. This statement is supported by ANOSIM, which detected a significant association between accessory gene content and species (ANOSIM R = 0.4525, p < 2 × 10⁻⁴).

We next investigated the correlation between the core gene tree and accessory gene tree within Las and Lso, because those groups had >2 accessions to compare. In the Lso clade, the position of 6 of the 8 Lso accessions remained the same across the two trees, suggesting that the core gene nucleotide evolution and accessory gene composition recapitulate similar evolutionary histories (Figure 6). Isolates CLsoZC1 and ISR100 experienced positional changes between the two trees. In the core gene tree, CLsoZC1 diverged from the 7 other Lso accessions; however, in the accessory gene tree, CLsoZC1 formed a group with 4 other Lso isolates. The Lso isolate ISR100 also changed position, as it was originally grouped with LsoNZ1, HenneA, R1, and RSTM in the phylogenetic tree built from the nucleotide core gene alignment, however, in the accessory gene tree, ISR100 formed a group with FIN111 and FIN114. Unlike Lso, the Las group showed extensive changes between the two trees. Most notably, we observed the presence of three clusters within the Las clade based on the core gene alignment, but this clustering pattern was not seen in the tree constructed from accessory gene composition. Altogether, the results suggest that accessory gene content varied by species, leading to similar phylogenetic inference across species based on core genes or gene presence/absence. However, there is enough gene content variation within Liberibacter species to create slight differences in the phylogenetic patterns within species.

To investigate the strength and pattern of selection on individual genes, we estimated the ratio of nonsynonymous mutations to synonymous mutations, referred to as ω or dN/dS. By calculating this ratio, we can identify genes with evidence of positive selection (ω > 1.0) which suggests that these genes may be relevant towards pathogenicity and virulence (Mitchell et al., 2012; Aleru and Barber, 2020). Using the codeml models available through the program PAML, we applied a series of nucleotide substitution models on single-copy core and accessory genes identified in the pan-genome analysis. We tested for positive selection using two methods: either globally across the entire phylogeny (referred to as the global test) or across codon sites (referred to as sites models). In total, we performed these analyses on 1,698 qualifying genes (436 core genes and 1,262 accessory genes; see Methods).

Using the global test, we estimated the single, average ω value for each gene across the phylogeny of the genus. This method assumes that the value of ω is constant across all of the branches in the gene tree and across the individual codons in the nucleotide alignment. For the 436 core genes, estimates of ω ranged from 0.00649 to 0.26883 with an average of 0.09442 (Figure 7A). All of the core genes had estimates of ω less than 1.0, and the vast majority (385/436) had ω estimates significantly less than 1.0 (p < 0.01, FDR correction), indicating a pervasive signal of selective constraint on amino acid replacing mutations. For the 1,262 accessory genes, ω ranged from 0.001 to 2.85930 with a mean of 0.22465 (Figure 7B). Among these, we identified 36 accessory genes with estimates of ω > 1, suggesting the possibility of positive selection, but only two had ω estimates significantly higher than 1.0 (p < 0.01, FDR correction). We compared the values of average ω estimated for the core and accessory genes; the two groups were significantly different, with accessory genes having significantly higher ω on average (p < 2.2 × 10⁻¹⁶, T-test).

We also applied the sites models in codeml, which test for variation in ω among codons in the gene sequences and also identify whether any sites have a history of positive selection. This represents a less conservative method than the global test to search for positive selection. The sites models are a group of nested models, so we first compared sites model M0 against sites model M3. Sites model M0 represents the null hypothesis that there is a single ω value for all sites in the gene sequence, while sites model M3 permits ω to vary among sites. The genes that were significant for the likelihood ratio test between sites model M0 and M3 were then analyzed under the sites models M2a and M1a which tests for a history of positive selection across codon sites. All but one of the core genes rejected the null hypothesis that a single ω value exists across all sites in the sequence, suggesting variable ω across codons (p < 0.01, FDR correction). (The gene ftcR was the one non-significant exception.) However, we did not find any core genes that had evidence of a history of positive selection in any of their codons based on M2a-M1a comparisons. The accessory genes displayed contrasting dynamics, because only 64% of the accessory genes (806/1,262) rejected the null that a single ω value exists across all sites in the M0–M3 comparison. We investigated the 806 genes further using the M2a-M1a model and found that 52 accessory genes had significant evidence of a history of positive selection among sites in the sequence (p < 0.01, FDR correction). Among these 52 accessory genes, ten also had evidence of positive selection via the global test, as indicated by estimates of ω > 1 (Table 4).