Inferring the Minimal Genome of Mesoplasma florum by Comparative Genomics and Transposon Mutagenesis

Several Mesoplasma species have been isolated from plants or insects since the 1980s (25, 30–34). Of these, we have obtained 13 M. florum strains available from culture collections (Table 1). Two strains originated in Florida, and 11 were collected on a longitudinal transect from Maryland to Colorado (Fig. 1A). This variety of environments, host organisms, and climates was expected to result in diversity of gene content in the genomes analyzed. The complete genome sequences of M. florum strains L1 and W37 have previously been reported (26). The genomes of the remaining 11 M. florum strains were sequenced by using a combination of the Illumina, Pacific Biosciences, and Sanger technologies. In total, seven genomes (L1, W37, BARC 787, MQ3, CnuA-2, MouA-2, and W23) were fully assembled, resulting in circular chromosomes without any ambiguous positions. The genomes of BARC 786, BARC 781, GF, W17, W20, and W12 were not unequivocally resolved during their assembly (containing one to four gaps) because of a total of five distinct repeated elements, either duplications or tandem repeats (Table 1). By keeping at least one copy of each problematic region in each final assembly, we estimate that virtually all genes are represented in these genomes and that no more than a few kilobase pairs are missing from the chromosomes of these strains. The genome size of the 13 sequenced strains was found to vary between 738,512 bp (BARC 787) and 830,640 bp (20), with a mean of 794.5 ± 25.6 kbp. As expected from the genome sequences of other members of the class Mollicutes (35), the GC content of every strain analyzed was found to be relatively low, with an average of 27.1% ± 0.2% (Table 1).

To apply a uniform procedure and minimize potential biases between annotations, gene prediction was performed for all of the M. florum genomes studied, including those of the previously sequenced L1 and W37 strains (see Data Sets S1 and S2 in the supplemental material, respectively). Every genome was predicted to contain 29 tRNA genes, as well as two rRNA gene loci, each encoding 5S, 16S, and 23S rRNAs. Between 651 and 740 protein-coding genes were predicted, depending on the M. florum strain analyzed, for a total of 9,074 putative proteins in the combined 13 genomes. The genomes showed a protein-coding gene to kilobase pair ratio of ~0.88, which is typical of bacterial genomes (36). All proteins were clustered on the basis of sequence homology (37). Genes conserved among all M. florum strains and M. capricolum as the outgroup were used to construct phylogenetic trees based on the Kimura distance model (38) (Fig. 1B) and maximum likelihood (39) (Fig. 1C). The two trees showed practically identical branch lengths and very similar topologies, diverging only slightly around the L1/GF/MouA-2 subgroup. Overall, the trees revealed two main groups that are, in turn, subdivided into smaller phylogenetic clusters. Although some closely related strains were isolated from nearby locations or similar organisms, the geographic origin and potential host were not sufficient to explain the phylogeny observed.

The predicted proteins in the M. florum strains were grouped into a total of 1,150 homologous gene clusters (Data Set S3). A core set of 546 clusters was observed in all representatives, resulting in the conservation of approximately 80% of the protein-coding genes in every strain (Fig. 2A and B). A majority (75.5%) of the 604 remaining gene clusters, also called the accessory genome or pangenome, was found in no more than three strains. Additionally, the number of gene cluster families discovered kept increasing as a function of the number of strains analyzed (Fig. 2A), suggesting that M. florum has an open pangenome and that gene diversity in this species was not fully explored by investigating these 13 genomes. Considering that an average of 23.5 gene families were found to be unique in each M. florum strain, several additional gene cluster families should be found simply by incorporating more genomes into the analysis.

