Impacts of Natural Selection on Evolution of Core and Symbiotically Specialized Genes in the Polytypic Species Neorhizobium galegae

1. Introduction

Root nodule bacteria (rhizobia), N₂-fixing symbionts of leguminous plants, represent a convenient model for the evolutionary genetics of symbiotic organisms. As in other bacteria, rhizobia genomes are composed of conservative core and variable accessory parts [1]. Core genes are responsible for housekeeping functions (basic metabolism, cell development and reproduction, template processes) while the accessory genes encode for various adaptive functions including symbiotic interactions with the leguminous plants. Rhizobia genomes are subjected to the multilevel evolution based on modifications of: (i) symbiotically specialized (sym) genes resulted in formation of polytypic species composed of host-specific biovars; (ii) core genes resulted in formation of cryptic (twin) species [2].

Previously we studied the genome diversification in Rhizobium leguminosarum, a polytypic species which includes two biovars: viciae (symbionts of legumes from tribe Fabeae, genera Lathyrus, Lens, Pisum, Vavilovia, Vicia) and trifolii (symbionts of genus Trifolium from the tribe Trifolieae) [3]. Cross-inoculation between these biovars is limited and results in the non-N₂-fixing nodules, which are usually underdeveloped and morphologically abnormal [4]. In R. leguminosarum, genomes are composed of circular chromosomes and several plasmids, one of them (pSym) having size 200-500 kb harbors the majority of sym genes. They include nod genes encoding for synthesis of nodulation-inducing lipo-chito-oligosaccharidic Nod factors, and nif/fix genes encoding for the nitrogenase synthesis and operation. We suggested that in R. leguminosarum, evolution of sym genes is implemented under impacts of the host-induced natural selection [3,5] while mechanisms for the core gene evolution remain obscure.

In the presented paper we compare the evolutionary dynamics of sym and core genes in goats’ rue rhizobia (Neorhizobium galegae), a polytypic species differentiated into host-specific biovars orientalis and officinalis – N₂-fixing symbionts of Galega orientalis and G. officinalis, respectively. In contrast to R. leguminosarum, N. galegae biovars cross-inoculate their hosts readily resulting in the morphologically normal although non-N₂-fixing nodules. Majority of N. galegae sym genes are located on chromids having sizes over 1600 kb. These circular replicons have a plasmid type repABC system combined with the core genes which are typically located on bacterial chromosomes, including tRNA and rRNA encoding genes [6].

Previously we demonstrated [7] that populations of N. galegae bv. orientalis collected at the North Caucasian region are more polymorphic for sym and core genes than N. galegae bv. officinalis populations. This difference apparently reflects the diversity of respective host plant species which is sufficiently higher in G. orientalis than in G. officinalis. Difference between two N. galegae biovars for nif/fix genes was much more pronounced than for nod genes since the host specificity of compared biovars pertains N₂ fixation, not nodulation activity [7].

The presented paper demonstrates that in N. galegae, impacts of driving/purifying selection (dN/dS) on gene evolution differ sufficiently in the core and sym genes and are biovar-specific. Strict phylogenetic congruence of core and sym genes was revealed in N. galegae, reflecting the location of sym genes on non-transmissible chromids. In spite of this congruence, mechanisms of evolution are different in core and sym genes, as it was demonstrated by analysis of correlations between p-distance and dN/dS values. These correlations as well as analysis of core gene ontology groups allow us to suppose that in bv. officinalis evolution is more dependent on adaptation to endosymbiotic niches than in bv. orientalis.

Analyses of the presented data completes some important gaps in our knowledge on the tradeoff between polymorphism of genes and impacts of natural selection on their evolution. Specifically, we suggest that in ecologically versatile organisms, such as N. galegae bv. orientalis, driving selection enhances the nucleotide gene polymorphism since the newly emerged alleles may coexist stabile with the preexisted ones. However, in ecologically restricted organisms such as N. galegae bv. officinalis, driving selection does not influence or even decreases the gene polymorphism since new alleles supported by this selection replace (push away) the preexisted alleles from the genepool.

2. Results

2.1. Gene polymorphism and natural selection

