Pseudogenization of the Chaperonin System in ‘Candidatus Phytoplasma pruni’ Reveals Insights into the Role of GroEL/Cpn60 in Phytopathogenic Mollicutes

1. Introduction

Mollicutes, so named because of their lack of a cell wall (Latin mollis, soft, + cutis, skin), are a major group within the “strong-walled” Gram-positive phylum Firmicutes. This class of bacteria is believed to have evolved around 600 million years ago [1] through a reductive process whereby the microorganisms reduced the metabolic capacities and structural features encoded on their genomes as they adopted an obligately intracellular lifestyle. Extant Mollicutes are universally dependent on eukaryotic host cells to complement their limited metabolic functions, although a few, notably Mycoplasmas, can be cultured in complex artificial medium and are a common contaminant of eukaryotic cell cultures [2]. Mollicutes are commonly associated with pathogenesis in plants (Spiroplasma, Phytoplasma) [3,4] and animals (Mycoplasma) [5]. Due to their reduced genomes, Mollicutes are among the bacteria with the smallest genome size - M. genitalium possesses one of the smallest known bacterial genomes at 580 kbp [6]. Their lack of a cell wall results in a generally pleomorphic shape, with the exception of the helical Spiroplasmas [2]. Other unique features of Mollicute genomes include a low G/C content, and the use of the codon UGA for tryptophan rather than a stop codon as in other bacteria [7]. The single exception to the latter characteristic is the Family Acholeplasmataceae, which consists of the two genera, Acholeplasma and Phytoplasma. Since these taxa use UGA as a stop codon, which is characteristic of non-Mollicute bacteria, they are presumed to be more closely related to the ancestral taxon [1,8].

The protein chaperonin system (GroE) consists of the two proteins GroEL (synonym Cpn60) and GroES (synonym Cpn10) [9]. This system is canonically involved in the prevention of aggregate formation and proper folding of cellular proteins as a part of the basic protein biosynthetic machinery, and for this reason these two genes are anticipated to be present in all prokaryotic and eukaryotic cells [10]. While only a subset of cellular proteins require the GroE system for their tertiary structure formation, these include proteins that are essential for cell viability, so that a loss of this system is thought to result in a cell that is nonviable [11]. Mollicutes, however, are again unique in that they are the only cells known generally to lack a GroE system, although it is present in a subset of these microorganisms [12]. Because the GroE system is an ancestral property (i.e., it is present in all non-Mollicute bacteria), it is thought that Mollicutes lost these genes during genome reduction [11,12]. The retention, or re-acquisition, of GroE is polyphyletic, with two groups of Mollicutes featuring these genes: a smaller clade consisting of Mycoplasma gallisepticum, M. genitalium, and M. pneumoniae, and a much larger group including nearly all of the Acholeplasmas and Phytoplasmas [11]. In addition, there are at least three clear instances of the acquisition of groEL by horizontal gene transfer (HGT) in Mollicutes: the gene of Mycoplasma penetrans is more closely related to that of Helicobacter pylori than other Mycoplasmas [12]; and two species of Spiroplasma, S. turonicium and S. kunkelii, appear to have acquired groEL genes by HGT [11].

These observations have led to speculation regarding the role of these proteins within these intracellular microorganisms, particularly since a very wide array of non-canonical “moonlighting” roles has been ascribed to GroEL in prokaryotic and eukaryotic cells [13]. Notably, Clark and Tillier [12] proposed that, among Mollicutes that encode GroEL in their genomes, the protein may play a role in pathogenesis by acting as an adhesin/invasin, suggesting that the re-acquisition or retention of the GroE chaperonin system may be associated with virulence. This notion is supported by a wide variety of observations that non-Mollicute bacteria can localize GroEL to the cell surface, or even secrete the protein into the extracellular space, clearly indicating a non-canonical, virulence-associated role for the protein in certain bacteria. For example, H. pylori localizes its GroEL to the cell surface, where it interacts with Toll-like receptors to induce interleukin-8 production in epithelial cells, indicating a role in gastric inflammation [14]. In addition, Legionella pneumophila, a non-Mollicute obligate intracellular pathogen, also displays GroEL at the cell surface, and intracellular forms of the pathogen are enriched for GroEL at the cell surface; moreover, L. pneumophila GroEL-coated latex beads are internalized by HeLa cells, where they are maintained as endosomes, mimicking an infective status [15,16]. Brucella abortus, another intracellular pathogen, uses its GroEL to bind to host prion protein and facilitate cell invasion [17]. Other non-Mollicute intracellular pathogens also use the GroEL encoded by one of the paralogous genes to gain access to their host cells and trigger pathogenesis, such as Chlamydia pneumoniae [18] and Mycobacterium tuberculosis [19]. These observations, along with many others in different pathosystems, provide strong reasons to consider carefully the role of GroEL/Cpn60 and the GroE system in pathogenic Mollicutes that have retained or re-acquired these genes.