To obtain a functional overview of the M. florum genomes, the genes from the core genome and pangenome were classified into clusters of orthologous groups (COG) functional categories (40). Most (74.7%) of the core genes could be associated with a COG category, while most (66.7%) of the accessory genes could not. Overall, more than half (53.4%) of the COG associations belonged to the translation (J), carbohydrate metabolism (G), replication/recombination/repair (L), general function prediction only (R), and amino acid transport metabolism (E) categories. The core genome was significantly enriched (Fisher test with Bonferroni correction, P < 0.05) for functions related to translation (COG category J) compared to the frequency of this same category in the entire genome (Fig. 2C). The distribution of COG categories in the pangenome was more diversified, and no category was found to be statistically significantly enriched or deprived when all M. florum strains were considered. Genes linked to carbohydrate metabolism (mainly phosphotransferase system [PTS] components and β-glucosidases) and genome maintenance, mostly in restriction/modification systems (COG G and L categories, respectively) varied greatly in number, depending on the strains, although they were generally more abundant in the accessory genome (Fig. 2D).

We next analyzed the genomic organization in M. florum strains by characterizing the relative positions of the core genes. We found that although some rearrangements can be observed, the genomes are mostly syntenic with a large number of genes conserved between the different strains and variable regions located within the same relative genomic loci in every strain (Fig. 3). Furthermore, the gene order conservation (GOC) scores (41) of the core genes calculated for every combination of two strains averaged 0.98 ± 0.06, where a GOC value of 1 indicates that all core genes are found in the same order between two strains. We also noted the presence of large weakly conserved or nonconserved regions, consistent with the presence of genomic islands. Regions with at least three consecutive noncore genes were further studied to determine their potential status as genomic islands resulting from horizontal gene transfer. An average of 13.9 ± 3.4 such regions were identified in each genome, representing 19.4% of all noncore gene clusters. One such cluster, detected only in strains W37, W12, and W20, was located near the 660-kbp position and contained 16 hypothetical proteins, as well as 1 protein predicted to be part of a phage tail fiber. Using BLASTP, we searched for protein homologs in other organisms. Of those, 12 had hypothetical protein homologs in spiroplasmas and entomoplasmas, which are the genera closest to mesoplasmas among the members of the class Mollicutes. We also detected an 11-gene cluster ranging from position 14,058 to position 24,472 in strain BARC 781. According to the BARC 781 genome annotation, this cluster contained genes related to type IV secretion systems and a predicted mobile element protein, whereas in the other strains this region was enriched in PTS component proteins. Interestingly, the genomic location between translation initiation factor 2 and dnaJ was also found to be highly variable, even among closely related strains. Depending on the strain, this locus contained between 1 and 21 genes encoding proteins annotated as PTS components, hypothetical proteins, phage-related proteins, restriction-modification systems, and transcription regulators. These results suggest that horizontal gene transfer events have occurred within the M. florum species, shaping its genomic landscape.

In addition to comparative genomic data analysis, we also performed transposon mutagenesis to identify essential genes in M. florum L1. A collection of 2,806 mutants in which insertions occurred, on average, every ~280 bp across the genome resulted in the interruption of 430 of 720 genes (Data Set S4). No transposon was observed in the remaining 290 genes, which could be essential genes or simply have been missed given the transposon insertion density. The probability of observing no transposon insertions within a gene was calculated on the basis of the transposon insertion density and the length of each individual coding sequence (CDS) as previously described (16, 42), assuming that the probability of having N insertions in a gene of length L follows a Poisson distribution (Data Set S1). Although the transposon insertion density in our work is superior to what was reported in other gene inactivation studies involving Mollicutes (12, 43), the average probability that a gene could have been missed by chance in our experiment is ~10%, corresponding to ~69 genes for the entire M. florum L1 genome. However, this distribution is strongly skewed toward small genes, with more than half (37/69) of the potentially missed genes having a length of <400 bp. For simplicity, all 290 genes that were not interrupted by a transposon were nonetheless considered putatively essential despite these limitations.