We demonstrat that nucleotide polymorphism in polytypic species N. galegae depends on driving/purifying (dN/dS-measured) natural selection which operates in a biovar-specific manner and differs in the core and sym genes (Table 1). The maximal impacts of driving selection (dN/dS > 1) were revealed in the high-polymorphic core genes of bv. orientalis, the minimal impacts – in the low-polymorphic sym genes of bv. officinalis. A complete phylogenetic congruence between sym and core genes was revealed in N. galegae using the Average Nucleotide Identity (ANI) analysis (Figure 1) reflecting the parallel and possibly inter-dependent processes of microevolution and cryptic speciation dependent on modifications of sym and core genes, respectively.

Analysis of correlations between nucleotide polymorphism (p-distance) and the driving/purifying selection (dN/dS) impacts suggests (Table 2) that this selection implements different roles in evolution of core and sym genes which depend on the N. galegae biovars. For core genes, driving selection results in a marked increase of polymorphism in bv. oreintalis (indicated by a highly significant positive correlation between p-distance and dN/dS), but this increase is much less evident in bv. officinalis (indicated by a significantly lower although positive correlation between p-distance and dN/dS). For sym genes, natural selection does not influence polymorphism in bv. orientalis (no correlation between p-distance and dN/dS values) but results in a decreased polymorphism in bv. officinalis (negative correlation between these values).

Importantly, fraction of polymorphic genes in the total genepools is higher in bv. orientalis than in bv. officinalis for sym genes: 38.5±7,8% and 7.7±4,3%, respectively (t_St = 3.46; P₀ < 0,01). However, for core genes these fractions do not differ: 74.6±0.70% and 75.9±0.82%, respectively (t_St = 1,08; P₀ > 0,10). Importantly, fractions of polymorphic genes is higher for core than for sym genes in both biovars suggesting a strong purifying selection for sym genes.

Differential impacts of natural selection on polymorphism of core and sym genes are evident in bv. officinalis (r values differ significantly), not in bv. orientalis (r values do not differ) suggesting that adaptive impacts of these genes are biovar-dependent (Table 2).

2.2. Gene ontology analysis

Previously we suggested that in R. leguminosarum, evolution of sym genes occurs under the impacts of host-controlled natural selection, while the factors responsible for evolution of core genes remain unclear [3,8]. For addressing these factors, we used the Gene Ontology tools providing the functional annotation of core genes that vary sufficiently for nucleotide polymorphism (p-distance) and for susceptibility to driving/purifying selection (dN/dS) [9]. The polymorphic core genes were distributed into a range of Gene Ontology Groups (GOGs) contrasting for p-distance or for dN/dS in which deviations of these parameters from their average values of GO-enrichment exceed the standard deviation (1.509 for dN/dS and 1.519 for p-distance). This approach allowed us to distribute 782 core genes which are polymorphic in both N. galegae biovars into 76 GOGs comprising four clusters with contrasting p-distance and dN/dS values: (i) GOGs are elevated over average in both biovars, orientalis and officinalis (Ori+Off+); (ii) GOGs are elevated over average in bv. orientalis but are below average in bv. officinalis (Ori+Off–); (iii) GOGs are elevated over average in bv. officinalis but are below average in bv. orientalis (Ori–Off+); (iv) GOGs are below average in both biovars (Ori–Off–). For the statistical analysis, these four clusters were established independently for p-distance and dN/dS values and are presented as CP-I…CP-IV and CS-I…CS-IV, respectively.

We demonstrated that clusters CP-IV and CS-IV in which both analyzed parameters (p-distance and dN/dS) are below average in both N. galegae biovars (Ori–Off–) are most numerous suggesting that purifying selection (dN/dS < 1) resulting in a decreased nucleotide polymorphism represents an important factor of core gene evolution (Table 3). However, tradeoff between gene polymorphism and natural selection differs greatly in the analyzed clusters: representations (%) in the total pool of 782 polymorphic genes for Ori+Off+ and Ori–Off+ clusters are higher for p-distance than for dN/dS (CP-I > CS-I; CP-III > CS-III), while in Ori+Off– and Ori–Off– clusters, gene representations are higher for dN/dS than for p-distance (CS-II > CP-II; CS-IV > CP-IV) (Figure 2).

