Shu-Ting Cho
Hung-Jui Kung
Weijie Huang
Saskia A. Hogenhout
Chih-Horng Kuo- 1Institute of Plant and Microbial Biology, Academia Sinica, Taipei, Taiwan
- 2Department of Crop Genetics, John Innes Centre, Norwich, United Kingdom
Phytoplasmas are plant-pathogenic bacteria that impact agriculture worldwide. The commonly adopted classification system for phytoplasmas is based on the restriction fragment length polymorphism (RFLP) analysis of their 16S rRNA genes. With the increased availability of phytoplasma genome sequences, the classification system can now be refined. This work examined 11 strains in the 16SrI group within the genus ‘Candidatus Phytoplasma’ and investigated the possible species boundaries. We confirmed that the RFLP classification method is problematic due to intragenomic variation of the 16S rRNA genes and uneven weighing of different nucleotide positions. Importantly, our results based on the molecular phylogeny, differentiations in chromosomal segments and gene content, and divergence in homologous sequences, all supported that these strains may be classified into multiple operational taxonomic units (OTUs) equivalent to species. Strains assigned to the same OTU share >97% genome-wide average nucleotide identity (ANI) and >78% of their protein-coding genes. In comparison, strains assigned to different OTUs share < 94% ANI and < 75% of their genes. Reduction in homologous recombination between OTUs is one possible explanation for the discontinuity in genome similarities, and these findings supported the proposal that 95% ANI could serve as a cutoff for distinguishing species in bacteria. Additionally, critical examination of these results and the raw sequencing reads led to the identification of one genome that was presumably mis-assembled by combining two sequencing libraries built from phytoplasmas belonging to different OTUs. This finding provided a cautionary tale for working on uncultivated bacteria. Based on the new understanding of phytoplasma divergence and the current genome availability, we developed five molecular markers that could be used for multilocus sequence analysis (MLSA). By selecting markers that are short yet highly informative, and are distributed evenly across the chromosome, these markers provided a cost-effective system that is robust against recombination. Finally, examination of the effector gene distribution further confirmed the rapid gains and losses of these genes, as well as the involvement of potential mobile units (PMUs) in their molecular evolution. Future improvements on the taxon sampling of phytoplasma genomes will allow further expansions of similar analysis, and thus contribute to phytoplasma taxonomy and diagnostics.
Introduction
Phytoplasmas are a group of insect-transmitted plant-pathogenic bacteria that reduce yields of diverse crops worldwide (Lee et al., 2000; Hogenhout et al., 2008; Bertaccini and Lee, 2018). These obligate parasites are related to the animal-pathogenic mycoplasmas (Chen et al., 2012), and both groups of these wall-less bacteria were assigned to the class Mollicutes under the phylum Tenericutes (Brown, 2010). However, unlike mycoplasmas and several other Mollicutes lineages, pure culture of phytoplasmas outside of their hosts has not been achieved yet. Due to this reason, all phytoplasmas were assigned to a “Candidatus” (“Ca.”) genus (The IRPCM Phytoplasma/Spiroplasma Working Team, 2004).
For the purposes of pathogen detection or research on their basic biology, a reliable system for strain identification and classification is critical. For example, it is important to understand if the diseases affecting different plants, transmitted by different vectors, and/or occurring in different locations are caused by the same or different phytoplasma strains. In addition, to infer the evolutionary history and biodiversity of phytoplasmas, defining species or other equivalent biological units are of fundamental importance. The commonly used system for phytoplasma classification was first established in the 1990s (Lee et al., 1993, 1998). This system utilizes a defined set of 17 restriction enzymes to perform restriction fragment length polymorphism (RFLP) analysis on a 1.25 kb PCR product derived from their 16S rRNA gene (Gundersen and Lee, 1996). Under this system, a set of reference strains were chosen, and a new strain could be assigned to the same 16S rRNA gene RFLP (16Sr) group and subgroup as a reference when the similarity coefficient is higher than 0.85 and 0.97, respectively (Zhao et al., 2009; Zhao and Davis, 2016). With the development of an easy-to-use web-based tool iPhyClassifier (Zhao et al., 2009), this system has been well-adopted by the phytoplasma research community.
However, despite the popularity, this 16Sr classification system has several shortcomings. To begin with, the RFLP approach considers only those nucleotides located within the restriction sites, rather than the entire sequence. Additionally, the similarity coefficients of the RFLP patterns do not provide information about phylogenetic relationships. Furthermore, intragenomic sequence variations of 16S rRNA genes may cause problems in classification. Most importantly, although 16S rRNA gene is considered as a universal marker for bacterial classification, the use of a single locus is inherently unreliable. Rather, adoption of multilocus sequence analysis (MLSA) is strongly recommended (Glaeser and Kämpfer, 2015). Unfortunately, selection of suitable markers for MLSA is not straightforward. Due to the highly variable evolutionary rates among different bacteria and among different genes (Kuo and Ochman, 2009), the markers that work well for one group do not necessary work well for others. Although extensive efforts have been devoted to the development of MLSA markers for phytoplasmas (Martini et al., 2019), genome-level assessment for the relative performance of these markers is lacking. Moreover, while the use of genome analysis to replace MLSA provides a potential solution and is gaining popularity (Konstantinidis and Tiedje, 2005; Konstantinidis et al., 2017; Jain et al., 2018; Parks et al., 2018), genome analysis is still challenging for uncultured bacteria such as phytoplasmas. For these reasons, development and evaluation of MLSA markers are still important.
To develop and evaluate MLSA markers, the availability of genome sequences provides a comprehensive guide for PCR primer design. For example, based on the complete genome sequence of “Candidatus Phytoplasma asteris” OY-M (Oshima et al., 2004), 18 PCR primer sets were designed to cover different regions of the chromosome (Kakizawa and Kamagata, 2014). Moreover, by developing a multiplex-PCR method for these primer sets and establishing the amplification patterns in a collection of reference strains, this previous study provided a fast and cost-effective method for strain identification (Kakizawa and Kamagata, 2014). As the number of available genome sequences increases, the genome-enabled MLSA maker development could be further improved in two critical aspects. First, the genome-scale molecular phylogeny provides a reliable reference for evaluating the performance of each candidate markers, such that the best-performing ones could be selected. Second, with the comprehensive genomic information of all reference strains, the chromosomal location of each candidate markers and the exact primer sequences could be evaluated using bioinformatic methods. In this regard, phytoplasmas belonging to the 16SrI group, also often known as the aster yellows (AY) group or “Ca. P. asteris” (i.e., a provisional species defined as encompassing all known subgroups within the 16SrI group) (Lee et al., 2004), represent the best system to test the concept of genome-assisted MLSA marker development. Among the >30 16Sr groups that have been described (Zhao and Davis, 2016), only eight have genome sequences available (Cho et al., 2019b). Moreover, the genome sequencing efforts have been highly focused on the 16SrI group (Table 1), which has a worldwide geographical distribution, wide range of plant hosts, and great impact on agriculture (Lee et al., 2004). Finally, this group is also the one that received most research attention on the molecular mechanisms of plant-microbe interactions; most of the phytoplasma effectors that have been characterized were initially studied in 16SrI phytoplasmas (Bai et al., 2009; Hoshi et al., 2009; MacLean et al., 2011; Sugio et al., 2011b; Tomkins et al., 2018; Wang et al., 2018b; Huang and Hogenhout, 2019). Previous studies of this group identified extensive levels of genome divergence and suggested that “Ca. P. asteris” may be classified into multiple operational taxonomic units (OTUs) equivalent to species (Firrao et al., 2013; Cho et al., 2019a). In this work, we expanded the sampling of available genome sequences to perform in-depth analysis of genome comparisons. Specifically, we aimed to provide quantitative guidelines for defining putative species boundaries, which could better inform phytoplasma taxonomy. Additionally, we aimed to develop MLSA markers that could facilitate future genetic characterization of the phytoplasmas belonging to the 16SrI or AY group.
Table 1. List of the genome sequences analyzed.
Materials and Methods
Data Sets and Analysis Methods
The 11 genome sequences included in this study are listed in Table 1. These data sets represent all phytoplasma genomes available from GenBank (Benson et al., 2018) as of January 2020 and are recognized as belonging to the 16SrI group, AY group, or “Ca. P. asteris” (Lee et al., 2004). Two other “Ca. Phytoplasma” species have been described as being affiliated with the 16SrI group, including “Ca. P. japonicum” (Sawayanagi et al., 1999) and “Ca. P. lycopersici” (Arocha et al., 2007), but were not included in this study due to the lack of genome sequence availability. Notably, “Ca. P. japonicum” was first reported as a member of the 16SrI-D subgroup (Sawayanagi et al., 1999) and later listed under the 16SrI group in the current taxonomy (The IRPCM Phytoplasma/Spiroplasma Working Team, 2004). However, the iPhyClassifier (Zhao et al., 2009) analysis of its 16S rRNA gene sequence (GenBank accession AB010425) assigned it to the 16SrXII-D subgroup, suggesting that its classification remained to be investigated.
To confirm the 16Sr group and subgroup assignments, we extracted the 16S rRNA gene sequences from the genomes and analyzed those sequences using the iPhyClassifier (Zhao et al., 2009). Most phytoplasma genomes contain two copies of 16S rRNA genes and each copy was processed individually. The strain OY-V lacks any 16S rRNA gene in the current assembly and therefore was excluded from this analysis. The strain NJAY has four partial sequences for its 16S rRNA genes and all of these four sequences were examined.
All bioinformatic tools for processing these data sets were used with the default settings unless stated otherwise. For correlation tests, the cor.test function implemented in R v3.4.4 (R Core Team, 2019) was used to examine the Pearson’s product moment correlation coefficient.
For strain TW1 (GenBank genome accession GCA_003181115.1), we downloaded the raw sequencing results from the NCBI Sequence Read Archive (SRA) to examine the phytoplasma sequences present in the sequencing libraries. The raw reads were mapped to two other reference genomes (i.e., AYWB and OY-M; see Table 1) to detect the presence of 16SrI-A and 16SrI-B type phytoplasmas, respectively. For the Oxford Nanopore MinION reads (SRA accession SRR7548026), the mapping was performed using Minimap2 v2.15 (Li, 2018). For the Illumina MiSeq paired-end reads (SRA accession SRR7548027), the mapping was performed using the Burrows-Wheeler Alignment (BWA) tool v0.7.17 (Li and Durbin, 2009). The mapping results were programmatically checked using SAMtools v1.9 (Li et al., 2009) to identify the sequence variations.
Genome Comparisons
For whole-genome comparison, FastANI v1.1 (Jain et al., 2018) was used to calculate the proportion of genomic segments mapped and the average nucleotide identity (ANI) of those segments in each genome-pair. For each pair, reciprocal comparisons between the two genomes were conducted (i.e., X as the query against Y as the reference and Y as the query against X as the reference). To calculate the proportion of genomic segments mapped, the numbers of segments mapped from those two reciprocal comparisons were combined as the numerator, and the total numbers of segments from those two genomes were combined as the denominator. The reciprocal ANI values were averaged. For example, in the AYWB-NJAY comparison, 183 out of the 239 segments of AYWB could be mapped to NJAY, and all 187 segments of NJAY could be mapped to AYWB. Based on these results, these two strains share (183 + 187) / (239 + 187) = 370 / 426 = 86.9% of their genomic segments. The ANI values were 99.5 and 99.8% in the two reciprocal comparisons, and an average ANI of 99.7% was reported.
For gene-centric investigation, the homologous gene clusters among all genomes were identified using OrthoMCL v1.3 (Li et al., 2003) based on the procedures described in our previous studies (Chung et al., 2013; Lo et al., 2013). Briefly, all coding sequences (CDS), including those annotated as pseudogenes, were included in the all-against-all similarity searches using BLASTN v2.6.0 (Camacho et al., 2009) with an e-value cutoff of 1e–15. The BLASTN results were used as the input for OrthoMCL, which normalized the within- and between-genome comparisons and identified homologs according to the principle of reciprocal best hits. Similar to the procedures for whole-genome comparisons, the proportion of homologs shared and the sequence identity of shared homologs were considered as two separate metrics of genome similarities. The proportion of homologs shared was calculated based on averaging (i.e., the number of homologs shared was used as the numerator, and the average number of homologs present in those two genomes was used as the denominator). The sequence identities of coding regions were calculated based on multiple sequence alignment of those single-copy genes present in all of the genomes compared. The homologs were aligned using MUSCLE (Edgar, 2004) and the sequence identities were calculated using the DNADIST program in PHYLIP v3.697 (Felsenstein, 1989); as expected, reciprocal comparisons all produced the same result for sequence identities.
In addition to the pairwise comparisons, the gene content of all genomes was compared using a procedure based on principal coordinates analysis (Lo et al., 2018). For this procedure, the homologous gene clustering result was converted into a matrix of 11 genomes by 1,233 gene clusters; the value in each cell indicates the gene copy number. This matrix was converted into a Jaccard distance matrix among genomes using the VEGAN package in R, then processed using the PCOA function in the APE package (Popescu et al., 2012).
To compare the genes for putative secreted proteins among these phytoplasmas, we examined all coding sequences using a procedure described in a previous study (Cho et al., 2019b). Briefly, SignalP v5.0 (Armenteros et al., 2019) was used to check the presence of signal peptide based on the Gram-positive bacteria model. The program TMHMM v2.0 (Krogh et al., 2001) was used to verify that these putative secreted proteins have no transmembrane domain. For those four experimentally characterized effectors (i.e., SAP05, SAP11/SWP1, SAP54/PHYL1, and TENGU), all putative homologs (i.e., those belong to the same homologous gene cluster as defined by OrthoMCL or have a BLASTN e-value of lower than 1e–100) were manually added to the list of regardless of the SignalP prediction results.
The genome alignment of representative strains with high quality assemblies were performed using genoPlotR v0.8.9 (Guy et al., 2010). The syntenic regions were identified using BLASTN v2.6.0 (Camacho et al., 2009), the cutoff values for sequence similarity were set to e-value = 1e–15 and alignment length = 2,000 bp.
Molecular Phylogenetics
For phylogenetic inference, MUSCLE v3.8.31 (Edgar, 2004) was used to generate multiple sequence alignments and PhyML v3.3 (Guindon and Gascuel, 2003) was used to infer maximum likelihood phylogenies. The visualization was performed using JalView v2.11 (Waterhouse et al., 2009) and FigTree v1.4.4. In each inference, PHYLIP v3.697 (Felsenstein, 1989) was used to generate 1,000 replicates for bootstrap analysis.
Development of Molecular Markers
To identify molecular markers that may be suitable for classification, those multiple sequence alignments of single-copy genes shared by all genomes were examined manually. Genes were ranked by the number of variable sites that could be used for distinguishing different OTUs. Subsequently, the list was prioritized based on the criterion of having multiple markers located in distinct regions of those phytoplasma chromosomes. The PCR primers (Table 2) were designed manually with the aim of producing ∼700–800 bp products to cover those variable regions, such that each product could be sequenced in one single Sanger sequencing reaction.
Table 2. List of molecular markers for phytoplasma classification.
To validate the primer design, a DNA sample derived from a periwinkle infected by the phytoplasma strain DY2014 (Cho et al., 2019b) was used as the template for PCR tests. All kits were used according to the manufacturer’s instructions. The PCR were performed using the KAPA HiFi HotStart ReadyMix (Roche, United States). The cycling conditions were: (1) an initial activation step at 95°C for 3 min; (2) 35 cycles of 98°C for 20 s, 55°C for 30 s, and 72°C for 60 s; and (3) a final extension step of 72°C for 5 min. The annealing temperature was lowered to 53°C for the newly designed primers for secY (Table 2), while the other four primer sets used 55°C. Agarose gel electrophoresis was conducted to confirm that all PCR tests produced one single band of expected size. For further confirmation, the PCR products were purified using the QIAquick Gel Extraction Kit (Qiagen, Germany), then sequenced using the BigDye Terminator v3.1 Cycle Sequencing Kit on an Applied Biosystems 3730XL DNA Analyzer (Thermo Fisher Scientific, United States). The sequencing results were compared to the published genome sequence for final confirmation.
The procedure for calculating the substitution rates was based on that described in our previous studies (Kuo and Ochman, 2009; Kuo et al., 2009). Briefly, the homologous protein sequences were aligned by using MUSCLE v3.8.31 (Edgar, 2004). The protein alignments were converted into codon-based nucleotide alignment using PAL2NAL v14 (Suyama et al., 2006), which were then processed using the YN00 method (Yang and Nielsen, 2000) implemented in PAML v4.9h (Yang, 2007). To identify the homologs of selected marker genes in other more divergent phytoplasmas, the full gene sequences were used as queries to run BLASTN searches against available genomes. After the homolog identification, multiple sequence alignments were performed using MUSCLE (Edgar, 2004) and the newly designed primers (Table 2) were examined manually to check if these primers could work on other phytoplasmas.
Results and Discussion
Genome Characteristics
Among the 11 phytoplasma genomes included in this study, four have the complete chromosomal sequence available while the other seven are draft assemblies with varying levels of completeness (Table 1 and Figure 1). For those incomplete genomes, the true genome sizes are difficult to determine. Nonetheless, the largest (i.e., OY-M; 853 kb) and the smallest (i.e., MBSP-M3; 576 kb) genomes in this data set are both completed assemblies, confirming the high level of genome size variation among these closely related strains. The total number of CDS is strongly correlated with the genome size (r = 0.92, p = 2.7e–5). However, the number and proportion of pseudogenes varied widely among these genomes. For example, the strains AYWB and NJAY have a genome-wide ANI of 99.7% (Figure 2A), yet the proportion of pseudogenes are 17 and 41% (Table 1), respectively. In another case, DY2014 and OY-V have a genome-wide ANI of 99.4%, yet the proportion of pseudogenes are 8 and 25%, respectively. Although these differences may reflect true biological processes, such as elevated mutation accumulations that have occurred in the evolutionary history of some strains, it is also possible that these differences may be explained by artifacts. Regarding the sequencing methods, AYWB was based on Sanger sequencing (Bai et al., 2006) and DY2014 was based on Illumina (Cho et al., 2019b), while NJAY was based on one single Roche/454 library (Sparks et al., 2018) and OY-V was based on one Illumina paired-end library and one Roche/454 GS FLX mate-pair library (Kakizawa et al., 2014). Due to the low GC-content of these genomes (i.e., ∼27–28%), it is plausible that sequencing errors in homopolymeric regions, particularly those errors originated from the 454 sequencing technology, may have contributed to the high numbers of pseudogenes with frameshift mutations in NJAY and OY-V. Moreover, assemblies that are more fragmented may also have more partial genes located at the ends of contigs. Due to these concerns, all pseudogenes were included in the gene-centric investigation of this study.
Figure 1. Maximum likelihood phylogenies. The numbers on the internal branches indicate the level of bootstrap support based on 1,000 resampling; only values ≥60% are shown. (A) 16S rRNA genes. The alignment contains 1,541 aligned nucleotide sites. The strain name, locus tag of the 16S rRNA gene (in parentheses), and the 16SrI subgroup classification (in square brackets) are labeled. A question mark ‘?’ in the subgroup classification indicates that the sequence is classified as a new subgroup, the existing subgroup with the highest similarity is provided. (B) Single-copy coding genes shared by all strains. The concatenated alignment contains 303 genes and 291,990 aligned nucleotide sites. (C) 16S rRNA gene and the five markers developed in this study; see Table 2 for detailed information. The concatenated alignment contains 4,995 aligned nucleotide sites.
Figure 2. Pairwise genome similarities. The strains with the complete genome sequences available are highlighted in bold. (A) Similarity scores based on whole genome comparisons. Above diagonal: average nucleotide identity (%); diagonal: genome size (kb); below diagonal: genomic segments mapped (%). (B) Similarity scores based on genic regions. Above diagonal: genic sequence identity (%); diagonal: number of homologous gene clusters; below diagonal: homologs shared (%).
Classification of the 16SrI Phytoplasmas Analyzed
The 16SrI subgroup assignments of these strains were mostly consistent with those reported in literature but there were a few surprises (Figure 1A). First, WBD was reported as a 16SrI-B strain (Chen et al., 2014). However, our examination of the sequence record deposited in GenBank revealed that the two copies of 16S rRNA genes in this genome differ by 4-bp, and iPhyClassifier (Zhao et al., 2009) assigned these two sequences to 16SrI-R and 16SrI-S, respectively. Results from molecular phylogeny also support that these two sequences form a clade that is distinct from other 16SrI-B sequences. Second, DY2014 has two 16S rRNA genes that differ by 1-bp. While one was assigned to 16SrI-B, the other was identified as being a new subgroup with a RFLP pattern most similar to 16SrI-B. Finally, TW1 has two sequences that differ by 12-bp, with one assigned to 16SrI-A and the other assigned to 16SrI-B (Town et al., 2018).
These results highlighted several issues of this RFLP-based classification system for phytoplasmas. First, the intra-genomic variation of 16S rRNA genes may result in a strain being assigned to different subgroups depending on which homolog was examined. This issue has been reported for other phytoplasmas (Liefting et al., 1996) and a three-letter subgroup designation system has been proposed (Wei et al., 2008). For example, a paulownia witches’-broom (PaWB) phytoplasma was first classified as a 16SrI-D strain (Lee et al., 1998) but later found to harbor a 16SrI-B type sequence as well. To accommodate these findings, the subgroup status of this PaWB was redesignated as 16SrI-(B/D)D, with the first two letters in parentheses denoting the types of those two 16S rRNA genes, and the third letter denoting its 16Sr subgroup assignment (Wei et al., 2008; Zhao et al., 2009). Although this three-letter system provides comprehensive information, it may be unnecessarily complex and confusing. Moreover, it is unclear which one of the two types should be chosen for the formal subgroup assignment. In the case of PaWB, it was assigned to 16SrI-D based on precedence. However, in the case of WBD (i.e., originally described as a 16SrI-B strain), it is unclear if it should be redesignated as 16SrI-(R/S)R or 16SrI-(R/S)S. Second, this RFLP-based system considers only the sequence variations in restriction sites. By ignoring the sequence variations in other positions, this approach could alleviate the problem of intra-genomic variations to a certain degree (e.g., the two sequences in AYWB differ by 2-bp, yet both were assigned to 16SrI-A). However, this approach wasted much of the information contained in the full sequence and the uneven weighing of nucleotide positions could introduce strong biases. The strain DY2014 represents one extreme example to illustrate such biases, with a single-bp difference that is sufficient to result in different subgroup assignments between the two 16S rRNA genes. Finally, a classification system that depends solely on the 16S rRNA genes while ignoring other regions of the genome does not provide sufficient resolution.
Regardless of the group/subgroup assignments of these 16S rRNA genes, all these sequences share >99% identity (Supplementary Table S1), which is above the general recommended thresholds for defining species. These thresholds include: (1) 97% identity that is commonly used for defining OTUs corresponding to species in microbiota surveys (Schloss and Handelsman, 2005), (2) 97.5% identity recommended for “Ca. Phytoplasma” species (The IRPCM Phytoplasma/Spiroplasma Working Team, 2004), and (3) 98.6% identity recommended for uncultivated bacteria (Konstantinidis et al., 2017). In other words, these strains may be considered as all belonging to “Ca. P. asteris” (Lee et al., 2004) based on their 16S rRNA gene sequences. However, delineating species-level OTUs based on one single gene may be problematic; more in-depth investigation of their genomic divergence is necessary.
Genomic Divergence Among 16SrI Phytoplasmas and Putative Species Boundaries
To further investigate the divergence of these phytoplasmas, we utilized two approaches to quantify their genomic similarities. In the first approach, the whole genome sequences were used to calculate the proportion of genomic segments shared, as well as the ANI values of these segments. A previous study that examined >90,000 prokaryotic genomes found that the ANI values exhibit a bimodal distribution, and the within-species comparisons almost always have >95% ANI (Jain et al., 2018). We found that based on this criterion, the 11 strains examined could be separated into three clusters that roughly correspond to their 16SrI subgroup assignments (Figure 2A). The two 16SrI-A strains, AYWB and NJAY, share 86.9% of their genomes and these shared regions are nearly identical (i.e., ANI = 99.7%). Because the NJAY assembly is highly fragmented and ∼10% smaller compared to the complete genome of AYWB, the low proportion of genomic segments shared may be explained by missing data in the NJAY assembly. The strain WBD shares only ∼93% ANI with all other 16SrI phytoplasmas compared and represents a distinct lineage. All of the remaining strains, including those assigned to 16SrI-B and TW1, share >97% ANI and could considered as a third cluster. This clustering result is consistent with the maximum likelihood phylogeny inferred using 303 single-copy genes shared by all strains (Figure 1B and Supplementary Table S2).
In the second approach, we considered only the putative coding regions (i.e., intact CDS and putative pseudogenes) of these genomes for calculating genomic similarities (Figure 2B). Compared to the genome-wide ANI approach, this gene-centric method is more laborious and could potentially be affected by annotation quality. Nevertheless, this approach has a higher confidence in inferring the homology prior to calculating sequence identities and provides a complementary method to the ANI approach. Although the exact genome similarities values differ slightly, the patterns found by using these two approaches are consistent.
Based on these results regarding genomic similarities and molecular phylogeny, it appeared that these 16SrI phytoplasmas could be classified into three species-level OTUs. These findings are consistent with a previous study that the conventional 97% identity threshold for delineating species in bacteria based on 16S rRNA gene is too low (Edgar, 2018). Similar to a previous study that examined diverse prokaryotes (Jain et al., 2018), the sequence identity values among these phytoplasmas also exhibited a bimodal distribution (Figure 3). Excluding the comparisons involving TW1, the within-OTU comparisons all have >97% genome-wide ANI, while between-OTU comparisons all have <94% ANI (Figures 2A, 3A). For sequence identities that were calculated based on only the genic regions, in which the homology could be inferred more confidently, a similar separation was found (i.e., within-OTU: >98.9%, between-OTU: <96%; see Figures 2B, 3B). These results suggest that those strains assigned to the same species-level OTU may have on-going homologous recombination that maintained high sequence similarities in their shared genomic regions. In contrast, certain genetic barriers exist to lower homologous recombination between different OTUs. In other words, these results provide further support to the hypothesis that the biological species concept, which is based on the barriers to homologous gene exchange, could be applied to bacteria as well (Bobay and Ochman, 2017). Alternatively, other hypotheses that may explain the discontinuity in genome similarities include ecological sweeps or stochastic neutral processes (Jain et al., 2018), and the exact biological mechanisms remained to be investigated.
Figure 3. Correlations between metrics of genome similarities. (A) Similarity scores based on whole genome comparisons. (B) Similarity scores based on genic regions. Filled symbols indicate the pairwise comparisons involving two complete genomes, empty symbols indicate the pairwise comparisons involving one or two draft genomes. The strain TW1 is a special case; its similarity scores in comparison with 16SrI-A and 16SrI-B strains are indicated by red and blue, respectively.
In addition to sequence identity, another important measurement for genetic similarity is the proportion of genomic regions shared. For this second measurement, a wide range was observed for both within- and between-OTU comparisons (Figures 2, 3). The wide spread was mostly explained by the inclusion of two highly fragmented draft genomes (i.e., OY-V and CYP; see Table 1), which resulted in underestimates that are artifacts stemming from missing data. However, because the non-redundant chromosomal regions that are more conserved have a higher probability of being included in a draft assembly, the inclusion of draft genomes may also result in overestimates of the genomic regions shared. Due to these uncertainties, we restricted the comparisons to those four strains with the complete genome sequences available (i.e., AYWB, OY-M, MBSP-M3, and De Villa). For this reduced data set with high confidence, strains assigned to the same species-level OTU share ∼78–93% of their chromosomal segments and ∼79–89% of their protein-coding genes, whereas strains assigned to different OTUs share ∼74–78% of their chromosomal segments and ∼70–75% of their protein-coding genes.
For an alternative approach of comparing the gene content divergence among these phytoplasmas, we conducted a principal coordinates analysis that considered gene copy numbers in addition to the patterns of gene presence or absence (Figure 4). With the exception of TW1, the results were consistent with the molecular phylogeny (Figure 1B) and pairwise genome similarities (Figures 2, 3). The first coordinate explained ∼40% of the variance and showed clear separation for the three species-level OTUs. The second coordinate explained ∼23% of the variance and showed that those 16SrI-B strains could be separated into three subgroups.
Figure 4. Principal coordinates analysis of gene content. Each genome is represented by one dot (color-coded according to the 16Sr subgroup assignment). The distance between dots corresponds to their differences in gene content. The % variance explained by each axis is provided in parentheses. The strains with the complete genome sequences available are highlighted in bold.
Taken together, these results supported the existence of multiple species-level OTUs within the previously described “Ca. P. asteris” (Lee et al., 2004) and provided quantitative guidelines for inferring the putative species boundaries based on overall genome similarities. Based on these findings, the critical questions to ask are the genotypic, phenotypic, and ecological differentiations among these OTUs. Unfortunately, given the limited number of strains with genome sequence available, it is not possible to obtain meaningful answers yet. For example, other 16SrI-A strains were known to infect different dicots or monocots, and were found in North America and Europe (Lee et al., 2004). However, no genomic information is available for these strains with different hosts and/or geographic distributions. It is important to characterize these strains and validate their OTU assignment before meaningful comparisons could be made between different OTUs. Similarly, more strains that belong to the same OTU with WBD must be sampled and characterized before the biological characteristics of this OTU could be described appropriately. These improvements in the understanding of genotypic and phenotypic differentiations among phytoplasmas are critical in informing future taxonomy revisions.
The Special Case of the TW1 Genome
Based on the aforementioned framework for phytoplasma classification, the strain TW1 represents a strange and difficult case. The high level of intra-genomic variation between its two copies of 16S rRNA genes has been noted in the initial characterization (Town et al., 2018). In our whole-genome analysis, this strain shares a higher proportion of its chromosomal segments and protein-coding genes with those 16SrI-A strains, yet has higher sequence identities with those 16SrI-B strains (Figures 2, 3). Because of these properties, this strain was clustered with those 16SrI-A strains based on gene content (Figure 4), while having closer phylogenetic relationships with those 16SrI-B strains (Figure 1B).
In the absence of an intuitive hypothesis that could explain the evolutionary processes leading to these conflicting findings, we examined the possibilities that these findings were results of artifacts. This TW1 genome was co-assembled using one long-read library (based on Oxford Nanopore MinION) and one short-read library (based on Illumina MiSeq) (Town et al., 2018). When the raw reads from these two sequencing libraries were examined separately, we found that the long-read library likely contains a 16SrI-A type phytoplasma. When the raw reads were mapped to the AYWB chromosome, all regions were covered and 3,725 sequence polymorphisms were found; when mapped to OY-M, 52,649 bp had no coverage and 31,401 sequence polymorphisms were found. In contrast, the short-read library likely contains a 16SrI-B type phytoplasma. When the raw reads were mapped to OY-M, 103,922 bp had no coverage and 9,489 sequence polymorphisms were found; when mapped to AYWB, 110,458 had no coverage and 29,323 sequence polymorphisms were found.
Based on these findings, a possible explanation is that during the assembly process, the long reads produced several 16SrI-A type scaffolds, which resulted in the high similarities of the TW1 genome content to 16SrI-A as observed (Figures 2–4). Subsequently, during the polishing stage, the exact sequences of these scaffolds were modified based on the short-read library, thus explaining the high sequence similarity of this TW1 genome to those 16SrI-B strains. In other words, in its current form, this TW1 genome assembly is likely an artifact that was incorrectly assembled by combining two sequencing libraries with phytoplasmas belonging to different OTUs, rather than a true representative of an existing phytoplasma strain. Due to the concern that the inclusion of this genome may have introduced biases, we repeated all phylogenetic inferences and excluded the TW1 sequences prior to multiple sequence alignment. The resulting tree topologies were consistent in terms of OTU assignments; the minor differences were related to the placement of CYP (Figure 1 and Supplementary Figure S1).
Development of Molecular Markers for Multilocus Sequence Analysis
To develop molecular markers that may be useful for the MLSA of these phytoplasmas and their relatives, we examined the multiple sequence alignments of those 303 single-copy genes shared by all of the genomes analyzed. The selection criteria were based on: (1) the number of informative sites that may be used to distinguish the three species-level OTUs identified, (2) the possibility of designing PCR primers for a product in the size range of ∼700–800 bp and covers a large number of informative sites, and (3) the chromosomal locations. The first two criteria were aimed to make the Sanger sequencing of these markers as cost-effective as possible, while the third criterion was aimed to make the MLSA robust against recombination.
Based on these criteria, we selected five markers to supplement the 16S rRNA gene. These five makers all have high densities of informative sites, particularly when compared to those markers that have been developed previously (Table 2 and Supplementary Figure S1). The substitute rates of these marker genes are relatively high but do not deviate strongly when compared to other shared single-copy genes (Supplementary Figure S2).
Based on the functional annotation, two of these selected markers (i.e., DegV family protein and TIGR00282 family metallophosphoesterase) are not generally considered as house-keeping genes. These results reflected our philosophy of marker selection (i.e., prioritize the candidates with the highest resolving power for the target taxa regardless of other irrelevant attributes). Similar to our criteria and results, a pioneer study on genome-enable MLSA marker design in phytoplasma selected several genes that were annotated as hypothetical proteins (Kakizawa and Kamagata, 2014). In terms of phylogenetic distribution, we found that four of these five selected genes are present as single-copy genes in other phytoplasmas, including “Ca. P. australiense” of the 16SrXII group (Tran-Nguyen et al., 2008), “Ca. P. aurantifolia” of the 16SrII group (Chung et al., 2013), “Ca. P. pruni” of the 16SrIII group (Lee et al., 2015), “Ca. P. ziziphi’ of the 16SrV group (Wang et al., 2018a), “Ca. P. oryzae” of the 16SrXI group (GenBank accession NZ_JHUK00000000), and “Ca. P. mali” of the 16SrX group (Kube et al., 2008). The only exception is degV, which is absent in two of the draft genomes (i.e., ‘Ca. P. aurantifolia’ and ‘Ca. P. oryzae’; may be located in the unassembled regions) and has two tandem copies in the complete genome sequence of ‘Ca. P. mali’. Nevertheless, the sequence divergence levels are too high for all these genes, such that the primers designed in this study (Table 2) may not work for these more divergent phytoplasmas outside of the 16SrI group. These results are consistent with our expectations that different MLSA markers are necessary for different phylogenetic depths.
In terms of chromosomal locations, these markers are well separated in all genomes examined despite the extensive genome rearrangements observed (Figure 5). Although the extent of recombination in phytoplasmas is not well understood, recent studies on Xylella fastidiosa (i.e., another insect-transmitted plant pathogen that is restricted to vascular tissues) revealed that recombination affected ∼0.5–15% of the genome, depending on the exact subspecies (Potnis et al., 2019; Vanhove et al., 2019). The recombined fragments could have sizes up to 31 kb, with an average of 1 kb. Based on these estimates, it is unlikely that one single recombination event would affect more than one of these five markers. Thus, the results of MLSA should be robust against the interference of recombination.
Figure 5. Genomic locations of marker genes. Nucleotide sequence identity of the conserved regions between genomes are illustrated (red: match on the same strand; blue: match on the opposite strand).
When we used these five newly developed markers for molecular phylogenetic inference, the concatenated alignment (together with the 16S rRNA gene, which is expected to be sequenced in the initial characterization of newly collected strains) produced a phylogeny (Figure 1C and Supplementary Figure S1C) that is comparable to the result from genome-scale analysis (Figure 1B and Supplementary Figure S1B). This result demonstrated the usefulness of these markers for future genotyping work. With better sampling of phytoplasma strains and more accurate classification of those strains using these markers, informed decision could be made to select representatives for genome sequencing efforts, which could contribute to the study of phytoplasma genetic diversity and evolutionary history, as well as improve the taxonomy. Moreover, as the sampling focus changes in the future (e.g., a higher resolution is desired for one of the OTUs), the list of conserved single-copy genes and the substitution rates of those genes (Supplementary Table S2) could provide a guide for developing more suitable MLSA markers.
Phylogenetic Distribution of Phytoplasma Effector Genes
With the well-resolved phylogeny (Figure 1B) as a framework, we investigated the distribution of putative effector genes among these phytoplasmas. The results indicated a high level of variation in the gene counts (Figure 6 and Supplementary Table S3). One clade (i.e., MBSP-M3, De Villa, and LD1) have only ∼10–14 genes that encode putative secreted proteins. These numbers are much lower than the ∼36–45 found in their sister clade (i.e., OY-M, OY-V, and DY2014). The low numbers are not artifacts of incomplete genome sequences. In fact, the two strains with the lowest numbers (i.e., MBSP-M3 and De Villa) both have the complete chromosomal sequences available. These findings suggest that there may have been lineage-specific reduction or expansion in effector gene counts, although the underlying ecological factors are unclear.
Figure 6. Phylogenetic distribution of putative effector genes. The cladogram that illustrates the phylogenetic relationships among these strains is based on Figure 1B (i.e., maximum likelihood phylogeny based on 303 shared single-copy genes). The number in parentheses after the strain name indicates the number of genes that encode putative secreted proteins in each genome (including putative pseudogenes). The strains with the complete genome sequences available are highlighted in bold. On the right-hand side, the table lists the locus tags of homologs for those four phytoplasma effectors that have been characterized. Those locus tags with an asterisk (*) were annotated as putative pseudogenes.
For those 16SrI phytoplasma effectors that have been characterized experimentally, namely SAP05 (Gamboa et al., 2019; Huang and Hogenhout, 2019), SAP11/SWP1 (Bai et al., 2009; Sugio et al., 2011a; Lu et al., 2014; Chang et al., 2018; Wang et al., 2018b, c; Lu et al., 2014), SAP54/PHYL1 (MacLean et al., 2011; Maejima et al., 2014; Orlovskis and Hogenhout, 2016), and TENGU (Hoshi et al., 2009; Sugawara et al., 2013; Minato et al., 2014), only TENGU is conserved among all strains. Although some cases of gene absence may be artifacts of incomplete genome assemblies, the strain De Villa (which has the complete chromosomal sequence available) lacks all other three effector genes, indicating that none of those three effectors is essential.
The variable pattern of effector gene presence and absence may be linked to the molecular evolution of potential mobile units (PMU) in these phytoplasma genomes (Bai et al., 2006). These mobile genetic elements were known to facilitate horizontal gene transfers between either closely or distantly related strains and may be gained or lost rapidly (Chung et al., 2013; Ku et al., 2013; Orlovskis et al., 2017; Cho et al., 2019b; Seruga Music et al., 2019). When we examined the chromosomal locations of those four effector genes, TENGU is not associated with PMUs, while the other three are. Finally, as noted in our previous study (Cho et al., 2019b), DY2014 and OY-V both have two copies of SAP11 homologs. One copy is located within a 16SrI-B type PMU, suggesting vertical inheritance, while the other copy is located within a 16SrI-A type PMU, suggesting horizontal acquisition. These observations provided further support to the importance of PMUs in phytoplasma effector evolution.
Conclusion
By quantifying the genetic divergence among a group of closely-related phytoplasmas, this work demonstrated a discontinuity in their similarities of genomic content and homologous sequences. The results suggested that these phytoplasmas could be classified into multiple distinct taxonomic units equivalent to species in other bacteria (Jain et al., 2018). Importantly, the widely used classification system for phytoplasmas that is based on the sequencing and RFLP analysis of their 16S rRNA genes does not provide sufficient resolution or accuracy. Rather, it is important to incorporate genomic information in the future revisions of phytoplasma taxonomy. A previous proposal of using 95% ANI as a cutoff to delineate bacterial species (Jain et al., 2018) appears to work well and warrants further consideration to be adopted. In recognition of the difficulties involved in genomic studies of uncultivated bacteria, this work also demonstrated the feasibility of using genome analysis to develop cost-effective MLSA markers, which will be useful for practical purposes. For future directions, it is critical to improve the taxon sampling of phytoplasma genomes strategically, such that similar approaches may be expand to genus-wide analysis to facilitate the studies of these important plant pathogens. Moreover, continuing effort in establishing axenic culture is critical for the study of these important plant pathogens. In the absence of axenic culture, it is important that the reference strains, particularly those with provisional species status, are made available to the scientific community, preferably in the form of a centralized micropropagation collection as suggested by the international phytoplasma working team (The IRPCM Phytoplasma/Spiroplasma Working Team, 2004). If live culture within plant or insect hosts is not possible, then the DNA samples should be made available. For genome sequencing projects, the raw sequencing results should be made publicly available together with the assembled genome sequences at the time of publication, such that the assemblies could be validated independently. The TW1 genome (Town et al., 2018) provided a good example to illustrate the importance of this point. Although the publication of an erroneous genome sequence is unfortunate, the raw reads made available by those authors allowed us to investigate and identify the issues; such openness is appreciated and important for science to move forward.
Author’s Note
This manuscript was posted as a preprint at bioRxiv (10.1101/2020.01.31.928135).
Data Availability Statement
All datasets generated for this study are included in the article/Supplementary Material.
Author Contributions
S-TC, H-JK, and C-HK performed the experiments and analyzed the data. WH and SH provided intellectual input. S-TC prepared the figures and Supplementary Materials. S-TC and C-HK wrote the manuscript. SH edited the manuscript. C-HK designed the experiments, acquired the funding, and supervised the project. All authors contributed to the article and approved the submitted version.
Funding
The funding was provided by Academia Sinica and the Ministry of Science and Technology of Taiwan (MOST 106-2311-B-001-028-MY3) to C-HK. SH and WH were supported by the Human Frontier Science Program grant RGP0024/2105 with additional support from Biotechnology and Biological Sciences Research Council grant BB/P012574/1 and the John Innes Foundation. The funders have no role in study design, data collection and interpretation, or the decision to submit this work for publication.
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Acknowledgments
We thank Ms. Mei-Jane Fang in the Genomic Technology Core Facility (Institute of Plant and Microbial Biology, Academia Sinica) for providing Sanger sequencing service. Comments received during the peer review process greatly improved this work. This manuscript is not meant to be a comprehensive review of all the available literature on phytoplasma classification; as such we apologize if we did not cite all publications on this topic in the manuscript.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmicb.2020.01531/full#supplementary-material
FIGURE S1 | Maximum likelihood phylogenies. Based on Figure 1, the TW1 sequences were excluded prior to multiple sequence alignment. The numbers on the internal branches indicate the level of bootstrap support based on 1,000 resampling; only values ≥60% are shown. (A) 16S rRNA genes. The strain name, locus tag of the 16S rRNA gene (in parentheses), and the 16SrI subgroup classification (in square brackets) are labeled. A question mark ‘?’ in the subgroup classification indicates that the sequence is classified as a new subgroup, the existing subgroup with the highest similarity is provided. (B) The 303 single-copy coding genes shared by all strains. (C) 16S rRNA gene and the five markers developed in this study; see Table 2 for detailed information.
FIGURE S2 | Multiple sequence alignments of selected marker genes. The primer sites are highlighted in orange. Shades of blue colors in the alignment indicate the levels of sequence conservation.
FIGURE S3 | Distribution of substitution rates among those single-copy genes shared by the 11 phytoplasma genomes analyzed. The synonymous substitution rate (Ks) and non-synonymous substitution rate (Ka) were calculated based on the homologs found in AYWB and MBSP-M3. Among those 303 single-copy genes, nine were annotated as pseudogenes and were omitted from the substitution rate calculation.
TABLE S1 | Pairwise nucleotide sequence identity of the molecular markers used in the phylogenetic inference. Numbers in the diagonal cells indicate the length of each unaligned sequence.
TABLE S2 | List of the 303 single-copy genes shared by all of the 11 phytoplasma genomes analyzed in this study. The gene identity, genomic location, and functional annotation of these genes are provided based on the information from AYWB. The substitution rate estimates are based on the comparison between AYWB and MBSP-M3.
TABLE S3 | List of the putative effector genes. Each gene is uniquely identified by its locus_tag. The “Cluster_id” field provides the assignment of homologous gene cluster; genes sharing the same Cluster_id are considered as homologs. Note that SAP11/SWP1 homologs are separated into two clusters (i.e., 503 and 598) due to high levels of sequence divergence. Additionally, one SAP05 homolog is in a separate cluster (i.e., 1117) due to truncation. The gene names and product description are based on existing GenBank annotation; additional information (e.g., homology information with known effector or pseudogene status) are provided in the “Note” field.
References
Armenteros, J. J. A., Tsirigos, K. D., Sønderby, C. K., Petersen, T. N., Winther, O., Brunak, S., et al. (2019). SignalP 5.0 improves signal peptide predictions using deep neural networks. Nat. Biotechnol. 37, 420–423. doi: 10.1038/s41587-019-0036-z
Arocha, Y., Antesana, O., Montellano, E., Franco, P., Plata, G., and Jones, P. (2007). ‘Candidatus Phytoplasma lycopersici’, a phytoplasma associated with ‘hoja de perejil’ disease in Bolivia. Int. J. Syst. Evol. Microbiol. 57, 1704–1710. doi: 10.1099/ijs.0.64851-0
Bai, X., Correa, V. R., Toruño, T. Y., Ammar, E.-D., Kamoun, S., and Hogenhout, S. A. (2009). AY-WB phytoplasma secretes a protein that targets plant cell nuclei. Mol. Plant Microbe Interact. 22, 18–30. doi: 10.1094/MPMI-22-1-0018
Bai, X., Zhang, J., Ewing, A., Miller, S. A., Jancso Radek, A., Shevchenko, D. V., et al. (2006). Living with genome instability: the adaptation of phytoplasmas to diverse environments of their insect and plant hosts. J. Bacteriol. 188, 3682–3696. doi: 10.1128/JB.188.10.3682-3696.2006
Benson, D. A., Cavanaugh, M., Clark, K., Karsch-Mizrachi, I., Ostell, J., Pruitt, K. D., et al. (2018). GenBank. Nucleic Acids Res. 46, D41–D47. doi: 10.1093/nar/gkx1094
Bertaccini, A., and Lee, I.-M. (2018). “Phytoplasmas: an update,” in Phytoplasmas: Plant Pathogenic Bacteria - I: Characterisation and Epidemiology of Phytoplasma – Associated Diseases, eds G. P. Rao, A. Bertaccini, N. Fiore, and L. W. Liefting (Singapore: Springer), 1–29. doi: 10.1007/978-981-13-0119-3_1
Bobay, L.-M., and Ochman, H. (2017). Biological species are universal across life’s domains. Genome Biol. Evol. 9, 491–501. doi: 10.1093/gbe/evx026
Brown, D. R. (2010). “Phylum XVI. Tenericutes Murray 1984a, 356VP (Effective publication: Murray 1984b, 33.),” in Bergey’s Manual of Systematic Bacteriology, eds N. R. Krieg, J. T. Staley, D. R. Brown, B. P. Hedlund, B. J. Paster, N. L. Ward, et al. (New York: Springer), 567–723. doi: 10.1007/978-0-387-68572-4_5
Camacho, C., Coulouris, G., Avagyan, V., Ma, N., Papadopoulos, J., Bealer, K., et al. (2009). BLAST+: architecture and applications. BMC Bioinform. 10:421. doi: 10.1186/1471-2105-10-421
Chang, S. H., Tan, C. M., Wu, C.-T., Lin, T.-H., Jiang, S.-Y., Liu, R.-C., et al. (2018). Alterations of plant architecture and phase transition by the phytoplasma virulence factor SAP11. J. Exp. Bot. 69, 5389–5401. doi: 10.1093/jxb/ery318
Chen, L.-L., Chung, W.-C., Lin, C.-P., and Kuo, C.-H. (2012). Comparative analysis of gene content evolution in phytoplasmas and mycoplasmas. PLoS One 7:e34407. doi: 10.1371/journal.pone.0034407
Chen, W., Li, Y., Wang, Q., Wang, N., and Wu, Y. (2014). Comparative genome analysis of wheat blue dwarf phytoplasma, an obligate pathogen that causes wheat blue dwarf disease in China. PLoS One 9:e96436. doi: 10.1371/journal.pone.0096436
Cho, S.-T., Kung, H.-J., and Kuo, C.-H. (2019a). Genomic similarities among 16SrI phytoplasmas and implications on species boundaries. Phytopathog. Mollic. 9:73. doi: 10.5958/2249-4677.2019.00037.9
Cho, S.-T., Lin, C.-P., and Kuo, C.-H. (2019b). Genomic characterization of the periwinkle leaf yellowing (PLY) phytoplasmas in Taiwan. Front. Microbiol. 10:2194. doi: 10.3389/fmicb.2019.02194
Chung, W.-C., Chen, L.-L., Lo, W.-S., Lin, C.-P., and Kuo, C.-H. (2013). Comparative analysis of the peanut witches’-broom phytoplasma genome reveals horizontal transfer of potential mobile units and effectors. PLoS One 8:e62770. doi: 10.1371/journal.pone.0062770
Coetzee, B., Douglas-Smit, N., Maree, H. J., Burger, J. T., Krüger, K., and Pietersen, G. (2019). Draft genome sequence of a “Candidatus Phytoplasma asteris”-related strain (aster yellows, subgroup 16SrI-B) from South Africa. Microbiol. Resour. Announc. 8:e00148-19. doi: 10.1128/MRA.00148-19
Dumonceaux, T. J., Green, M., Hammond, C., Perez, E., and Olivier, C. (2014). Molecular diagnostic tools for detection and differentiation of phytoplasmas based on chaperonin-60 reveal differences in host plant infection patterns. PLoS One 9:e116039. doi: 10.1371/journal.pone.0116039
Edgar, R. C. (2004). MUSCLE: multiple sequence alignment with high accuracy and high throughput. Nucl. Acids Res. 32, 1792–1797. doi: 10.1093/nar/gkh340
Edgar, R. C. (2018). Updating the 97% identity threshold for 16S ribosomal RNA OTUs. Bioinformatics 34, 2371–2375. doi: 10.1093/bioinformatics/bty113
Firrao, G., Martini, M., Ermacora, P., Loi, N., Torelli, E., Foissac, X., et al. (2013). Genome wide sequence analysis grants unbiased definition of species boundaries in “Candidatus Phytoplasma.”. Syst. Appl. Microbiol. 36, 539–548. doi: 10.1016/j.syapm.2013.07.003
Gamboa, C., Cui, W., Quiroga, N., Fernández, C., Fiore, N., and Zamorano, A. (2019). Identification of 16SrIII-J phytoplasma effectors using a viral vector. Phytopathog. Mollic. 9:229. doi: 10.5958/2249-4677.2019.00115.4
Glaeser, S. P., and Kämpfer, P. (2015). Multilocus sequence analysis (MLSA) in prokaryotic taxonomy. Syst. Appl. Microbiol. 38, 237–245. doi: 10.1016/j.syapm.2015.03.007
Guindon, S., and Gascuel, O. (2003). A simple, fast, and accurate algorithm to estimate large phylogenies by maximum likelihood. Syst. Biol. 52, 696–704. doi: 10.1080/10635150390235520
Gundersen, D. E., and Lee, I.-M. (1996). Ultrasensitive detection of phytoplasmas by nested-PCR assays using two universal primer pairs. Phytopathol. Mediterr. 35, 144–151.
Guy, L., Roat Kultima, J., and Andersson, S. G. E. (2010). genoPlotR: comparative gene and genome visualization in R. Bioinformatics 26, 2334–2335. doi: 10.1093/bioinformatics/btq413
Hogenhout, S. A., Oshima, K., Ammar, E.-D., Kakizawa, S., Kingdom, H. N., and Namba, S. (2008). Phytoplasmas: bacteria that manipulate plants and insects. Mol. Plant Pathol. 9, 403–423. doi: 10.1111/j.1364-3703.2008.00472.x
Hoshi, A., Oshima, K., Kakizawa, S., Ishii, Y., Ozeki, J., Hashimoto, M., et al. (2009). A unique virulence factor for proliferation and dwarfism in plants identified from a phytopathogenic bacterium. Proc. Natl. Acad. Sci. U.S.A. 106, 6416–6421. doi: 10.1073/pnas.0813038106
Huang, W., and Hogenhout, S. A. (2019). “Phytoplasma effectors have converged onto degrading plant transcription factors with fundamental roles in plant development and defense to insects,” in IS-MPMI XVIII Congress Abstracts of Poster Presentations, (Glasgow: IS-MPMI), 158. doi: 10.1094/MPMI-32-10-S1.1
Jain, C., Rodriguez-R, L. M., Phillippy, A. M., Konstantinidis, K. T., and Aluru, S. (2018). High throughput ANI analysis of 90K prokaryotic genomes reveals clear species boundaries. Nat. Commun. 9:5114. doi: 10.1038/s41467-018-07641-9
Kakizawa, S., and Kamagata, Y. (2014). A multiplex-PCR method for strain identification and detailed phylogenetic analysis of AY-group phytoplasmas. Plant Dis. 98, 299–305. doi: 10.1094/PDIS-03-13-0216-RE
Kakizawa, S., Makino, A., Ishii, Y., Tamaki, H., and Kamagata, Y. (2014). Draft genome sequence of “Candidatus Phytoplasma asteris” strain OY-V, an unculturable plant-pathogenic bacterium. Genome Announc. 2:e00944-14. doi: 10.1128/genomeA.00944-14
Konstantinidis, K. T., Rosselló-Móra, R., and Amann, R. (2017). Uncultivated microbes in need of their own taxonomy. ISME J. 11, 2399–2406. doi: 10.1038/ismej.2017.113
Konstantinidis, K. T., and Tiedje, J. M. (2005). Genomic insights that advance the species definition for prokaryotes. Proc. Natl. Acad. Sci. U.S.A. 102, 2567–2572. doi: 10.1073/pnas.0409727102
Krogh, A., Larsson, B., von Heijne, G., and Sonnhammer, E. L. L. (2001). Predicting transmembrane protein topology with a hidden Markov model: application to complete genomes. J. Mol. Biol. 305, 567–580. doi: 10.1006/jmbi.2000.4315
Ku, C., Lo, W.-S., and Kuo, C.-H. (2013). Horizontal transfer of potential mobile units in phytoplasmas. Mob. Genet. Elements 3:e26145. doi: 10.4161/mge.26145
Kube, M., Schneider, B., Kuhl, H., Dandekar, T., Heitmann, K., Migdoll, A., et al. (2008). The linear chromosome of the plant-pathogenic mycoplasma ‘Candidatus Phytoplasma mali.’ BMC Genomics 9:306. doi: 10.1186/1471-2164-9-306
Kuo, C.-H., Moran, N. A., and Ochman, H. (2009). The consequences of genetic drift for bacterial genome complexity. Genome Res 19, 1450–1454. doi: 10.1101/gr.091785.109
Kuo, C.-H., and Ochman, H. (2009). Inferring clocks when lacking rocks: the variable rates of molecular evolution in bacteria. Biol. Direct 4:35. doi: 10.1186/1745-6150-4-35
Lee, I.-M., Davis, R. E., and Gundersen-Rindal, D. E. (2000). Phytoplasma: phytopathogenic mollicutes. Annu. Rev. Microbiol. 54, 221–255. doi: 10.1146/annurev.micro.54.1.221
Lee, I.-M., Gundersen-Rindal, D. E., Davis, R. E., and Bartoszyk, I. M. (1998). Revised classification scheme of phytoplasmas based on RFLP analyses of 16S rRNA and ribosomal protein gene sequences. Int. J. Syst. Evol. Microbiol. 48, 1153–1169. doi: 10.1099/00207713-48-4-1153
Lee, I.-M., Gundersen-Rindal, D. E., Davis, R. E., Bottner, K. D., Marcone, C., and Seemüller, E. (2004). ‘Candidatus Phytoplasma asteris’, a novel phytoplasma taxon associated with aster yellows and related diseases. Int. J. Syst. Evol. Microbiol. 54, 1037–1048. doi: 10.1099/ijs.0.02843-0
Lee, I.-M., Hammond, R., Davis, R., and Gundersen, D. (1993). Universal amplification and analysis of pathogen 16S rDNA for classification and identification of mycoplasmalike organisms. Phytopathology 83, 834–842. doi: 10.1094/phyto-83-834
Lee, I.-M., Shao, J., Bottner-Parker, K. D., Gundersen-Rindal, D. E., Zhao, Y., and Davis, R. E. (2015). Draft genome sequence of “Candidatus Phytoplasma pruni” strain CX, a plant-pathogenic bacterium. Genome Announc. 3:e01117-15. doi: 10.1128/genomeA.01117-15
Lee, I.-M., Zhao, Y., and Bottner, K. D. (2006). SecY gene sequence analysis for finer differentiation of diverse strains in the aster yellows phytoplasma group. Mol. Cell. Probes 20, 87–91. doi: 10.1016/j.mcp.2005.10.001
Li, H. (2018). Minimap2: pairwise alignment for nucleotide sequences. Bioinformatics 34, 3094–3100. doi: 10.1093/bioinformatics/bty191
Li, H., and Durbin, R. (2009). Fast and accurate short read alignment with Burrows–Wheeler transform. Bioinformatics 25, 1754–1760. doi: 10.1093/bioinformatics/btp324
Li, H., Handsaker, B., Wysoker, A., Fennell, T., Ruan, J., Homer, N., et al. (2009). The sequence alignment/map format and SAMtools. Bioinformatics 25, 2078–2079. doi: 10.1093/bioinformatics/btp352
Li, L., Stoeckert, C. J., and Roos, D. S. (2003). OrthoMCL: identification of ortholog groups for eukaryotic genomes. Genome Res. 13, 2178–2189. doi: 10.1101/gr.1224503
Liefting, L. W., Andersen, M. T., Beever, R. E., Gardner, R. C., and Forster, R. L. (1996). Sequence heterogeneity in the two 16S rRNA genes of Phormium yellow leaf phytoplasma. Appl. Environ. Microbiol 62, 3133–3139. doi: 10.1128/aem.62.9.3133-3139.1996
Lim, P. O., and Sears, B. B. (1992). Evolutionary relationships of a plant-pathogenic mycoplasmalike organism and Acholeplasma laidlawii deduced from two ribosomal protein gene sequences. J. Bacteriol. 174, 2606–2611. doi: 10.1128/jb.174.8.2606-2611.1992
Lo, W.-S., Chen, L.-L., Chung, W.-C., Gasparich, G. E., and Kuo, C.-H. (2013). Comparative genome analysis of Spiroplasma melliferum IPMB4A, a honeybee-associated bacterium. BMC Genomics 14:22. doi: 10.1186/1471-2164-14-22
Lo, W.-S., Gasparich, G. E., and Kuo, C.-H. (2018). Convergent evolution among ruminant-pathogenic Mycoplasma involved extensive gene content changes. Genome Biol. Evol. 10, 2130–2139. doi: 10.1093/gbe/evy172
Lu, Y.-T., Li, M.-Y., Cheng, K.-T., Tan, C. M., Su, L.-W., Lin, W.-Y., et al. (2014). Transgenic plants that express the phytoplasma effector SAP11 show altered phosphate starvation and defense responses. Plant Physiol. 164, 1456–1469. doi: 10.1104/pp.113.229740
MacLean, A. M., Sugio, A., Makarova, O. V., Findlay, K. C., Grieve, V. M., Tóth, R., et al. (2011). Phytoplasma effector SAP54 induces indeterminate leaf-like flower development in Arabidopsis plants. Plant Physiol. 157, 831–841. doi: 10.1104/pp.111.181586
Maejima, K., Iwai, R., Himeno, M., Komatsu, K., Kitazawa, Y., Fujita, N., et al. (2014). Recognition of floral homeotic MADS domain transcription factors by a phytoplasmal effector, phyllogen, induces phyllody. Plant J. 78, 541–554. doi: 10.1111/tpj.12495
Martini, M., Lee, I.-M., Bottner, K. D., Zhao, Y., Botti, S., Bertaccini, A., et al. (2007). Ribosomal protein gene-based phylogeny for finer differentiation and classification of phytoplasmas. Int. J. Syst. Evol. Microbiol. 57, 2037–2051. doi: 10.1099/ijs.0.65013-0
Martini, M., Quaglino, F., and Bertaccini, A. (2019). “Multilocus genetic characterization of phytoplasmas,” in Phytoplasmas: Plant Pathogenic Bacteria - III: Genomics, Host Pathogen Interactions and Diagnosis, eds A. Bertaccini, K. Oshima, M. Kube, and G. P. Rao (Singapore: Springer), 161–200. doi: 10.1007/978-981-13-9632-8_9
Minato, N., Himeno, M., Hoshi, A., Maejima, K., Komatsu, K., Takebayashi, Y., et al. (2014). The phytoplasmal virulence factor TENGU causes plant sterility by downregulating of the jasmonic acid and auxin pathways. Sci. Rep. 4:7399. doi: 10.1038/srep07399
Mitrović, J., Kakizawa, S., Duduk, B., Oshima, K., Namba, S., and Bertaccini, A. (2011). The groEL gene as an additional marker for finer differentiation of ‘Candidatus Phytoplasma asteris’-related strains. Ann. Appl. Biol. 159, 41–48. doi: 10.1111/j.1744-7348.2011.00472.x
Orlovskis, Z., Canale, M. C., Haryono, M., Lopes, J. R. S., Kuo, C.-H., and Hogenhout, S. A. (2017). A few sequence polymorphisms among isolates of maize bushy stunt phytoplasma associate with organ proliferation symptoms of infected maize plants. Ann. Bot. 119, 869–884. doi: 10.1093/aob/mcw213
Orlovskis, Z., and Hogenhout, S. A. (2016). A bacterial parasite effector mediates insect vector attraction in host plants independently of developmental changes. Front. Plant Sci. 7:885. doi: 10.3389/fpls.2016.00885
Oshima, K., Kakizawa, S., Nishigawa, H., Jung, H.-Y., Wei, W., Suzuki, S., et al. (2004). Reductive evolution suggested from the complete genome sequence of a plant-pathogenic phytoplasma. Nat. Genet. 36, 27–29. doi: 10.1038/ng1277
Pacifico, D., Galetto, L., Rashidi, M., Abbà, S., Palmano, S., Firrao, G., et al. (2015). Decreasing global transcript levels over time suggest that phytoplasma cells enter stationary phase during plant and insect colonization. Appl. Environ. Microbiol. 81, 2591–2602. doi: 10.1128/AEM.03096-14
Parks, D. H., Chuvochina, M., Waite, D. W., Rinke, C., Skarshewski, A., Chaumeil, P.-A., et al. (2018). A standardized bacterial taxonomy based on genome phylogeny substantially revises the tree of life. Nat. Biotechnol. 36, 996–1004. doi: 10.1038/nbt.4229
Popescu, A.-A., Huber, K. T., and Paradis, E. (2012). ape 3.0: New tools for distance-based phylogenetics and evolutionary analysis in R. Bioinformatics 28, 1536–1537. doi: 10.1093/bioinformatics/bts184
Potnis, N., Kandel, P. P., Merfa, M. V., Retchless, A. C., Parker, J. K., Stenger, D. C., et al. (2019). Patterns of inter- and intrasubspecific homologous recombination inform eco-evolutionary dynamics of Xylella fastidiosa. ISME J. 13, 2319–2333. doi: 10.1038/s41396-019-0423-y
R Core Team (2019). R: A Language and Environment for Statistical Computing. Vienna: R Foundation for Statistical Computing.
Sawayanagi, T., Horikoshi, N., Kanehira, T., Shinohara, M., Bertaccini, A., Cousin, M.-T., et al. (1999). ‘Candidatus Phytoplasma japonicum’, a new phytoplasma taxon associated with Japanese Hydrangea phyllody. Int. J. Syst. Evol. Microbiol. 49, 1275–1285. doi: 10.1099/00207713-49-3-1275
Schloss, P. D., and Handelsman, J. (2005). Introducing DOTUR, a computer program for defining operational taxonomic units and estimating species richness. Appl. Environ. Microbiol. 71, 1501–1506. doi: 10.1128/AEM.71.3.1501-1506.2005
Schneider, B., and Gibb, K. S. (1997). Sequence and RFLP analysis of the elongation factor Tu gene used in differentiation and classification of phytoplasmas. Microbiology 143, 3381–3389. doi: 10.1099/00221287-143-10-3381
Seruga Music, M., Samarzija, I., Hogenhout, S. A., Haryono, M., Cho, S.-T., and Kuo, C.-H. (2019). The genome of ‘Candidatus Phytoplasma solani’ strain SA-1 is highly dynamic and prone to adopting foreign sequences. Syst. Appl. Microbiol. 42, 117–127. doi: 10.1016/j.syapm.2018.10.008
Sparks, M. E., Bottner-Parker, K. D., Gundersen-Rindal, D. E., and Lee, I.-M. (2018). Draft genome sequence of the New Jersey aster yellows strain of ‘Candidatus Phytoplasma asteris.’. PLoS One 13:e0192379. doi: 10.1371/journal.pone.0192379
Sugawara, K., Honma, Y., Komatsu, K., Himeno, M., Oshima, K., and Namba, S. (2013). The alteration of plant morphology by small peptides released from the proteolytic processing of the bacterial peptide TENGU. Plant Physiol. 162, 2005–2014. doi: 10.1104/pp.113.218586
Sugio, A., Kingdom, H. N., MacLean, A. M., Grieve, V. M., and Hogenhout, S. A. (2011a). Phytoplasma protein effector SAP11 enhances insect vector reproduction by manipulating plant development and defense hormone biosynthesis. Proc. Natl. Acad. Sci. U.S.A. 108, E1254–E1263. doi: 10.1073/pnas.1105664108
Sugio, A., MacLean, A. M., Kingdom, H. N., Grieve, V. M., Manimekalai, R., and Hogenhout, S. A. (2011b). Diverse targets of phytoplasma effectors: from plant development to defense against insects. Annu. Rev. Phytopathol. 49, 175–195. doi: 10.1146/annurev-phyto-072910-095323
Suyama, M., Torrents, D., and Bork, P. (2006). PAL2NAL: robust conversion of protein sequence alignments into the corresponding codon alignments. Nucl. Acids Res. 34, W609–W612. doi: 10.1093/nar/gkl315
The IRPCM Phytoplasma/Spiroplasma Working Team (2004). ‘Candidatus Phytoplasma’, A TAXON FOR THE WALL-LESS, NON-helical prokaryotes that colonize plant phloem and insects. Int. J. Syst. Evol. Microbiol. 54, 1243–1255. doi: 10.1099/ijs.0.02854-0
Tomkins, M., Kliot, A., Marée, A. F., and Hogenhout, S. A. (2018). A multi-layered mechanistic modelling approach to understand how effector genes extend beyond phytoplasma to modulate plant hosts, insect vectors and the environment. Curr. Opin. Plant Biol. 44, 39–48. doi: 10.1016/j.pbi.2018.02.002
Town, J. R., Wist, T., Perez-Lopez, E., Olivier, C. Y., and Dumonceaux, T. J. (2018). Genome sequence of a plant-pathogenic bacterium, “Candidatus Phytoplasma asteris” strain TW1. Microbiol. Resour. Announc. 7:e01109-18. doi: 10.1128/MRA.01109-18
Tran-Nguyen, L. T. T., Kube, M., Schneider, B., Reinhardt, R., and Gibb, K. S. (2008). Comparative genome analysis of “Candidatus Phytoplasma australiense” (subgroup tuf-Australia I; rp-A) and “Ca, Phytoplasma asteris” strains OY-M and AY-WB. J. Bacteriol. 190, 3979–3991. doi: 10.1128/JB.01301-07
Vanhove, M., Retchless, A. C., Sicard, A., Rieux, A., Coletta-Filho, H. D., Fuente, L. D. L., et al. (2019). Genomic diversity and recombination among Xylella fastidiosa subspecies. Appl. Environ. Microbiol. 85:e02972-18. doi: 10.1128/AEM.02972-18
Wang, J., Song, L., Jiao, Q., Yang, S., Gao, R., Lu, X., et al. (2018a). Comparative genome analysis of jujube witches’-broom phytoplasma, an obligate pathogen that causes jujube witches’-broom disease. BMC Genomics 19:689. doi: 10.1186/s12864-018-5075-1
Wang, N., Li, Y., Chen, W., Yang, H. Z., Zhang, P. H., and Wu, Y. F. (2018b). Identification of wheat blue dwarf phytoplasma effectors targeting plant proliferation and defence responses. Plant Pathol. 67, 603–609. doi: 10.1111/ppa.12786
Wang, N., Yang, H., Yin, Z., Liu, W., Sun, L., and Wu, Y. (2018c). Phytoplasma effector SWP1 induces witches’ broom symptom by destabilizing the TCP transcription factor BRANCHED1. Mol. Plant Pathol. 19, 2623–2634. doi: 10.1111/mpp.12733
Waterhouse, A. M., Procter, J. B., Martin, D. M. A., Clamp, M., and Barton, G. J. (2009). Jalview Version 2—a multiple sequence alignment editor and analysis workbench. Bioinformatics 25, 1189–1191. doi: 10.1093/bioinformatics/btp033
Wei, W., Lee, I.-M., Davis, R. E., Suo, X., and Zhao, Y. (2008). Automated RFLP pattern comparison and similarity coefficient calculation for rapid delineation of new and distinct phytoplasma 16Sr subgroup lineages. Int. J. Syst. Evol. Microbiol. 58, 2368–2377. doi: 10.1099/ijs.0.65868-0
Yang, Z. (2007). PAML 4: phylogenetic analysis by maximum likelihood. Mol. Biol. Evol. 24, 1586–1591. doi: 10.1093/molbev/msm088
Yang, Z., and Nielsen, R. (2000). Estimating synonymous and nonsynonymous substitution rates under realistic evolutionary models. Mol. Biol. Evol. 17, 32–43. doi: 10.1093/oxfordjournals.molbev.a026236
Zhao, Y., and Davis, R. E. (2016). Criteria for phytoplasma 16Sr group/subgroup delineation and the need of a platform for proper registration of new groups and subgroups. Int. J. Syst. Evol. Microbiol. 66, 2121–2123. doi: 10.1099/ijsem.0.000999
Zhao, Y., Wei, W., Lee, I.-M., Shao, J., Suo, X., and Davis, R. E. (2009). Construction of an interactive online phytoplasma classification tool, iPhyClassifier, and its application in analysis of the peach X-disease phytoplasma group (16SrIII). Int. J. Syst. Evol. Microbiol. 59, 2582–2593. doi: 10.1099/ijs.0.010249-0
Keywords: plant pathogen, comparative genomics, average nucleotide diversity, taxonomy, multilocus sequence analysis, effector
Citation: Cho S-T, Kung H-J, Huang W, Hogenhout SA and Kuo C-H (2020) Species Boundaries and Molecular Markers for the Classification of 16SrI Phytoplasmas Inferred by Genome Analysis. Front. Microbiol. 11:1531. doi: 10.3389/fmicb.2020.01531
Received: 31 January 2020; Accepted: 12 June 2020;
Published: 10 July 2020.
Edited by:
Haiwei Luo, The Chinese University of Hong Kong, ChinaReviewed by:
Christophe Garcion, INRA Centre Nouvelle-Bordeaux-Aquitaine, FranceXavier Foissac, Institut National de Recherche pour l’Agriculture, l’Alimentation et l’Environnement (INRAE), France
Kenro Oshima, Hosei University, Japan
Amit Yadav, National Centre for Cell Science, India
Copyright © 2020 Cho, Kung, Huang, Hogenhout and Kuo. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Chih-Horng Kuo, Y2hrQGdhdGUuc2luaWNhLmVkdS50dw==