Phytoplasmas, plant pathogenic microorganisms that are transferred between host plants through an insect vector [20], are among the two groups of Mollicutes that are typically considered to have retained or re-acquired a GroE system [11]. However, an exception to this was noted when the genomes of four members of the ribosomal group 16SrIII, classified as ‘Candidatus Phytoplasma pruni’ were sequenced in 2012 [21]. The genomes of phytoplasmas causing diseases known as Italian clover phyllody, Poinsettia branch-inducing phytoplasma disease, milkweed yellows, and Vaccinium witches’ broom were notably lacking a GroE chaperonin system. The lack of groEL genes in a clearly pathogenic branch of the phytoplasmas calls into question the role that GroEL may play in the pathogenesis caused by these microorganisms. However, all of these genomes are highly fragmented, with at least 150 contigs present in each. Phytoplasma genomes have proven to be difficult to sequence prior to the introduction of long-read sequencing technologies, since the inability to culture them means that the genomes are always metagenome assembled genomes. Despite technological developments for sample preparation including antibody-based enrichment [22], many phytoplasma genomes deposited in public repositories remain fragmented and incomplete. This hinders the search for groEL genes in phytoplasma genomes and can mask the possibility of groEL gene acquisition by HGT. Thus the reported lack of a GroE system in this singular group of phytoplasmas remains an unresolved question.

We address this issue by examining all reported genomes of ‘Candidatus Phytoplasma pruni’ for the presence of groEL/groES genes or pseudogenes, and provide two new high-quality genome assemblies for 16SrIII phytoplasmas in subgroups A and F. We report the pseudogenization of groEL in some of the reported 16SrIII genomes and a complete lack of groES in all of these genomes, and have confirmed through analysis of all reported high-quality phytoplasma genomes that the 16SrIII group constitutes the only phytoplasmas that truly lack a chaperonin system. This observation casts doubt on the previously speculated moonlighting role of GroEL in the attachment and invasion of host cells and suggests that the protein may play a canonical role in the folding of client proteins in these microorganisms.

2. Materials and Methods

2.1. Strain Sources, DNA Isolation, and Genome Sequencing