Vertical axis shows the representations (in %, with standard errors) of each cluster in the total pool of 782 analyzed genes (data from Table 3 are used, sizes of columns are given in Table S1 in Supplement).

Analysis of COG composition enables us to reveal several regularities in the functional diversity of core genome which pertain its operational (involved in cellular metabolism) and informational (involved in template processes) components. Specifically, analysis of GOGs identified on the basis of p-distance demonstrates (Table S2 in Supplement) that low-polymorphic genes responsible for metabolism of N-compounds (nucleosides, amino acids) are assigned to CP-II, while the low-polymorphic genes for lipid and oligosaccharide metabolism are in CP-III. This functional segregation of core genes may reflect a predominant dependence on symbiotic N₂ fixation typical of bv. officinalis: its host, G. officinalis grows at the limited locations under unfavorable edaphic conditions restricting the bacteria survival in soil and increasing their dependence on the nodular niches. This dependence is presumable much lower for bv. orientalis which is distributed in favorable locations occupied by its host (G. orientalis) and therefore this biovar persist in soil more stably than bv. officinalis.

As expected, genes which are less variable in both biovars (Ori–Off– clusters) are associated with highly conservative template processes. Specifically, genes encoding for translation are revealed in CP-IV while genes for replication, transcription, translation and DNA repair are in CS-IV (Tables S2–S4 in Supplement).

In order to address the tradeoff between nucleotide polymorphism and impacts of natural selection on the core gene evolution, we analyzed distributions of GOGs among the clusters identified using p-distance or dN/dS values (Table 3). For assessing statistically the coincidence of these distributions we calculated separately for two N. galegae biovars the frequencies of GOGs with elevated (over average) values of dN/dS or of p-distance within the GOGs with elevated or decreased values of p-distance or dN/dS, respectively. Analysis of data on the ecologically versatile bv. orientalis demonstrates that: (i) frequency of GOGs with elevated dN/dS values is significantly higher among GOGs with elevated p-distance than among GOGs with decreased p-distance (Figure 3a); (ii) frequency of GOGs with elevated p-distance is higher among GOGs with elevated dN/dS than among groups with decreased dN/dS (Figure 3b). Therefore, for ecologically versatile bv. orientalis, in both reciprocal comparisons of GOGs contrasting for p-distance and dN/dS values, “High-in-High” frequencies exceed “Low-in-High” frequencies suggesting that driving selection is responsible for an elevated core gene polymorphism. However, in the ecologically restricted bv. officinalis, no such difference was revealed suggesting that tradeoff between the core gene polymorphism and the natural selection impacts on the gene evolution depends on the adaptive strategies of compared N. galegae biovars.

3. Discussion

Root nodule bacteria (rhizobia) represent the genetically best studied group of symbiotic microorganisms. Being highly effective producers of N compounds for terrestrial ecosystems, these bacteria are characterized by exclusive ecological and agronomic importance. Moreover, rhizobia are used as a model to address a range of general evolutionary problems including the genomic mechanisms of micro- and macroevolution [8], and the tradeoff between genetic polymorphism and natural selection [10]. As it was demonstrated previously, purifying (stabilizing) selection usually results in a decreased population/gene polymorphism [11] while disruptive and negative frequency-dependent selection result in an extended (balanced) polymorphism [12,13]. The data on influence of driving (directed) selection on the genetic polymorphism are contradictory: it may be increased [14,15,16,17], conserved or even decreased [18,19,20] by this selection.

Application of rhizobia as a model to study impacts of natural selection on genetic diversity is based on the complicated genomic structures in which the core parts encoding for housekeeping functions differ in their natural histories from the accessory parts including the symbiotically specialized (sym) genes. Convenient models for analyzing the rhizobia genome dynamics are represented by polytypic species Neorhizobium galegae composed of two host-specific biovars – bv. orientalis and bv. officinalis (symbionts of Galega orientalis and G. officinalis), and Rhizobium leguminosarum, composed of bv. viciae (symbionts of plants from tribe Fabeae, genera Lathyrus, Lens, Pisum, Vavilovia, Vicia), bv. trifolii (symbionts of genus Trifolium from tribe Trifolieae), bv. phaseoli (symbionts of Phaseolus vulgaris from tribe Phaseoleae).

