A comparative genomic analysis of gene position conservation

In order to infer chromosome topologies and their association with gene dosage gradients down the oriC−ter axis, we established a comparative genomic framework. Our approach is based on the frequencies of translocations between distal segments of the chromosome, as defined by sequence comparisons between pairs of closely related organisms. This is based on the supposition that spatially proximal regions of the chromosome would recombine more frequently, and recombination combined with selection would result in the establishment of translocation events between such loci. This approach is similar to the experimental mapping of the topology of the E. coli chromosome using recombination frequencies between distal loci (Valens et al. 2004), with the difference that the translocations we identify are products of recombination followed by selection and probably drift over long timescales.

We obtained complete genome sequences and gene annotations for ∼3000 fully sequenced bacterial genomes from the NCBI. Because gene dosage gradients down the oriC−ter axis are dependent on growth rates, we sought to annotate each organism in our collection of genomes with growth rate information. To do so, we assembled a dataset of growth rates in the form of minimum doubling times for over 100 bacteria from Freilich et al. (2009).

Although growth rate information is sparse relative to genome sequence data, we know that the number of rRNA genes encoded by a bacterium is a predictor of its growth rate (Condon et al. 1995; Klappenbach et al. 2000), with fast-growing bacteria encoding more rRNA operons than slow-growing ones. Consistent with this knowledge, we observed a negative correlation (ρ_Spearman= −0.76) between minimum doubling times and rRNA copy numbers for these ∼100 bacteria (Figure 1A).

• Download figure
• Open in new tab
• Download powerpoint

Figure 1

Comparative Genomics Framework. (A) Plot representing the relationship between minimum doubling time (hr) and 16S rRNA copy number (ρ_Spearman = −0.76, P-value <10⁻¹⁰). Minimum doubling time values >50 hr are not represented in the plot. (B) Flowchart showing the selection of 262 pairs of bacteria used in this study. (C) Schematic showing the conservation of gene position (solid lines) and gene translocations (dashed lines) between pairs of bacteria. (D) Plot representing the correlation between residual ρ_GOC and R_f (ρ_Spearman = 0.31, P-value = 2.2 × 10⁻⁷). Inner panel shows a statistically significant difference between residual ρ_GOC for slow (R_f ≤1) and fast (R_f > 1) growing bacteria (P-value = 3.3 × 10⁻⁶, Wilcoxon test). Asterisks indicate P-value < 10⁻³.

Based on the fit between minimum doubling times and rRNA copy numbers, we estimated the growth rate expected for each rRNA copy number observed in our dataset. However, growth rate itself is not necessarily a strong predictor of differences in gene dosage between oriC-proximal and ter-proximal loci. Instead, a better correlate is the R-factor (R_f), defined as the number of replication cycles per cell division. Using a previous study from the Rocha group as a guide (Couturier and Rocha 2006), we calculated the R_f for each bacterium as described in Materials and Methods. We then defined organisms that are expected to initiate replication more than once per cell cycle on average (R_f > 1) as fast-growing, and the rest as slow-growing. We also estimated R-factor independently using a larger, more recent dataset (Vieira-Silva and Rocha 2010) of minimum doubling times for ∼200 bacteria, and found estimates of R-factor derived from this dataset to be consistent with those used in this study (Figure S1).

We organized the ∼3000 genomes into 262 genome pairs on the basis of phylogenetic distance and growth rate. We considered pairs of bacteria with ≥97% 16S rRNA sequence identity (Figure S2) and selected those with similar growth rates (same 16S rRNA copy numbers) (Figure 1B). We then removed: (a) pairs of strains of the same species, which typically showed very few long-range translocations (Figure S3); (b) redundant pairs to minimize overrepresentation of certain phylogenetic groups: for example, there are several hundred comparisons involving E. coli and Salmonella, and of these one was picked at random; and (c) pairs of obligate parasites with genomes smaller than 1.3Mb (Koonin and Wolf 2008), whose reduced genomes are known to be stable with exceptionally high gene order conservation (Tamas 2002; Mira et al. 2002). Thus, our final selection consisted of 232 individual bacterial species forming 262 pairs of bacteria with the same 16S rRNA copy number and ≥97% 16S rRNA sequence identity. A phylogenetic tree of these bacteria, based on their 16S rRNA sequences, is shown in Figure S4.

