Genome-wide functional divergence after the symbiosis of proteobacteria with insects unraveled through a novel computational approach

Abstract

Symbiosis has been among the most important evolutionary steps to generate biological complexity. The establishment of symbiosis required an intimate metabolic link between biological systems with different complexity levels. The strict endo-cellular symbiotic bacteria of insects are beautiful examples of the metabolic coupling between organisms belonging to different kingdoms, a eukaryote and a prokaryote. The host (eukaryote) provides the endosymbiont (prokaryote) with a stable cellular environment while the endosymbiont supplements the host's diet with essential metabolites. For such communication to take place, endosymbionts' genomes have suffered dramatic modifications and reconfigurations of proteins' functions. Two of the main modifications, loss of genes redundant for endosymbiotic bacteria or the host and bacterial genome streamlining, have been extensively studied. However, no studies have accounted for possible functional shifts in the endosymbiotic proteomes. Here, we develop a simple method to screen genomes for evidence of functional divergence between two species clusters, and we apply it to identify functional shifts in the endosymbiotic proteomes. Despite the strong effects of genetic drift in the endosymbiotic systems, we unexpectedly identified genes to be under stronger selective constraints in endosymbionts of aphids and ants than in their free-living bacterial relatives. These genes are directly involved in supplementing the host's diet with essential metabolites. A test of functional divergence supports a strong relationship between the endosymbiosis and the functional shifts of proteins involved in the metabolic communication with the insect host. The correlation between functional divergence in the endosymbiotic bacterium and the ecological requirements of the host uncovers their intimate biochemical and metabolic communication and provides insights on the role of symbiosis in generating species diversity.

Figure 1. Constraints operating in endosymbiotic bacteria of aphids (A) and carpenter ants (B) in comparison with their free-living relative bacteria Escherichia coli and Salmonella typhimurium.

We compared the constraints operating in protein-coding genes between endosymbiotic and free-living bacteria by dividing the non-synonymous-to-synonymous rates ratio of endosymbionts (ω_e) by that of their free-living relatives (ω_f) and we called this ratio R(ω) [R(ω) = ω_e/ω_f; represented in the Y-axis). We plotted genes according to their position in the bacterial chromosome (X-axis). We also indicate R(ω) = 1 since this is the value at which genes have not changed their selective constraints.

Figure 2. Distribution of highly constrained genes among the functional categories in Buchnera sp. (blue bars) and Blochmannia sp. (red bars).

The different functional categories as explained by the Cluster of Orthologous Groups (COG) are represented in the X-axis. The height of the bar represents the relative contribution of each class (i) of size (t), to the total number of genes under strong selective constraints (n_i: R(ω) = ω_e/ω_f<1) when considering the whole dataset (T). This normalized number hence was calculated as φ = (n_i/t) * (t/T). Classes showing significant enrichment with highly constrained genes under a hypergeometric distribution are labeled by (*, P<0.05; **, P<10⁻²; ***, P<10⁻³). We also labeled those functional classes significantly underrepresented by highly constrained genes using green stars.

Figure 3. Distribution of genes under functional divergence among the functional categories in Buchnera sp. (blue bars) and Blochmannia sp. (red bars).

The different functional categories as explained by the Cluster of Orthologous Groups (COG) are represented in the X-axis. The height of the bar represents the relative contribution of each class (i) of size (t), to the total number of genes under functional divergence (n_i: R(ω) = ω_e/ω_f<1) when considering the whole dataset (T). This normalized number hence was calculated as φ = (n_i/t) * (t/T). Classes showing significant enrichment with genes under functional divergence under a hypergeometric distribution are labeled by (*, P<0.05; **, P<10⁻²; ***, P<10⁻³). We also labeled those functional classes significantly underrepresented by highly constrained genes using green stars.

Figure 4. Distribution of genes under functional divergence among the metabolic pathways significantly enriched or impoverished with these genes in Buchnera (A) and Blochmannia (B).

The different metabolic classes are color-coded. Dotted line separates metabolic pathways enriched with functionally divergent genes (above the line) from those impoverished with these genes (below the line). The height of the bar represents the relative contribution of each class (i) of size (t), to the total number of genes under functional divergence (n_i: R(ω) = ω_e/ω_f<1) when considering the whole dataset (T). This normalized number hence was calculated as φ = (n_i/t) * (t/T).

Figure 5. Genome wide identification of functional divergence.

Proteins are identified to be under functional divergence if they show amino acid sites presenting significant evidence of shifts in the evolutionary rates in cluster 1 (cluster under study) compared to cluster 2 (background cluster). To measure functional divergence at site i, we first calculate all pair-wise BLOSUM transition values in the pair-wise comparison of the sequences in the tree. Sequences in cluster 1 are compared to the outgroup and the BLOSUM transition values between outgroup and cluster 1 generate a distribution that is compared to that generated when comparing sequences of cluster 2 to the outgroup. The change in the physico-chemical properties of amino acids from cluster 1 to cluster is indicated by colored squares. If the transition scores are significantly more radical when comparing the outgroup to cluster 1 at that amino acid site of the protein than when we compared the outgroup to cluster 2 then we consider the site to be under functional divergence. The significance of the transition scores in cluster 1 is calculated by comparing the distribution of scores in cluster 1 to that in cluster 2 and significance is considered at the 1% confidence level.