Previously we demonstrated that in R. leguminosarum, sym and core genes differ greatly in nucleotide polymorphism (p-distance) and are phylogenetically non-congruent, suggesting that evolution of these genes is independent [3]. This independency may be resulted from intensive horizontal gene transfer in rhizobia populations wherein two gene groups recombine randomly due to location of sym genes on mobile plasmids [1].

In order to reveal impacts of natural selection on the rhizobia gene polymorphism we used the set of N. galegae strains originated from North Caucasus region. As we demonstrated earlier [7], diversity of nucleotide sequences (measured as p-distance) in N. galegae is higher for core genes than for sym genes and is biovar-dependent: bv. orientalis is more polymorphic than bv. officinalis for both gene groups (Table 1). This difference may be due to contrasting levels of diversity in the respective host plant species. Specifically, North Caucasus is the longstanding center of origin for G. orientalis while colonization of this area by G. officinalis is more recent [21]. Previously we quantified diversity of two Galega species in North Caucasus using the nucleotide polymorphism analysis in a range of genes followed by genomic fingerprinting and confirmed the morphological data suggesting a higher G. orientalis diversity with respect to G. officinalis [22,23]. An important source of genetic diversity in N. galegae may be represented by translocations of the Insertion Sequences (IS) which are more abandoned in bv. orientalis than in bv. officinalis [24].

In this paper, we demonstrated that in N. galegae, core and sym genes are phylogenetically congruent (Figure 1) apparently due to their restricted recombination based on location of sym genes on non-mobile chromids. Nevertheless, some evolutionary important parameters of diversity are different in these genes: tradeoff between nucleotide polymorphism and evolutionary impacts of natural selection depend on the gene group (core or sym) and on the N. galegae biovar (orientalis or officinalis) (Table 1). Analyses of the total gene pools (Table 2) as well as of Gene Ontology Groups (GOGs) (Table 3, Figure 3), suggest that driving selection pressures result in an increased polymorphism of core genes in bv. orientalis, not in bv. officinalis. We suggest that in bv. orientalis, maintenance of novel core gene alleles by driving selection (dN/dS > 1) may be combined with preservation of preexisted alleles due to a broad ecological amplitude of this biovar, therefore, its genetic polymorphism is elevated. However, in bv. officinalis, the newly emerged gene alleles possibly substitute the preexisted ones (due to a restricted ecological amplitude of this biovar), therefore, gene polymorphism in this biovar is not changed.

In accordance to contrasting ecological affinities of the Galega species, a range of differences between their symbionts were revealed: (i) low polymorphic GOGs are affiliated with N metabolism (apparently responsible for symbiotic adaptations) in bv. officinalis and with the synthesis of surface polysaccharides (responsible for adaptations to edaphic stresses) in bv. orientalis (Tables S1–S3 in Supplement); (ii) sym genes evolve under purifying selection (dN/dS < 1) impacts in bv. officinalis while a neutral evolution (dN/dS ≈ 1) was revealed for these genes in bv. orientalis; (iii) evolution of core genes occurs mostly under impacts of driving selection in bv. orientalis while this evolution is neutral in bv. officinalis (Table 1, Table 2 and Table 3, Figure 3). From these data, we can suppose that operation of sym genes is most critical for bv. officinalis because at North Caucasian region this biovar persists under unfavorable soil conditions and should survive mostly due to colonization of endosymbiotic niches. However, bv. orientalis persists under more favorable conditions, as compared to bv. officinalis, therefore adaptations to edaphic factors dependent on core genes are highly important for bv. orientalis.

Interestingly, N. galegae differs in its evolutionary dynamics from the previously studied [3] R. leguminosarum species: sym and core genes in N. galegae are more similar in their diversity parameters than in R. leguminosarum. This difference between two polytypic species may be due to a more restricted recombination of sym and core genes in N. galegae with respect to R. leguminosarum (Table 4). Comparative analysis of these species contributes sufficiently to understanding of tradeoff between microevolution, speciation and macroevolution and between genetic polymorphism and natural selection.

