Exploring the symbiotic pangenome of the nitrogen-fixing bacterium Sinorhizobium meliloti

General features of AK83 and BL225C genomes

Strain AK83 has the highest percentage of ORFans (orphan ORFs) (5.63%), defined as those genes with no detectable similarity with other genes in any other organism [20], and the lowest number of ORFs coding for proteins with homology to COGs, InterPro, GO and Rhizobase entries, while strain Rm1021 shows the highest number of transposases and Insertion Sequences (152) compared to AK83 and BL225C (135 and 76, respectively). Other relevant features of general importance are those related to the environmental sensing and transport: we noticed that BL225C and AK83 strains have more Type IV secretion systems-related proteins than the reference strain Rm1021. However, a similar number of ABC transporters-related proteins and of Type III secretion systems-related proteins were found across the three genomes, tough no fully functional Type III secretion systems were detected, as previously noticed for Rm1021 strain [9]. Finally, two-component signal transduction systems related proteins are slightly higher in Rm1021 and BL225C strains than in AK83 (129, 134 and 125, respectively).

Defining core and accessory S. meliloti genome

By comparing the 19095 CDSs, found in the three genomes, a set of 7824 orthologous groups was identified; a subset of 5124 was conserved across all the three genomes and accordingly defined as the core genome of S. meliloti species. The remaining 2700 orthologous groups were defined as members of the accessory genome for these three genomes. The strain with more unique genes is AK83, with 843 exclusive groups, while BL225C and Rm1021 have 469 and 602 exclusive groups, respectively (Figure 1). In Additional file 1 the full list of core and accessory proteins is reported.

When the very recently published SM11 genome [10] was added to our proteome set of the three strains (AK83, BL225C, Rm1021), the core genome contained 5075 orthologous groups, with the loss of 49 groups, suggesting a certain stability of the core genome size in this species, while the accessory genome comprised 3810 orthologous groups.

In order to define possible differences in functions encoded by the core and/or the accessory genomes and by the different strains, each protein was assigned to a COG category (Additional file 2) and the abundance of each COG category was plotted (Figure 2, Additional file 3). Statistically significant differences between core and accessory genome were found only for COG category L (DNA replication, recombination and repair) and for proteins with no assigned COG (X): in these two categories, the accessory genome is enriched, especially in the category X. Similar enrichment in CDSs with no assigned function has been previously reported in the accessory genome of other organisms [21] as well as in the S. meliloti plasmid pSmeSM11a [22]. For other COG categories, no statistically significant difference was found, though a higher representation of all assigned functions was found in the core. Finally no significant difference of COG categories between the three strains was found (Additional file 3).

Structural genomics

While BL225C contains three replicons as Rm1021, AK83 is composed by five circular replicons, corresponding to the chromosome, pSymA and pSymB, which are also present in the genome of Rm1021 and BL225C and two new small replicons 1 and 2, respectively 0.26 Mbp and 0.07 Mbp in size as schematized in Figure 3. The genomic structures of the fully assembled complete genomes were compared with a full-scale genomic alignment (see materials and methods). The chromosome and the pSymB chromid are characterized by a high resistance to genome rearrangements, with an almost perfect shared synteny, with only few insertions in the chromosome of strain AK83 and few rearranged regions of the chromid pSymB. The other replicons showed indeed lower degrees of synteny: in particular Plasmid2 of strain AK83 had no region of similarity with the other replicons of strain Rm1021 and BL225C (and with other plasmids available in the NCBI database). Concerning the symbiosis-required megaplasmid pSymA, a very low degree of synteny was observed indicating an increased rate of rearrangements for pSymA, and indeed evidences for rearrangements were noticed, since at least three fragments of Plasmid1 showed an high degree of similarity with the Rm1021 pSymA (ca. 47 kbp), while just one fragment was found to be highly similar to the pSymA of BL225C (ca. 27 kbp). These data suggest that AK83 Plasmid1 may be derived from or represent an evolutionary step toward the megaplasmid pSymA. Interestingly, a fragment 8 kbp long of Plasmid1 showed similarity with another symbiotic-related plasmid (Sinorhizobium fredii NGR234 plasmid pNGR234b). On the other hand, megaplasmid pSymA of strain BL225C compared to Rm1021 showed a higher number of syntenic regions and a lower number of inversions than AK83 (6 regions over 27 for the former and 11 regions over 24 for the latter). The genomic structure of the newly sequenced S. meliloti strain SM11 was also analyzed (data not shown); as expected SM11 pSmeSM11c (replicon carrying symbiotic functions, analogous to Rm1021 pSymA) is the most diverse replicon in comparison with replicons of strains AK83 and BL225C. No significant homologies were found between the small plasmids of strain AK83 (Plasmid 1 and Plasmid 2) and strain SM11 genome, as well as no significant hits between strain SM11 small plasmids (pSmeSM11a and pSmeSM11b) and AK83 and BL225C genomes. However, the presence on strain SM11 pSmeSM11a plasmid of two regions of 15 and 10 kbp syntenic to megaplsmid pSymA of strain Rm1021 was confirmed [22].