The gene interruption and conservation data for the L1 strain were combined to identify genes that are most likely to be important for M. florum (Fig. 4A). Most of the genes presumed to be essential in the L1 strain on the basis of transposon mutagenesis were also conserved across all M. florum strains, with only 25 that were not associated with the core genome (Fig. 4B). A putative function could be attributed to only 5 of these 25 noncore but essential genes, of which 4 could be associated with a COG category and the remaining 20 were annotated as hypothetical proteins (Data Set S1). The 110 genes interrupted by transposons and absent from the core genome represent interesting first-step candidates for genome streamlining.

Three reduced versions of the M. florum L1 genome were designed on the basis of the gene conservation and transposon mutagenesis data (Fig. 5A). The first scenario excluded all noncore genes that were also interrupted by a transposon (610 genes remaining), the second contained only the core genes (585 genes from the 546 conserved gene cluster families in M. florum L1), and the third included only the putatively essential genes (290 genes). Both rRNA loci were included in these genome configurations, although only a single copy could be sufficient to sustain growth (2, 44). In each scenario, the approximate genome size was estimated by removing the CDSs of the candidate genes, while all non-CDSs were kept since most of the promoters and regulatory sequences in M. florum remain to be identified. A fourth genome reduction strategy was prepared by including all of the M. florum L1 protein-coding genes having an ortholog in M. mycoides JCVI-syn3.0, currently the closest approximation of a minimal genome. Interestingly, orthologs were identified for 401 of the 585 M. florum core genes. The M. mycoides JCVI-syn3.0-inspired genome reduction approach also included 5 genes identified as nonessential and noncore, while 57 genes marked as essential in M. florum L1 were absent.

We compared the COG category occurrences in the genome reduction strategies, as well as in the original M. florum L1 and M. mycoides JCVI-syn3.0 genomes (Fig. 5B). Globally, the proportions of the COG categories were similar in the genome reduction scenarios. Genes related to translation (COG category J) were more represented than any other. Genes related to lipid metabolism (I) made up the second most conserved category, although no gene predicted to perform lipid degradation or elongation was detected in M. florum. This category instead contained multiple acyltransferases, which could be used to anchor lipoproteins to the membrane (35), as well as a number of NADH and short-chain dehydrogenases. Overall, the carbohydrate metabolism (G) category contained the largest number of accessory genes since even in the most conservative proposition, 33.8% of these genes were identified as dispensable. Of those, 72.7% were PTS transport components or 6-phospho-β-glucosidases whose functions were duplicated.

M. florum strains were grown at 34°C in ATCC 1161 medium (27). Genomic DNA (gDNA) extraction was performed with the Quick-gDNA MiniPrep kit (Zymo Research) for the preparation of Illumina libraries and with Puregene Yeast/Bact. kit B (Qiagen) for PCR-free Sanger sequencing and Pacific Biosciences libraries. Extractions were made in accordance with the manufacturer’s instructions, except that M. florum cells were washed in a resuspension buffer (8 mM HEPES, 272 mM sucrose; pH 7.4) prior to gDNA extraction.

Genomes were assembled as previously described for the M. florum W37 genome (26). For each strain, two Illumina libraries were prepared, one with 200- to 250-bp inserts and the second with 450- to 750-bp inserts, both sequenced in paired-end reads of 144 and 100 bp, respectively. Pacific Biosciences RS libraries (obtained by C₂ chemistry) with inserts of >5 kbp were also prepared. Error correction of the Pacific Biosciences reads was performed by using the Illumina reads, and all sequences were subsequently assembled with Roche gsAssembler version 2.6 and Ray version 2.1.0 (48). The two assemblies were merged and manually inspected before being completed (or completely scaffolded; Table 1) by Sanger sequencing of selected PCR products (the primers used are available on request). PCRs were performed with VeraSeq DNA polymerase (Enzymatics). The routine PCR conditions used were 30 s at 95°C; 30 cycles of 10 s at 95°C, 30 s at the appropriate annealing temperature, and 30 s/kbp at 72°C; and 2 min at 72°C. PCR products were purified by solid-phase reversible immobilization bead capture with Agencourt AMPure XP magnetic beads (Beckman Coulter, Inc.) (49). Sanger sequencing reads were generated by the Plateforme de Séquençage et de Génotypage of the Centre de Recherche du Centre Hospitalier de l’Université Laval. Pacific Biosciences sequencing was performed by the Yale Center for Genome Analysis and at the Centre d’innovation Génome Québec et Université McGill. All of the genomes were annotated with the RAST server (50) and FIGfam Release 70.