Specifically, we demonstrate that in both rhizobia species, N. galegae and R. leguminosarum core genes responsible for speciation and macroevolution differ greatly in their natural histories from sym genes responsible for microevolution. Different genetic mechanisms were proposed for micro- and macro-evolution by J. Philiptschenko (1927) [25] who was the first to define these processes and correlated them with changes of the eukaryotic nuclear and cytoplasmic genes, respectively. Later, R. Goldschmidt (1949) [26,27] suggested that microevolution is based on the minor adaptive changes (“micro-mutations”) which can not initiate the speciation and macroevolutionary processes dependent on “macro-mutations” (responsible for generation of “hopeful monsters”). The similar approach was proposed in the punctuated evolution concept [28,29], which may be apparently used to address the symbiosis evolution since hosting of symbiotic microbes by eukaryotic organisms represents the rapid evolutionary bursts in contrast to gradual evolution suggested by the models gradualist evolution based on natural selection [2,30].

The other important issue of the rhizobia evolutionary genetics pertains the tradeoff between driving selection and gene polymorphism which may be increased by this selection in an ecologically versatile organism (such as N. galegae bv. orientalis) allowing a broad allelic diversity in the analyzed genes. However, in an ecologically restricted organism (such as N. galegae bv. officinalis), gene polymorphism is not changed or is even decreased under impacts of driving selection since co-existence of different gene alleles is presumable blocked. An extended bioinformatics analysis is required to analyze a relationship between adaptive potentials of organisms and impacts of natural selection on their polymorphism expressed in the diverse rhizobia species and in other symbiotic organisms.

4. Materials and Methods

4.1. Collection of strains and DNA sequencing

During expedition to the North Caucasus in 2003, a number of soil samples was collected, from which rhizobia strains of bv. orientalis and bv. officinalis were isolated [22,23]. A total of 14 rhizobia strains were isolated from soil samples in a microvegetation experiment using nodules of Galega orientalis and G. officinalis according to standard protocol [31]. They include strains of bv. officinalis: NG_35_off (JANFGK000000000), NG_37_off (VYYB00000000), NG_46_off (JANFGL000000000), NG_47_off (JAMQCN000000000), NG_58_off (JANFGM000000000), NG_77_off (JANFGN000000000), NG_81_off (JANFGO000000000), NG_110_off (VZUM00000000)) and of bv. orientalis: NG_35_ori (JANFGP000000000), NG_46_ori (JANFGQ000000000), NG_58_ori (VZUN00000000), NG_77_ori (JANFGR000000000), NG_87_ori (VZUL00000000), and NG_110_ori (JANEZU000000000).