Syringa (lilac) infected with ‘Ca. P. pruni’ (16SrIII) strain 2A1 was maintained at the Centre for Plant Health, Canadian Food Inspection Agency, North Saanich, British Columbia, Canada. Symptomatic milkweed (Asclepias sp.) was observed at the Canadensis Botanical Garden, Central Experimental Farm of Agriculture and Agri-Food Canada, Ottawa, Ontario, Canada (https://maps.app.goo.gl/urU9SzTiAYX2Zba59). Phytoplasma infection of the milkweed tissue was confirmed using a qPCR assay targeting 16S rRNA genes as described [23]. The phytoplasma was typed using nested PCR amplification of 16S genes (F2nR2) [24,25], and cloned amplicons were sequenced by Sanger (Eurofins). The RFLP type was determined for six clones using the iPhyclassifier [26]. For both plant types, leaf tissue (midrib and petiole) was cut from infected leaves and homogenized in liquid nitrogen using a mortar and pestle that was pre-treated with DNAaway (Thermofisher). DNA was extracted from the plant material using a Wizard High Molecular Weight DNA extraction kit (Promega). DNA yield was determined using a Qubit fluorometer (ThermoFisher), and molecular weight was assessed using agarose gel electrophoresis. This sample was named MYp-CanS4.

Bacterial DNA was enriched from the leaf DNA extract using a NEBNext Microbiome Enrichment Kit (New England Biolabs) as previously described [27]. Both short-read Illumina and long-read Nanopore sequencing technologies were used for genome sequencing. For Illumina, 0.5 μg of enriched DNA was used as the input material for each DNA library preparation. The sequencing libraries were generated using a NEBnext UltraII FS DNA lib prep kit E7805 >100 ng, according to the manufacturer’s instructions, with the following specifications: 10 min fragmentation, size selection using 30 µl / 10 µl beads for 200-475 bp, 5 cycles of amplification, and use of dual index barcodes using NEBnext Multiplex Oligos for Illumina (NEB, USA). The DNA libraries were sequenced on an Illumina NovaSeq 6000 platform (Illumina, San Diego, USA), and 150 bp paired-end reads were generated. For Oxford Nanopore Technology (ONT) sequencing, 0.12 µg for bacterial-enriched lilac and 0.525 µg for bacterial-enriched milkweed DNA were used as the input material for the library preparations. The sequencing library was prepared using an ONT Ligation Sequencing Kit (LSK-SQK114), followed by sequencing on a MinION Mk1C device using an R10.4.1 flowcell (ONT, Oxford, UK).

Illumina reads were quality trimmed using trimmomatic [28], by removing bases from the beginning of reads with a quality score below 3, from the end of reads with a quality score below 20, and cutting any reads with a quality score below 15 over a sliding window of 4 bases. Reads shorter than 36 bases after trimming were removed. The Illumina reads from the lilac-2A1 sample were then mapped (without merging, to minimize reads loss) using a reference database prepared from ‘Ca. P. pruni’ strain PR2021 (GCA_029746895) with bowtie2 [29]. The MYp-CanS4 Illumina reads were mapped using bwa [30] with a reference database prepared from the previously reported milkweed yellows genome (GCA_000309485). Nanopore reads were not mapped but were filtered to exclude reads under 1000 bp using filtlong ( https://github.com/rrwick/Filtlong). The top 90% of reads were retained, based on a composite quality score reflecting the average base quality and read length using filtlong. Illumina and nanopore reads were co-assembled using unicycler [31]. Assemblies were examined for the possible presence of plasmids using plasmer [32], and assembly completeness was assessed using Busco [33] with the default prokaryotic reference database (bacteria_odb12). Assembled genomes were annotated locally using bakta [34], as well as by NCBI using the Prokaryotic Genome Annotation Pipeline (PGAP) [35].

A separate hybrid assembly was prepared for the MYp-CanS4 sample by first removing host reads (including chromosomes, mitochondria, and chloroplast) using bmtagger [36]. Chromosomal sequences for Asclepias syriaca were downloaded from milkweedbase.org, while chloroplast (NC_022432.1) and mitochondrial (NC_022796.1) genomes were downloaded from NCBI. Nanopore reads were filtered to exclude those <1000 bp using filtlong, then mapped to remove host reads using minimap2 [37]. These reads were then used to generate a unicycler Illumina-nanopore co-assembly from all non-host reads.

2.2. Examination of Previously Reported Phytoplasma Genomes

To retrieve all phytoplasma genomes from GenBank, the NCBI Genomes database (https://www.ncbi.nlm.nih.gov/datasets/genome/) was queried using the search term, “phytoplasma”. This search (2025/01/01) retrieved 271 genomes, from which a single example of each unique taxID was selected. In cases where multiple taxID were represented, the genome with the highest completeness score reported was selected. Genomes with a reported CheckM completeness score below 90% were excluded from analysis. This resulted in 65 phytoplasma genomes in the final dataset (Table S1).

Genomes corresponding specifically to group 16SrIII were downloaded from NCBI using the information provided by Fernandez et al. [38] (Table 1). This ensured that no currently reported genomes within this group were missed in our analysis.

Genes annotated as groEL, groES, and 16S rRNA were selected from each genome outside of group 16SrIII. In addition, taxonomic markers secY, secA, tuf, and nusA [39] were retrieved from Peanut Witches’ Broom phytoplasma (PnWB; GCA_000364425.1) and from selected 16SrIII genomes.

Table 1. ‘Candidatus Phytoplasma pruni’ group 16SrIII genomes included in the analysis reported in this work. This table is based on and modified from Fernandez et al. 2024 [38].

Phytoplasma	strain	16SrIII- subgroup	Genome level	Host	location	accession	genome size, kb	Reference
Cicuta Witches broom	CicWB	16SrIII-J	contig (16)	Conium maculatum	Argentina	GCA_035853675.1	758	[38]
Phytoplasma Vc33	Vc33	16SrIII-J	contig (36)	Catharanthus roseus	Chile	GCA_001623385.2	687	[40]
China tree decline	ChTDIII	16SrIII-B	contig (67)	Melia azedarach	Argentina	GCA_013391955.1	791	[41]
Italian clover phyllody	MA	16SrIII-B	contig (197)	Catharanthus roseus	Italy	GCA_000300695.1	597	[21]
Poinsettia branch-inducing	JR	16SrIII-A	contig (185)	Euphorbia pulcherrima	USA	GCA_000309465.1	631	[21]
Milkweed Yellows Phytoplasma	MW	16SrIII-F	contig (158)	Catharanthus roseus	Canada	GCA_000309485	584	[21]
Vaccinium Witches’ Broom	VAC	16SrIII-F	contig (272)	Vaccinium myrtillus	Italy	GCA_000309405.1	648	[21]
‘Candidatus Phytoplasma pruni’	PR2021	16SrIII-A	chromosome (1)	Euphorbia pulcherrima	Taiwan	GCA_029746895.1	710	[42]
‘Candidatus Phytoplasma pruni’	CX	16SrIII-A	contig (46)	Catharanthus roseus	USA	GCA_001277135.1	599	[43]
Milkweed Yellows Phytoplasma	MYp- CanS4	16SrIII-F	contig (7)	Asclepias syriaca	Canada	GCA_050286905.1	694	this work
‘Candidatus Phytoplasma pruni’	2A1	16SrIII-A	contig (2)	lilac	Canada	GCA_033391615.1	625	this work

Table 1. ‘Candidatus Phytoplasma pruni’ group 16SrIII genomes included in the analysis reported in this work. This table is based on and modified from Fernandez et al. 2024 [38].

Phytoplasma	strain	16SrIII- subgroup	Genome level	Host	location	accession	genome size, kb	Reference
Cicuta Witches broom	CicWB	16SrIII-J	contig (16)	Conium maculatum	Argentina	GCA_035853675.1	758	[38]
Phytoplasma Vc33	Vc33	16SrIII-J	contig (36)	Catharanthus roseus	Chile	GCA_001623385.2	687	[40]
China tree decline	ChTDIII	16SrIII-B	contig (67)	Melia azedarach	Argentina	GCA_013391955.1	791	[41]
Italian clover phyllody	MA	16SrIII-B	contig (197)	Catharanthus roseus	Italy	GCA_000300695.1	597	[21]
Poinsettia branch-inducing	JR	16SrIII-A	contig (185)	Euphorbia pulcherrima	USA	GCA_000309465.1	631	[21]
Milkweed Yellows Phytoplasma	MW	16SrIII-F	contig (158)	Catharanthus roseus	Canada	GCA_000309485	584	[21]
Vaccinium Witches’ Broom	VAC	16SrIII-F	contig (272)	Vaccinium myrtillus	Italy	GCA_000309405.1	648	[21]
‘Candidatus Phytoplasma pruni’	PR2021	16SrIII-A	chromosome (1)	Euphorbia pulcherrima	Taiwan	GCA_029746895.1	710	[42]
‘Candidatus Phytoplasma pruni’	CX	16SrIII-A	contig (46)	Catharanthus roseus	USA	GCA_001277135.1	599	[43]
Milkweed Yellows Phytoplasma	MYp- CanS4	16SrIII-F	contig (7)	Asclepias syriaca	Canada	GCA_050286905.1	694	this work
‘Candidatus Phytoplasma pruni’	2A1	16SrIII-A	contig (2)	lilac	Canada	GCA_033391615.1	625	this work

2.3. Examination of 16SrIII Genomes for Sequences Related to groEL and groES

Each of the downloaded genomes was used to prepare a local BLAST database. The 16SrIII genome 2A1 was queried using tBLASTn, using a GroEL amino acid sequence from ‘Ca. P. asteris’ (strain AYWB; cpnDB ID b8392; 16SrI) as input. This retrieved a short DNA sequence (165 bp) that encoded a predicted protein of 55 amino acids with a BLASTp match to phytoplasma GroEL. This sequence was then used to query all other 16SrIII genomes using BLASTn. This resulted in significant matches in some genomes, but only short, nonspecific matches (under 20 bp) in others. The predicted amino acid sequences of each of the matching regions was compared to the GenBank database using BLASTp, and those with a strongly significant match to GroEL (E values between 10^-5 and 10^-50) were considered to be potential pseudogenes. A similar approach was used for groES.

2.4. Examination of the Gene Neighborhood of groEL in Group 16SrII and 16SrIII Phytoplasmas

The annotated genomes of ‘Ca. P. pruni’ strain CX (Ga0100078) and ‘Ca. P. aurantifolia strain PnWB (Ga0248296) at the Genomes Online Database (https://gold.jgi.doe.gov/) were used to examine the gene neighborhood of the intact GroE system in PnWB compared to the putative pseudogene identified in strain CX. Genes immediately downstream of the annotated groEL gene or gene fragment from each strain were downloaded and used to identify the corresponding genes in each of the group III phytoplasmas that were suspected to harbor groEL pseudogenes. The locations and annotations of the genes corresponding to the immediate gene neighbors of the GroE system in PnWB were noted in each strain, and intergenic distances calculated. Both gene and predicted amino acid sequences of these genes were downloaded and used for sequence similarity determination.

3. Results and Discussion

3.1. Genome Sequencing of ‘Ca. P. pruni’ Strains 2A1 and MYp-CanS4

The genome sequence of strain 2A1 was 625 kb in length (Table 1), with a G/C content of 27.32%. The genome consisted of two scaffolds with lengths of 336 kb and 289 kb. 594 protein-coding genes were predicted in the genome, 437 of which had function prediction. In addition, the genome contained two copies of 16S rRNA, as is typical for phytoplasmas [44]. Both genes typed at the iPhyclassifier as 16SrIII-A (F=1.0), indicating a lack of the 16S rRNA gene heterogeneity that has been observed in some phytoplasmas in the 16SrIII group [45].

Symptomatic milkweed at the Canadensis Botanical Garden was confirmed to be infected by phytoplasma using qPCR targeting 16S rRNA genes. The 16S rRNA genes were amplified by PCR, sequenced, and examined using RFLP analysis. This showed that strain MYp-CanS4 typed as 16SrIII-F (F=1.0), again with no evidence of gene heterogeneity. This is consistent with previous RFLP typing results of milkweed yellows phytoplasma [21]. The geographic location of the MYp detected in this study is consistent with the original report of MYp in Tichbourne, Ontario in the early 1990s [46], indicating that this phytoplasma has been circulating in Eastern Canada for many years.

The total length of the MYp-CanS4 genome was 694 kb, which is considerably longer than the 584 kb that was previously reported (Table 1). The genome was assembled into 7 contigs with an N50 of 388 kb, which is also an improvement over the N50 of 8 kb in the previously reported MYp genome [21]. The MYp-CanS4 genome, like all phytoplasma genomes, featured a low G/C percentage (26.7%). Annotation of the MYp-CanS4 genome using PGAP revealed the presence of 718 genes in total, including 678 coding sequences, 6 rRNA genes (2 each of 5S, 16S, 23S), and 24 pseudogenes. Both copies of the 16S rRNA-encoding genes typed at the iPhyClassifier as 16SrIII-F, consistent with the PCR results and with the known classification of MYp as a 16SrIII-F phytoplasma [47]. A plasmid of 4.4 kb was predicted by plasmer (contig 6). In accordance with this prediction, BLAST analysis of this contig showed strong similarity (93% sequence identity, e = 0.0) to the annotated plasmid of strain PR2021, pPR2021. Examination of the completeness of all the reported 16SrIII genomes from Table 1 using Busco revealed that the highest completeness score, 67.2%, was observed for the genomes PR2021 (the only chromosome-level assembly), MYp-CanS4, 2A1, and ChTDIII. In contrast, the previously reported MYp assembly had a completeness score of 66.4%, and the lowest score (59.5%) was observed for Vc33. These relatively low Busco completeness scores are due to the reduced genomes of phytoplasmas, which are missing some of the reference genes used by Busco to determine completeness.

3.2. Examination of Phytoplasma Genomes for GroE Chaperonin System

In total, 65 unique taxonomic IDs corresponding to a wide variety of phytoplasmas were examined (Table S1). Phylogenetic analysis of the 16S rRNA genes from these phytoplasmas provided results that were consistent with previous 16S gene-based analysis [11,48], with 3 basal clades represented. One major group consisted mostly of 16Sr groups I and XII, and this group was more closely related to the Acholeplasma taxa compared to the other phytoplasma strains. A second group consisted of a diverse range of phytoplasma taxa, with a complex branching pattern. The third major clade consisted of two main branches, including the 16SrIII group (‘Ca. P. pruni’), and a second branch consisting primarily of group 16SrII strains, along with ‘Ca. P. melaleucae’ (16Sr XXV-A) (Figure 1). The 16SrIII group was most closely related to the 16SrII group, as reported previously [48].

All but 6 of the genomes examined (59) had intact GroE systems, with both groEL and groES found in each genome (Table S2). All of the genomes that lacked a GroE system belonged to the 16SrIII group. The genomes with intact groEL/groES included a wide taxonomic variety, with 31 different phytoplasma species and 16 distinct 16Sr groups represented (Table S2). This analysis indicates that, among all high-quality phytoplasma genomes that have been sequenced to date, ‘Ca. P. pruni’ (16SrIII) is unique in lacking a GroE system. All other phytoplasma taxonomic groups possess both groEL and groES, which are located next to each other in their respective genomes in all cases. While it was previously noted that four phytoplasma genomes from 16SrIII lack a GroE system [21], this observation can now be extended to 11 genomes. Furthermore, while the GroE system was known to be well represented within phytoplasma groups other than 16SrIII [49,50,51,52], it was previously not known if 16SrIII is indeed the only exception. The results shown in Table S1 provide strong evidence that the GroE system is very widely distributed in phytoplasmas other than 16SrIII. Phylogenetic analysis of full-length groEL sequences representing all of the phylogenetic diversity that is currently present in public databases (NCBI) is shown in Figure 2. This analysis demonstrates the very wide sequence divergence of this phylogenetic marker, even in strains that are relatively closely related, underscoring the utility of this marker for differentiation of phytoplasma strains outside of the 16SrIII group [52].

3.3. Pseudogenization of groEL in 16SrIII Phytoplasmas

tBLASTn analysis of the genomes of all group III phytoplasmas showed that, consistent with the initial report, most of the sequenced genomes from this group contained no coding sequences with any significant similarity to amino acid sequences of phytoplasma GroEL. However, 5 of the 11 genomes did contain nucleotide sequences encoding predicted amino acid sequences that were strikingly similar to GroEL from other phytoplasmas (Table 2). These amino acid sequences were short (50, 55, or 62 residues) and could not encode a functional protein. One of these putative pseudogenes, encoded on the genome of ‘Ca. P. pruni’ strain CX (GCA_001277135.1) was annotated as a pseudogene (labelled “hypothetical protein” by the Prokaryotic Genome Annotation Pipeline at GenBank) [35], while the others were identified by local custom BLAST databases. Alignment of these predicted amino acid sequences to the 546 amino acid sequence of GroEL from PnWB phytoplasma (16SrII) revealed that the gene fragments encoded amino acids near the carboxy terminus of the protein (Figure 3). Phylogenetic analysis of the 16SrIII phytoplasma strains based on the sequences of five taxonomic markers (secY; secA; nusA; tuf; rp) revealed that the three types of identified pseudogenes clustered consistently with the phylogenetic placement of the phytoplasmas containing them (Figure 4). The putative pseudogenes encoding these amino acids are considered to be unitary pseudogenes, in that there is no evidence of gene duplication in the genome, so that the loss of these genes is associated with a complete loss of function [55]. No evidence of nucleic acid sequences encoding amino acids with similarity to GroES was found in any of the genomes.

3.4. Conserved Synteny of GroE System in 16SrII and Pseudogenes of 16SrIII

To determine if the pseudogenes were orthologous by synteny, rather than simply a case of coincidental similarity [55], the gene neighborhood of the intact GroE system in a close phylogenetic neighbor of 16SrIII, ‘Ca. P. aurantifolia’ (PnWB; 16SrII) was examined and compared to the locations of the putative pseudogenes.

The GroE system (groEL/groES) in PnWB is located immediately adjacent to two genes that are annotated as “multidrug resistance ABC transporter ATP-binding and permease protein” (Figure 5). These genes are identified as evbH (contig 9: 24713-26554, reverse strand) and evbG (contig9: 26551-28284, reverse strand). evbG is located 216 nucleotides downstream of groEL (contig9: 28499-30136, reverse strand), while groES is located 542 nucleotides upstream of groEL (contig9: 30677-30961, reverse strand) (Table 2).

In the 16SrIII phytoplasmas that contain putative groEL pseudogenes, there are two genes located immediately next to each pseudogene (Figure 5). These protein-coding genes are annotated variously at GenBank as “multidrug resistance ABC transporter ATP-binding and permease protein” (strain PR2021), or, “ABC transporter ATP-binding protein” (strains 2A1, CX, MA-ICP). These genes encode proteins that are similar in size to PnWB evbG and evbH – 576/577 amino acids for evbG, and 593/613 amino acids for evbH (Table 2). In each 16SrIII genome, the putative pseudogene is located immediately beside (159-429 nucleotides) the gene that is similar in size to evbG of PnWB, suggesting a similar gene neighborhood of the putative pseudogene in 16SrIII compared to the intact groEL gene in PnWB. (Table 2; Figure 5).

Table 2. Amino acid lengths and genome coordinates of groEL pseudogenes and genes immediately upstream in 16SrII and 16SrIII.

Phytoplasma

strain

16SrIII-subgroup

groEL (pseudo)gene length, bp

Amino acids

Coordinates (contig:bases)

evbG ortholog length, bp

Amino acids

Coordinates (contig:bases)

groEL-evbG distance, bp

evbH ortholog length, bp

Amino acids

Coordinates (contig:bases)

Cicuta Witches broom

CicWB

16SrIII-J

none

Phytoplasma Vc33

Vc33

16SrIII-J

none

China tree decline

ChTDIII

16SrIII-B

none

Italian clover phyllody

MA

16SrIII-B

153

50

33:3851-4003

1728

576

33:1856-3586

266

1782

593

33:78-1859

Poinsettia branch-inducing

JR

16SrIII-A

186

62

171:1963-2148

17031

5661

5:13448-15150

ND1

1782

593

5:11670-13451

Milkweed Yellows

MYp-CanS4

16SrIII-F

none

Vaccinium Witches’ Broom

VAC

16SrIII-F

none

‘Candidatus Phytoplasma pruni’

PR2021

16SrIII-A

186

62

686562-686747

1728

576

684481-686211

352

1782

593

682703-684484

‘Candidatus Phytoplasma pruni’

CX

16SrIII-A

165

55

8:22632-22796

1728

576

8:20474-22204

429

1782

593

8:18696-20477

‘Candidatus Phytoplasma pruni’

2A1

16SrIII-A

165

55

1:56157-56321

1728

576

1:56749-58479

159

1782

593

1:58476-60257

Peanut witches' broom phytoplasma

PnWB

16SrII-A

1638

546

9:28499-30136

1731

577

9:26551-28284

216

1842

613

9:24713-26554

1minus strand - contig 5 length is 15150 bp. Therefore, the junction of evbG and groEL pseudogene is disrupted and intergenic distance cannot be calculated

Table 2. Amino acid lengths and genome coordinates of groEL pseudogenes and genes immediately upstream in 16SrII and 16SrIII.

Phytoplasma

strain

16SrIII-subgroup

groEL (pseudo)gene length, bp

Amino acids

Coordinates (contig:bases)

evbG ortholog length, bp

Amino acids

Coordinates (contig:bases)

groEL-evbG distance, bp

evbH ortholog length, bp

Amino acids

Coordinates (contig:bases)

Cicuta Witches broom

CicWB

16SrIII-J

none

Phytoplasma Vc33

Vc33

16SrIII-J

none

China tree decline

ChTDIII

16SrIII-B

none

Italian clover phyllody

MA

16SrIII-B

153

50

33:3851-4003

1728

576

33:1856-3586

266

1782

593

33:78-1859

Poinsettia branch-inducing

JR

16SrIII-A

186

62

171:1963-2148

17031

5661

5:13448-15150

ND1

1782

593

5:11670-13451

Milkweed Yellows

MYp-CanS4

16SrIII-F

none

Vaccinium Witches’ Broom

VAC

16SrIII-F

none

‘Candidatus Phytoplasma pruni’

PR2021

16SrIII-A

186

62

686562-686747

1728

576

684481-686211

352

1782

593

682703-684484

‘Candidatus Phytoplasma pruni’

CX

16SrIII-A

165

55

8:22632-22796

1728

576

8:20474-22204

429

1782

593

8:18696-20477

‘Candidatus Phytoplasma pruni’

2A1

16SrIII-A

165

55

1:56157-56321

1728

576

1:56749-58479

159

1782

593

1:58476-60257

Peanut witches' broom phytoplasma

PnWB

16SrII-A

1638

546

9:28499-30136

1731

577

9:26551-28284

216

1842

613

9:24713-26554

1minus strand - contig 5 length is 15150 bp. Therefore, the junction of evbG and groEL pseudogene is disrupted and intergenic distance cannot be calculated

To determine if these similar genes are likely orthologs, we compared the predicted amino acid sequence similarity of these genes in PnWB and the 16SrIII strains. Within the 16SrIII genomes, the two genes immediately adjacent to the putative groEL pseudogenes had amino acid similarity scores (by clustalw) of 96-99%, suggesting that they are orthologous genes within these strains (Figure 5, Table S2). The amino acid similarity scores of these 16SrIII genes to the PnWB evbG and evbH were 47% and 48% (Figure 5, Table S2). These scores are similar to the amino acid similarities of the orthologous taxonomic markers secY, secA, nusA, and tuf in the 16SrII and 16SrIII genomes, which ranged from 38-56% (Table S3). In contrast, when all other genes in the PnWB chromosome with annotations containing the phrase “ABC transporter ATP-binding” were compared to the genes immediately upstream of the putative groEL pseudogenes in 16SrIII, the amino acid similarity scores were much lower – 15-20% (data not shown). These amino acid sequence similarity comparisons suggest that the two genes encoding ABC transporters in 16SrIII are indeed orthologs of evbG and evbH in PnWB, confirming the syntenic relationship of groEL genes and pseudogenes in these closely related phytoplasmas.

The genomic region upstream of the groEL pseudogenes in 16SrIII strain CX contained no evidence of the shorter groES genes, and displayed a dearth of annotated genes in general compared to the corresponding gene region in PnWB. Only a few short hypothetical proteins are annotated in the 16SrIII genome, on the opposite strand of the groEL pseudogenes and evbG/evbH (Figure S1). The reason for this genomic "desert” is unclear, but this area seems to have been subjected to an ancient pseudogenization that resulted in the loss of the GroE system, with only a small vestige of the groEL gene remaining in this region in extant strains. This was not the case for the 16SrIII strain 2A1, but there was no evidence of a gene with any similarity to groES in the area upstream of the groEL pseudogene (Figure S1).

3.5. No Evidence of HGT in Non-Fragmented 16SrIII Genomes

Due to the clear cases of groEL acquisition by HGT that were previously documented in Mollicutes [11,12], we examined 16SrIII genomes for the presence of any groEL genes, even those more distantly related to phytoplasmas. No evidence of HGT was observed in any of the genomes. Many of the previously reported 16SrIII genomes are highly fragmented (more than 150 contigs) based on short-read but deep and accurate sequencing chemistry (Illumina). In addition, the assembly process commonly involves a mapping step to exclude non-phytoplasma DNA because of the metagenomic nature of phytoplasma DNA samples. To address this possible source of error, we prepared a MYp assembly by only removing host reads (chromosomal, mitochondrial, and chloroplast), leaving all non-host reads (mostly bacterial, but possibly also including other host-associated microorganisms) available for assembly. Removal of host reads left 21,754,990/ 57,879,811 Illumina reads (37.6%), compared to 12,482,861 reads (21.6%) that remained when mapped to the previously reported MYp genome sequence. This suggests that many other microorganisms were associated with the milkweed leaf sample. When these host-depleted, assembled reads were queried for the presence of groEL-like sequences, 9 sequences were retrieved, all from bacteria (Table S2). Examination of the assembled contigs containing these groEL sequences (the chromosomal context) provided taxonomic identifications consistent with the groEL sequences in all cases (Table S4). These results exclude the possibility of a groEL gene being acquired by HGT, as its sequence would then be located within a context of phytoplasma-like DNA sequences.

Here we have provided two additional high-quality, non-fragmented assemblies of 16SrIII phytoplasmas, and we observed no evidence of groEL acquisition by HGT either in the genomes reported here or in previously published assemblies that are less fragmented (eg strain PR2021 is a single contiguous chromosome). Therefore, we conclude that 16SrIII genomes are unique among phytoplasmas in that they contain no functional GroE system.

4. Conclusions

GroEL/Cpn60 is a protein that is essential for cell viability due to its ability to interact with and facilitate the folding of a wide range of client proteins in the cell, including proteins of a variety of sizes (20-80 kDa) with functions in regulatory and structural processes [13]. The protein also possesses a very wide array of non-canonical functions, which are often associated with gene duplication events [13,56,57]. Under this scenario, preservation of the ancestral function of GroEL can be retained while the duplicated gene can accumulate mutations that may provide alternative functions, leading to functional adaptation of the protein [57]. While this model may partially explain the functional diversity of this essential protein, it does not explain how Mollicutes have adapted to genome reduction including the loss of the GroE system. Schwarz et al. concluded that groEL pseudogenization could be overcome by just a few amino acid changes within a limited range of client proteins that are essential for cell viability, implying that the viability barrier caused by the loss of groEL is not particularly high [11].

Nevertheless, a functional GroE system is the ancestral state for bacteria, given that all bacterial taxa except most of the Mollicutes possess it. Since the Acholeplasmataceae are considered to be a basal clade within the Mollicutes due to their phylogenetic position [1,58] and retention of the AUG stop codon [8], it is likely that the single-copy GroE system in the genera Phytoplasma and Acholeplasma represents a retention of the genes from the Mollicute progenitor, rather than a re-acquisition as has been speculated [11,12]. This is evident in the 16SrIII genomes that retain a vestige of the groEL gene, as its predicted amino acid sequence resembles that of other phytoplasmas. If the GroE system in Acholeplasmataceae represents a loss and re-acquisition, then the 16SrIII group would have lost, then re-gained, and lost this system again. While this is possible in the highly plastic genomes of the phytoplasmas [59], it seems unlikely. It is also notable that, while virtually all species of Acholeplasma have retained a GroE system, a single species, A. oculi, has lost these genes [11]. A single taxon within both of these genera seems to have lost the ancestral GroE system, while most members of this group have retained it.

Through genome sequencing and analysis of previously sequenced genomes, we have demonstrated that group 16SrIII phytoplasmas have lost the ancestral GroE system, which is retained in all other taxa of phytoplasmas. No evidence of the presence of groEL/groES genes by horizontal gene transfer was observed, indicating that these bacteria, like most other Mollicutes, have not re-acquired the functionality of GroE. Since groEL is absent in an obviously pathogenic group of phytoplasma, it seems unlikely that the protein is involved in pathogenesis as has been speculated from observations in other Mollicutes [12]. The lack of gene duplication within phytoplasma groEL genes also argues against a moonlighting role for the protein in host cell attachment and invasion and indicates that the protein most likely retains a canonical role in the folding of host cell proteins. Those Mollicutes that have lost the ancestral genes during genome reduction seem to have crossed the low evolutionary barrier identified by Schwarz et al. [11] and adapted their client proteins to be able to fold in the absence of this system.