The location of transposases and insertion sequences was analyzed. Indeed in most of the cases, transposase encoding genes were enriched in regions carrying the traces of genomic rearrangements (with a frequency of 1.53, 1.32 and 0.59 IS/10 kb in RM1021, AK83 and BL225 respectively), than in highly syntenic regions where the corresponding frequencies are 0.19, 0.11 and 0.06 IS/10 kb in the three genomes. As reported in Figure 3, 17 transposases were found in a 80 kbp insertion of the chromosome of strain Rm1021 (3'357'000 to 3'440'000), as well as in the biggest five non-syntenic regions of the AK83 chromosome. In the pSymA megaplasmids, 13 transposases were found in the non-syntenic region of strain Rm1021 (200'000 to 390'000), 14 and 15 transposases were found in strain AK83 in two non-syntenyc regions (300'000 to 457'000 and 1'102'000 to 1'198'000) and in strain BL225C 17 and 2 transposases were found in two non-syntenyc regions (630'000 to 934'000 and 144'000 to 1'482'000). Finally, 4 transposases were found flanking the two regions of strain AK83 Plasmid1 that are also present on megaplasmid pSymA in strain Rm1021.

Accessory genome and symbiosis-related functions

Looking at the genetic content of the symbiotic megaplasmid pSymA of Rm1021 in comparison with AK83 and BL225C, previous findings using CGH were confirmed [6, 23]: in fact, many genes harbored by Rm1021 pSymA were found missing in the genomes of the two other strains (Figure 4). Moreover, two regions of Rm1021 pSymA (200'000 to 390'000 and 679'000 to 708'000), predicted to be a large part of the so-called microaerophilic gene set [24], were not present in strain AK83 and to a lesser extent also in strain BL225C.

We then focused on genes involved in some aspects of the symbiotic process, that we identified using a data mining strategy combining different sources of information, using the following approach: for each orthologous group having a predicted link to a NodMutDB and/or Rhizobase member (see materials and methods) the related literature was retrieved and analyzed to speculate its actual role in symbiosis. This approach was combined with other annotation sources (such as KEGG and Interpro) in a dedicated data mining procedure (Additional file 4); together with this approach, also those proteins with names containing symbiosis-related terms (e.g. fix, nif, nod) were retrieved. In particular we were interested in genes that are differentially present in the three strains and that could be related to the different symbiotic phenotypes of these strains.

By using the strategy described above, we identified, among the accessory genes, those that have been implicated in the symbiotic process in rhizobia through experimental work. Within a total number of 290 orthologous groups retrieved, 61 of them were found to belong to the accessory genome (21%) (Additional file 5).

In conclusion, it can be noted that 48 orthologous groups are missing in the symbiosis defective AK83, which is, in fact, a large proportion of the accessory genes related to symbiosis (79%).

The symbiosis-related panregulon

To further investigate the genomic differences that may be related to the variable symbiotic phenotypes of AK83, BL225C and Rm1021 strains, the predicted regulons of a series of symbiotic transcription factors were analyzed. A set of eight transcriptional regulators related to symbiotic interaction was chosen, based on the knowledge of their DNA binding sites (Additional file 6) and on regulon information on closely related rhizobial species that may share the same binding site with S. meliloti. The set included transcriptional regulators involved in root exudates perception and early nodulation steps (NodD, NolR), microaerophilic adaptation (FixK), nitrogenase synthesis (FixJ, NifA), iron uptake (Fur), EPS biosynthesis (ChvI) and plant invasion competition (NesR) (Table 3). It should be noticed that the FixK regulated genes only partially overlap with those found to be expressed under microaerophilic conditions [24], which may be under control of other regulators. For each transcriptional regulator, genes putatively regulated and present in the genomes of strains Rm1021, AK83 and BL225C were sorted out by HMM scanning and the core (conserved in all strains) and the accessory (variable among strains) putative regulons were defined (Figure 5a). We defined the panregulon as the totality of gene families controlled by a specified transcription factor in a certain number of genomes, in analogy to the term pangenome [12], and which is formed by the core and the accessory regulons. The putative panregulon varies in sizes from 101 (FixK) to 6 (NolR) orthologous groups (i.e. genes); since some of the regulated targets could be part of an operon or be additional regulators, the actual size of the predicted panregulon is definitely under-estimated. NolR and NesR show a very little panregulon (with a core regulon of 3 and 5 orthologs, respectively); the accessory regulons of the other transcriptional regulators are very large and variable, accounting for 31-79% (average 55%) of all panregulons. The occurrence of wide accessory regulons could be due to the absence in one or two strains of the target genes or to the absence of the regulatory upstream sequences when the genes are present. The absence of target genes is the most frequent case ranging from 50% (NesR) to 100% (NolR) of the targeted genes, with an average value of 69% (Figure 5b). On the contrary the variability of DNA binding sites upstream CDSs (genes were still present in a given genome but were not putatively regulated by that transcription factor) is less frequent with an average value of 31%. The list of all genes sorted out as putatively regulated by the selected transcriptional regulators is reported in Additional file 6. The composition of the accessory regulons in terms of un-annotated targets was calculated counting the number of CDS with no COG classification (Additional file 7) resulting in an average percentage of 49%, while for the core regulons the percentages of un-annotated targets is lower (21%).

			Genes regulated
Transcription factor	Symbiotic process	Reference	Rm1021	AK83	BL225C
NodD1	Flavonoid perception	[40, 68]	10	13	12
ChvI	EPS biosynthesis	[45]	54	65	52
FixK	Microaerophilic adaptation	[41]	54	61	54
FixJ	Nitrogenase synthesis and functioning via nifA	[69]	26	31	35
NifA	Nitrogenase biosynthesis	[41]	35	48	42
Fur	Iron uptake	[43]	9	8	9
NolR*	Optimization of nodulation, bacterial growth on solid medium, survival under stress conditions, and conjugative transfer of plasmids	[70]	3	6	3
NesR	Competition for plant nodulation	[44]	6	5	7

* gene disrupted by a frameshift mutation in Rm1021 [46]

The genes of the predicted regulons were then divided into 5 functional groups: electron transport, symbiosis, nitrogen metabolism, others or un-annotated. In Figure 6a the eight regulons are showed; in many cases putatively regulated genes could be matched with experimental data: for NodD1 the experimental regulation of nodL, nodF, nodA, nodM and syrM [40] was confirmed by the analysis; for FixK the regulation of fixNOQP ₁ , fixNOQP ₂ , fixGIS, arcABC, napEFDABC, norBCE, cycB2 and degP4 was confirmed [41]; for NifA the regulation of nifHDKEX, fixABCX operons and the nifB genes [41] were also confirmed. Experimental confirmation of regulatory predictions by NolR was also found for nodA, nodD1, nodD2 and nodM [42], by Fur for sitA [43], by NesR for the metHK operon and the ahcY gene [44] and by ChvI for ropB [45]. Interestingly, even if the nolR gene is disrupted by a single base insertion in strain Rm1021 [46], the NolR binding sites in front of nodD1, nodD2 and nodA are maintained, suggesting that the inactivation may be relatively recent. Co-regulation by different factors was observed on several genes; in particular 7 genes were putatively regulated by more than one regulator (nodD1, nodD2, nodA, nodM, the AK83 copy of nodM and two other uncharacterized proteins) and a co-regulation by 4 out of the 8 selected regulators was also found (NolR on NodD1 and FixJ on FixK), indicating that some of the symbiotic regulators are linked together in a network that ensures the coordination of the expression of the genes required during infection and nitrogen fixation. The FixK regulator is predicted to control the highest number of genes involved in the symbiotic process, as well as those involved in nitrogen metabolism and electron transport needed in the microaerophilic environment of the bacteroid [24].

Concerning the variability in the predicted regulons (Figure 6b), it is evident that AK83, shows some differences in the regulatory networks: in AK83 purB is apparently not controlled by NifA; nodM is regulated by NolR only in strain AK83; in AK83 FixK controls napE and napF. It should be noted that these differences were not observable from patterns of gene presence/absence, illustrating the added value of regulon prediction in comparative genomics.

Bacterial strains and culture conditions

BL225C, isolated in Italy in an alfalfa field and AK83, isolated in the Aral sea region, were deposited at the German Collection of Microorganisms and Cell Cultures (DSMZ) with accession codes DSM23914 for strain BL225C and DSM23913 for strain AK83. AK83 strain is also present, as original specimen after initial isolation, in the culture collection of All-Russia Institute of Agricultural Microbiology (RIAM, St. Petersburg, Russia). S trains were cultured on solid or liquid TY medium [49] with 0.2 g/liter CaCO₃ at 30°C.

Whole-genome shotgun sequencing and draft annotation

Total DNA was isolated from S. meliloti AK83 and BL225C cultures with a CTAB method according to the recommended protocols by JGI-DOE http://my.jgi.doe.gov/general/. Genome sequencing was performed at the Joint Genome Institute (JGI) (Walnut Creek, California, USA) using a combination of Illumina [50] and 454 technologies [51]. The 454 Titanium standard data and the 454 paired end data were assembled together with Newbler, version 2.3. The Newbler consensus sequences were computationally shredded into 2 kb overlapping fake reads (shreds). Illumina sequencing data was assembled with VELVET, version 0.7.63 [52], and the consensus sequences were computationally shredded into 1.5 kb overlapping fake reads (shreds). The 454 Newbler consensus shreds, the Illumina VELVET consensus shreds and the read pairs in the 454 paired end library were integrated using parallel phrap, version 4.24 (High Performance Software, LLC). The software Consed [53–55] was used in the following finishing process. Illumina data was used to correct potential base errors and increase consensus quality using the software Polisher developed at JGI (Alla Lapidus, unpublished). Possible mis-assemblies were corrected using gapResolution (Cliff Han, unpublished), Dupfinisher [56], or sequencing cloned bridging PCR fragments with subcloning. The gaps between contigs in the genomes of strain AK83 and BL225C were closed by editing in Consed, by PCR and by Bubble PCR (J-F Cheng, unpublished) primer walks. For strain AK83 a total of 968 additional reactions and 11 shatter libraries were necessary to close gaps and to raise the quality of the finished sequence; the final assembly is based on 279.6 Mb of 454 draft data which provides an average 31.3 × coverage of the genome and 426 Mb of Illumina draft data which provides an average 62 × coverage of the genome. For strain BL225C a total of 801 additional reactions were necessary to close gaps and to raise the quality of the finished sequence. The final assembly is based on 290.2 Mb of 454 draft data which provides an average 27× coverage of the genome and 308 Mb of Illumina draft data which provides an average 44× coverage of the genome. For both genomes the gene model prediction and draft annotation was generated using Prodigal [57]. Sequences and annotation can be accessed through the JGI web site at the addresses http://genome.jgi-psf.org/sinma/sinma.home.html for AK83 and http://genome.jgi-psf.org/sinmb/sinmb.home.html for BL225C; both strains are being submitted in GenBank with the following master records [Genbank: NZ_AEDG01000000, Genbank: NZ_AEDH01000000] for BL225C and AK83, respectively.

Annotation

Annotation was performed again on the three genomes using Blast+ 2.2.23 [58] and InterproScan 4.6 [59]. A bidirectional best blast hit (BBH) approach was used to annotate all the predicted proteins in the three genomes using the following three databases: NCBI nr, downloaded on May 18, 2010, all the Rhizobase http://genome.kazusa.or.jp/rhizobase/ proteomes, downloaded on June 1^st, 2010 and the KEGG database http://www.genome.jp/kegg/, downloaded on May 26, 2010; for the first two databases an E-value threshold of 1e-10 was applied, while for the KEGG database a threshold of 1e-50 was applied. For the domain scan using InterproScan the Interpro database release 27 was used. All the results were linked to literature using the Interpro database, the UNIPROT database http://www.uniprot.org/ release 2010_06, the KEGG database, the Rhizobase Bibliome and the nodMutDB [5]. The number of predicted transposases and IS was inferred with a keyword search in the annotated protein set, while the number of putative Type III and Type IV secretion systems-related proteins, two component systems and ABC transporters was inferred searching specific interpro domains in each strain proteome.

Structural Genomics

Genomic alignment between the complete genomes was generated using the megablast algorithm [58] retaining only hits of more than 10 kbp; the starting point of the replicons was changed in order to generate a clearer syntenic map with the Artemis Comparison Tool (ACT) from the Artemis suite [60]. The similarity of the two smaller plasmids of strain AK83 with other known plasmids was inspected using the megablast algorithm on the entire nucleotide NCBI database and retaining only those hits bigger than 5 kbp. The structure of the reference genome was compared to the newly sequenced genomes using the nucleotide sequence of each protein with a nucleotide Megablast with a 50% identity threshold: results were visualized using DnaPlotter from the Artemis suite [60].

Orthology

Since the actual magnitude of a pangenome is computable only by sequencing each strain of the desired species [61], here we will refer to the term pangenome not as the full gene complement of the S. meliloti species, but only to the observable one.

The three proteomes were clustered into orthologous groups using a BBH approach through InParanoid 4.0 [62] and MultiParanoid [63]. The BLOSUM 80 matrix was used during the InParanoid run, while the "unique" flag was applied in the MultiParanoid run; MultiParanoid could be reliably used since the three genomes are evolutionary closely related. The pangenome size of the species Escherichia coli, Bacillus anthracis and Streptococcus pneumoniae were determined picking three complete proteomes at random from the NCBI public database and applying the same approach as S.meliloti; ten repetitions with different genomes were performed to calculate the average pangenome size.

Functional Enrichment

To elucidate if the accessory genome was enriched in a particular function, the proportions of the COG categories [64] in the core and accessory genome were compared; to give statistical significance to the difference an enrichment analysis was performed, in a similar way as in Brilli et al. [19]; one million random samplings were performed and the COG proportions of each sample was compared to a sample from the whole genome. P-values below 0.05 were considered significant.

Promoter prediction

Promoter prediction was performed by taking the nucleotide sequences in region -600 + 100 around the predicted gene start of all protein coding sequences in the three genomes. HMMer 3.0 [65] was used to build the promoter box HMMs (hmmbuild program) and to scan the promoter regions (hmmsearch program). The input alignments that generated the HMMs were retrieved from MEME [66] scans on sequences derived from literature. HMM scan was performed by switching off all the heuristic filters, collecting all the hits and calculating the score mean and standard deviation; after verification of the normality of the score distribution using Past [67] only those hits having a score greater than 3 standard deviation above the mean value were retained. For the prediction of the two FixJ DNA binding motifs, the results obtained were merged together.

Data storage and scripting

All the collect data from annotation, orthology and promoter prediction were stored in a MySQL relational database and linked together in a proteome-centric way. All the analysis were performed using ad-hoc Python scripts, taking advantage of the BioPython and SciPy packages.