Isolates were cultivated at 28°C and 220 rpm for 48 h in modified yeast mannitol broth (YMB) with 1% sucrose [32]. DNA was obtained by the lysozyme–SDS–phenol–chloroform extraction protocol, with minor modifications [33]. Sequencing of strains NG_37_off and NG_87_ori was performed on a Pacbio RSII instrument with P6 in two SMRT cells (Pacific Biosciences of California, Inc., Menlo Park, CA, United States). PacBio sequencing and subsequent error correction analysis and assembly were performed at Arizona Genomics Institute (US). The assembly the strains was carried out de novo using HGAP (https://github.com/jtchien0925/PacBio_HGAP_assembly). Sequencing of other 12 rhizobia strains, 7 strains of biovar officinalis and 5 strains of biovar orientalis, was performed on a MiSeq genomic sequencer (Illumina Inc., San Diego, CA, United States), according to the manufacturer’s protocol, using the MiSeq Reagent Kit, 600 Cycles (Illumina, United States) at the Genomics Core Facility, Siberian Branch, Russian Academy of Sciences (Institute of Chemical Biology and Fundamental Medicine, Novosibirsk). The assembly of the sequences was carried out using the CLC Workbench (https://digitalinsights.qiagen.com/products-overview/discovery-insights-portfolio/analysis-and-visualization/qiagen-clc-genomics-workbench/) by mapping on reference genomes NG_37_off and NG_87_ori each biovar respectively.

4.2. Finding core genes by global alignment

All genes of each genome were matched with the genes of other 13 genomes using the global alignment method. For this, BLAT was used (http://genome.ucsc.edu/cgi-bin/hgBlat). All paired genes were sorted in descending order of identity. The paired genes with a maximal identity of at least 70% in the DNA sequence were selected. After that a table of all genes and their presence in each of the 14 strains was generated. Only genes that were found in all 14 strains were selected for the core genome.

For bv. orientalis, the estimated number of core genes is approximately 5000, for bv. officinalis – 4200, while 3900 core genes are common for two biovars. For analyzing the variability and selection indices we used 3840 genes from bv. orientalis and 2734 genes from bv. officinalis (genes are common for two biovars with non-zero polymorphism and dN/dS).

4.3. Symbiotic genes

For both biovars we analyzed 16 nod genes (encoding for Nod factor synthesis) 8 nif genes (for nitrogenase synthesis) and 15 fix genes (for electron and energy supply of nitrogenase).

4.4. Gene alignment using Muscle

DNA sequences of 14 strain variants of each gene were aligned by MUSCLE (Multiple Sequence Comparison by Log-Expectation, https://www.ebi.ac.uk/Tools/msa/muscle/) using standard coding sequence alignment parameters.

4.5. Calculation of nucleotide polymorphism (p-distance)

DNA polymorphism of each gene was calculated based on the number of nucleotide substitutions for each pair of strains using standard metrics (https://www.megasoftware.net/mega1_manual/Distance.html). DNA regions with undetermined sequences (N, non-detected) and gaps were not taken into account. The number of substitutions was normalized by dividing the total length of the compared genes without gaps and undefined nucleotides. The matrixes (sized 14X14 by the number of strains) of p-distances of each gene were calculated. The average polymorphism of each gene was calculated using the average p-distance of all elements of the matrix excluding diagonal elements (distance of a gene with itself is zero).

4.6. Calculation of dN/dS index

The calculation of the dN/dS ratio of nonsynonymous (dN) to synonymous (dS) substitutions was performed according to the Jukes-Cantor (JC) model (https://bioinformatics.cvr.ac.uk/calculating-dnds-for-ngs-datasets/). In JC model, dN/dS for each codon was calculated separately and compared with the theorized ratio of substitutions. For example alanine is encoded by three different codons when there are nine possible single substitution of each codon and consequently its theoretical dN/dS ratio is 3/9 or 1/3. The obtained dN/dS value of each gene was normalized by dividing by the number of coding codons of the compared sequences. Then the matrixes (sized 14X14 by the number of strains) of dN/dS indexes of each gene were calculated. The average dN/dS index of each gene was calculated by the average dN/dS index of all elements of the matrix excluding diagonal elements.

4.7. Functional annotation of core genes, Gene Ontology (GO)

We used eggnog-mapper (https://github.com/eggnogdb/eggnog-mapper/issues/135) to annotate newly assembled genomes and assign genes to certain functional groups of Gene Ontology. A detailed transcript of each group was performed on the AmiGO 2 website [34]. AmiGO 2 is a project to create the next generation of AmiGO, the current official web-based toolkit for searching and browsing Gene Ontology’s database. The Gene Ontology Consortium (GOC) provides computable knowledge regarding the functions of genes and gene products.

4.8. Determination of the predominance of functional groups of genes GO (gene Ontology) and statistical significance

The prevalence of certain groups of genes was calculated as the ratio of the actual number of genes to their expected number based on the sample size and the total number of genes of the selected group,

P_enrich=N_obs/(S_sep*(N_sep/N_genome)),

where P_enrich is the predominance of a given group of genes, N_obs is the number of genes in the sample, S_sep is the number of genes in a given group (GO), N_sep is the sample size, N_genome is the total number of genes found in Gene Ontology.

The statistical significance of the predominance of certain groups of genes was obtained using a permutation test that simulates the same value of the size of the group and the total number of genes. The permutation test was performed 10,000 times, this was enough to calculate statistical significance 95% (P. val. <0.05).