For all the 262 pairs of bacteria, we identified gene orthologs by bidirectional best-hit phmmer (E-value <10⁻¹⁰). We found that ∼72% (median) of genes in one organism had an ortholog in its partner genome. The percentage of genes conserved in a pair of organisms decreased with phylogenetic distance between the compared genomes and showed little correlation with growth rate (Figure S5).

For all these genes, we computed the position of each gene on its chromosome. The position of a gene was defined by its distance from oriC along the path of the replication fork and normalized to genome size. Positions equidistant from oriC on either replichore were considered equivalent in these calculations. Thus, a difference in gene position would indicate a difference in dosage, depending on the organism’s R-factor. For each genome and its partner, we defined gene position vectors, in which each gene with an identifiable ortholog in its partner genome was represented by its distance from oriC. Then, for each genome pair, we computed the correlation coefficient between their gene position vectors (Figure 1C) and treated this correlation as an estimate of gene order conservation (ρ_GOC).

As expected for closely related organisms, the average gene order conservation across the 262 pairs of organisms studied here was high (median ρ_GOC ∼0.9). Loss of gene order conservation is dependent on phylogenetic distance (Suyama and Bork 2001; Tamames 2001). Hence, we calculated residual ρ_GOC after correcting for the dependence of ρ_GOC on phylogenetic distance (Figure S6). We found that the residual ρ_GOC increases with increasing R-factor (R_f) weakly but in a statistically significant fashion (P-value= 3.3 × 10⁻⁶, Wilcoxon test; Figure 1D), consistent with earlier calculations by Rocha and colleagues (Couturier and Rocha 2006).

To identify patterns of long-range translocations, which could inform us about chromosome topologies, we took the following coarse-grained approach. We binned each chromosome into four equally sized segments. The first segment, termed ’O’ for the presence of oriC, comprised 25% of the chromosome centered around oriC. The ’L’ (Left) and the ’R’ (Right) segments each covered 25% of the chromosome on either side of the O segment. The remaining 25%, including the terminus of replication, was termed the ’T’ (Terminus) segment (Figure 2A). Each gene with an ortholog in a partner genome was assigned to one of these four bins.

• Download figure
• Open in new tab
• Download powerpoint

Figure 2

Chromosomal bin-based conservation of function. (A) Schematic representing equal sized bacterial chromosomal segments. Red: Origin bin (O); green: Right bin (R); blue: Terminus bin (T); black: Left bin (L). (B) Boxplot showing the enrichment of translation and ribosome biogenesis genes (COG category J) in O, R, T, and L chromosomal bins. (The dashed gray line at y = 1.3 corresponds to a P-value of 0.05. Values ≥1.3 are statistically significant.). Inner panel shows the enrichment of COG J in the O bin of fast-growing bacteria (R_f > 1). (C) Boxplot showing the depletion of horizontally acquired genes in different chromosomal bins. Inner panel showing the depletion of horizontally acquired genes in the O bin of both slow- (R_f ≤ 1) and fast (R_f > 1) growing bacteria.

Analysis of gene functions, for each bin, at a broad level using the COG annotations revealed that the only statistically significant association between gene function and gene positioning was the enrichment of genes involved in translation and ribosome biogenesis in the O (origin) bin (Figure 2B and Figure S7). In line with previous literature (Couturier and Rocha 2006), we observed that the average distance of translation and ribosome biogenesis genes from oriC decreases with increase in R_f (Figure S8). The enrichment of ribosome biogenesis and translation genes in the O segment was particularly apparent in fast-growing bacteria (R_f > 1) but not in slower growing ones (Figure 2B). In contrast, there was a significant depletion of predicted horizontally acquired genes in the O bin, in both fast- and slow-growing bacteria, and more so in the former (Figure 2C). This is consistent with the increased chance of novel gene integration disrupting essential gene functions in this region of the chromosome. These results are also in agreement with observations made in E. coli (Zarei et al. 2013) wherein horizontally acquired genes were defined using methods different from what we have used here.

Long-range interreplichore translocations

The degree of conservation – defined as the proportion of genes with an ortholog in a closely related partner genome – decreased from the O bin toward the T bin for fast-growing organisms, with slow-growing organisms showing little effect (Figure S9).

We then asked whether genes present in a given segment of a chromosome would be conserved within the same segment in its closely related partner genome and whether this property would differ across the four segments defined here. We define P_i – the gene position conservation for segment i (where i stands for the four chromosomal segments O, T, R, or L) − as the proportion of genes in that segment, which are conserved in the corresponding segment in its partner genome (Figure 3A). The median P_O and P_T, for the O and the T bins, across our dataset of 262 genome pairs, was ∼80%, and the distributions of P_O and P_T were tight. In contrast, P_L and P_R showed a wide distribution, with the median around 50% (Figure 3A).