The annotated genomes were analyzed with get_homologues (51) (v.2.0) to identify the core genome and pangenome by the COGtriangle (37) method for bidirectional best-hit determination. The comparison of the M. florum strains and M. capricolum subsp. capricolum ATCC 2734 (NCBI RefSeq accession no. NC_007633.1) to identify the protein set used for the phylogenetic tree creation was performed by the same method, as was the comparison of M. florum L1 and the synthetic bacterium Syn3.0 (GenBank accession no. CP014940.1). Genome synteny was determined by using the core protein cluster coordinates and visualizing their locations across all of the genomes by using the GMV genome browser (v.1e-93) (52). GOC scores were determined by using the core genes as previously described (41). Genes separated by fewer than five core genes were considered contiguous (41). For the functional analysis, the latest COG database (53) was downloaded from the NCBI ftp website (ftp://ftp.ncbi.nih.gov/pub/COG/COG2014/data/). The proteins predicted in the 13 M. florum genomes were compared to the database with BLAST (54). The results were filtered on the basis of the E value with a stringent threshold of 1e-10. Available conversion tables from the NCBI ftp website were used to convert matching proteins to COG identification numbers to functional categories. Enriched categories were determined with the Fisher exact test and Bonferroni correction for multiple tests by comparing the frequency of genes in a COG category in the M. florum core or accessory genome with that in the entire genome.

An alignment of the amino acid sequence of 412 conserved proteins of the 13 sequenced M. florum strains and M. capricolum was made with ClustalO (v.1.2.1) (55). Unaligned and low-confidence regions were removed with Gblocks (v.0.91b) (56) to produce a sequence matrix of 138,476 amino acid sites. Both phylogenetic trees (Fig. 1B and C) were made from this alignment with SeaView (v.4.6.1) (57). The distance tree was generated by neighbor joining with BIONJ (58) and the Kimura distance model (38). The maximum-likelihood tree was generated with PhyML (v3.0) (39) and the LG evolutionary model (59). Bootstrap values were calculated by using 100 regular bootstrap replicates.

Tn5 transposomes were assembled in vitro with EZ-Tn5 transposase (Epicentre) without Mg²⁺ as recommended by the manufacturer. The transposon DNA was obtained from the digestion of the pTT01 plasmid with restriction enzyme PvuII and conferred tetracycline resistance. The transposomes were transformed via electroporation as described by Matteau et al. (27) with a voltage of 2.5 kV. Insertion mutants were selected on ATCC 1161 solid medium supplemented with 15 µg/ml tetracycline. Colonies were picked as they became visible, with certain mutants growing more slowly than the parental strain. The transposon insertion site was determined by Sanger sequencing of gDNA. Genes that contained at least one transposon insertion were considered nonessential.

The complete genome sequences and annotations of the 11 newly sequenced M. florum strains are available in GenBank under the following accession numbers: BARC 781, CP022511; BARC 786, CP022510; BARC 787, CP022514; CnuA-2, CP022513; GF, CP022509; MouA-2, CP022508; MQ3, CP022512; W17, CP022507; W20, CP022506; W23, CP022505; W12, CP022432. The comparison of NCBI annotations and gene predictions used in our study for strain L1 is presented in Data Set S1. Gene predictions for the 12 other M. florum strains can be found in Data Set S2. The composition of the M. florum COG can be found in Data Set S3. The transposon insertion sites observed in M. florum L1 are provided in Data Set S4.