• Download figure
• Open in new tab
• Download powerpoint

Figure 3

Dynamics of gene organization. (A) Schematic and boxplot showing the conservation of gene position (P_i) in O bin (O−O), T bin (T−T), R bin (R−R), and L bin (L−L). (B) Schematic representing O−T and T−O translocations (black arrows) and L−R and R−L translocations (gray arrows) and boxplot showing the proportion of these translocations (T_(i,j)).

Next, we asked whether a gene present in a particular bin in one organism translocates to the diametrically opposite bin (or a neighboring bin − Figure S10) in the genome of a closely related organism. We called this translocation rate (T_i-j, where i and j indicate the segments between which the translocation rate is measured) (Figure 3B). Consistent with our observation of high P_O and P_T, we found that T_O−T and T_T−O are very low with a median of ∼1% (Figure 3B). This could be a possible consequence of selection acting against translocations that would drastically disrupt gene dosage, and is consistent with the observation of few contacts between oriC and ter.

We next looked at translocations between the L and the R segments (T_L−R and T_R−L), which represent interreplichore translocations (Figure 3B). For ∼42% of genome pairs, T_L−R/R−L <0.1; for the remaining, T_L−R/R−L ∼0.45 (median). Therefore, in the majority of genome pairs compared here, genes present in one replichore in one genome are often present in the opposite replichore in its partner genome.

Note that in identifying interreplichore translocations, by considering only the L and the R bins, we neglected 25% of each replichore immediately adjacent to oriC and ter. This was done to ensure that our observations were not skewed by previously reported high inversion rates immediately around oriC and ter (Eisen et al. 2000; Suyama and Bork 2001; Tillier and Collins 2000). Note, however, that the proportion of interreplichore translocations within the O and the T segments correlate well with that between the R and the L segments. Translocations appear to be more frequent in the ter-proximal than in the oriC-proximal half of the chromosome (Figure S11). The difference in translocation rates between the oriC-proximal and the ter-proximal halves of a chromosome not only suggests a complementary strategy for studying long-range translocations, but also further corroborates the interplay of chromosome topology and gene organization with the process of replication.

To ensure that the interreplichore translocations we describe are not merely a result of discrepancies in gene annotations, we predicted genes using Glimmer (Delcher et al. 1999) for all genomes included in this study. We compared predicted gene start positions with those reported in NCBI and found that the two datasets are in agreement with each other (Figure S12).

In summary, we show that interreplichore translocations between the L and the R bins are probably more common than could be perceived from the literature (Eisen et al. 2000; Tillier and Collins 2000). Further, the variation in the distribution of L−R translocation rates is only weakly explained by phylogenetic distance (Figure S13 and Figure S14).

Interreplichore translocations as a possible readout of chromosome topology

Rates of recombination have been used as a proxy for measuring contact frequencies between distant regions of a chromosome. Together with the higher resolution analysis of contact frequencies by chromosome conformation capture methods and microscopy, these have so far underscored two broad chromosome topologies. First, in E. coli it has been shown that oriC is localized to the midcell and the left and the right arms of the chromosome stretch toward opposite poles with few interreplichore contacts (Nielsen et al. 2006; Wang et al. 2006). In line with this, we observed low interreplichore translocations (T_L−R/R−L ∼0.03) between E. coli and Salmonella (two closely related bacteria).

In contrast, in C. crescentus (Le et al. 2013; Umbarger et al. 2011) and M. smegmatis (Santi and McKinney 2015), oriC and ter are present at the two poles of the cell and the left and the right arms of the chromosome are in close contact with each other. Both these genomes showed high interreplichore translocations (T_L−R/R−L ∼0.50 and ∼0.70, respectively), consistent with their chromosome topology.

Next, we analyzed the chromosome contact map from a previously published Hi-C experiment for C. crescentus (Le et al. 2013). Pairs of loci translocated between the two replichores displayed significantly higher contact frequencies in Hi-C experiments when compared to pairs of interreplichore loci not showing translocations (Figure 4A and Figure S15). These indicate that high rates of translocation are localized to spatially proximal regions of the chromosome.

• Download figure
• Open in new tab
• Download powerpoint

Figure 4

Interreplichore translocations and chromosomal contact points. (A) Boxplot representing the normalized interreplichore contact frequencies as derived from Le et al. (2013) for C. crescentus NA1000 (NC_011916); for all the probed gene positions with or without interreplichore translocations (P-value < 10⁻¹⁰, Wilcoxon test). (B) Scatterplot representing 386 interreplichore (R → L) translocations between C. crescentus NA1000 (NC_011916) and C. segnis (NC_014100) in red and randomized dataset (as described in Materials and Methods) in blue. Inner panel shows the distribution of D_inter = | d_Ri/Lj – d_Li/Rj | (solid line) for interreplichore translocations between C. crescentus NA1000 (NC_011916) and C. segnis (NC_014100). Dashed line indicates distribution of D_inter for randomized data (as described in Materials and Methods).

Many interreplichore translocations in C. crescentus (Figure 4B) and in M. smegmatis (Figure S16A) are symmetric, i.e., a gene at a certain distance d_R from oriC in one replichore translocates to a position d_L in the other replichore such that D_inter = |d_R − d_L| is close to zero. This is again consistent with the topology of these chromosomes, and also with the idea that a pair of concurrent replication forks, which are typically equidistant from oriC, are hotspots for recombination. Such symmetric translocations are true, not only for C. crescentus and M. smegmatis, but are a general feature of genomes with high T_L−R/R−L values (Figure 5A), but less so for genome pairs with low T_L−R/R−L (T_L−R/R−L < 0.1) (Figure S16B). These are in line with the high ρ_GOC values across most genome pairs tested here, indicating that interreplichore translocations rarely disrupt the relative distance of genes from oriC, and thereby maintain their dosage.

• Download figure
• Open in new tab
• Download powerpoint

Figure 5

Selective value of chromosome topology. (A) Plot representing the distribution of D_inter = |d_R/L – d_L/R| for all bacteria showing high interreplichore translocations (T_L−R/R−L > 0.1). Dashed line indicates the distribution of D_inter for randomized datasets (as described in Materials and Methods), similar to the inner panel of Figure 4B. Inner panel shows the percentage of gene pairs with D_inter < 0.1 for all bacteria showing high interreplichore translocations, alongside random expectation (determined from randomized datasets as described in Materials and Methods). (B) Plot representing the proportion of leading-lagging strand flips post interreplichore translocations as a function of the proportion of interreplichore translocations (ρ_Spearman = 0.17, P-value < 10⁻¹⁰) in cyan. A randomized dataset (as described in Materials and Methods) is shown in gray. Adjacent bloxplot represents the distribution of leading-lagging strand flips in bacteria with high and low interreplichore translocations alongside random expectation (R_m) for each dataset (D) (low T_L−R/R−L, P-value = 2.9 × 10⁻⁶, Wilcoxon test and high T_L−R/R−L, P-value < 10⁻¹⁰,Wilcoxon test). Asterisks indicate P-value < 10⁻³.

We also found that interreplichore translocations do not impact strand bias in gene densities. We noticed that the translocated partners of leading strand genes tend to be encoded on the leading strand of the opposite replichore. This is more prominent for organisms showing high T_L−R/R−L as compared to low T_L−R/R−L (Figure 5B). This difference in strand maintenance between genomes with high and low T_L−R/R−L values might point to mechanistic differences between the two groups of organisms in the way these translocations occur.

Finally, we observed that interreplichore translocations rarely break operons. However, the contiguity of genes belonging to the same functional category is often broken by these translocations (Figure S17). In order to explicitly test the impact of these interreplichore translocations on gene expression states, we obtained publicly available microarray data for Shewanella oneidensis MR1 (high T_L−R/R−L; T_L−R/R−L ∼0.56) across many conditions from the M3D database (http://m3d.mssm.edu). Analysis of these microarray data showed that genes present on different replichores show a higher correlation in expression (ρ_Spearman > 0.5) when they are equidistant from oriC. In contrast, gene pairs showing little or negative correlations between their gene expression profiles tend to show greater differences in their distance from oriC. Although the effect is very weak, the trend is statistically significant in terms of the difference between medians (P-value < 10⁻¹⁰, Wilcoxon test; Figure S18). This suggests that interreplichore translocations between genes equidistant from oriC may not disrupt their gene expression states.

Thus, interreplichore translocations − potentially encouraged by the spatial proximity of the two replichores − neither disrupt gene dosage gradients down the oriC−ter axis nor do they impact strand biases in gene densities (Figure S22); they may have only a small effect on gene expression states.