The comparison of the gene expression profiles in exponential and stationary phase reveals that the oriC-proximal genes have a higher expression in exponential than in stationary phase, while the terminus-proximal genes have a lower expression level (Sobetzko et al., 2013). A plot representation of the 2013 E. coli CSH50 wild type data is shown in Figure 1A. This spatio-temporal gene expression pattern may reflect a cellular program to adapt to changing conditions. The pattern may emerge due to the strategic positioning of genes regulated by global transcription factors such as σ⁷⁰, σ³⁸ or abundant regulatory proteins like Fis, H-NS, IHF, the cAMP receptor protein (CRP), or the leucine-responsive regulatory protein (Lrp). These factors regulate hundreds of genes and therefore may have an impact on a global pattern. A new analysis of the most recent regulonDB data set of target gene positioning of global regulators confirms a gradient of regulated genes in the σ⁷⁰, σ³⁸, CRP, and Lrp regulons along the oriC-ter axis (see Figures 1B–D and Supplementary Figure S6) found in Sobetzko et al. (2012). Activity of these factors depends on the cellular state. While σ⁷⁰ is the dominant transcription factor during exponential phase, σ³⁸ competes with σ⁷⁰ for RNA polymerase (RNAP) during stationary phase (Jishage et al., 2002). The updated regulatory data of the regulonDB database confirmed genes specifically regulated by σ⁷⁰ are more abundant in proximity to oriC (Sobetzko et al., 2012) (see Figure 1B), which may contribute to the observed increase of oriC-proximal genes during exponential phase. Genes regulated specifically by σ³⁸, however, are more abundant at the terminus region (see Figure 1C), which can explain the observed increase of oriC-distal gene expression in stationary phase compared to exponential phase. CRP and Lrp, both important regulators during starvation and stationary phase (Schultz et al., 1988; Tani et al., 2002), regulate genes with a characteristic distribution gradient along the oriC-ter axis (see Figure 1D). In this case, repression of more oriC-proximal genes during stationary phase would also contribute to the observed increase in gene expression of oriC-proximal genes during exponential phase. In general, a combination of different regulatory proteins with non-uniform distribution of target genes along the oriC-ter axis contributes to the spatial gene expression pattern observed when comparing exponential and stationary phase. These patterns can be supported by DNA supercoiling sensitivity of promoters mediated by DNA structure and regulatory proteins (Blot et al., 2006; Sobetzko et al., 2012; Lal et al., 2016; Pineau et al., 2022).

Besides promoter regulation, which differs in its activity regarding growth phases, replication activity is another potential regulatory factor. During exponential growth, many bacteria perform multifork replication to ensure chromosome replication when the doubling time is shorter than the duration of replication (see Figure 1E). This resulting copy number effect leads to a higher expression of oriC-proximal genes during fast growth compared to slow growth or stationary phase. Toward the terminus region, this effect is gradually reduced. Marker-Frequency-Analysis (MFA) allows to visualize and quantify the copy number effect when using whole-genome DNA sequencing of exponential growing cells (see Figure 1F). For E. coli, a gradual decrease of reads along the oriC-ter axis is observed representing the average copy number of the sequenced culture (Skovgaard et al., 2011). However, both regulatory factors, promoter regulation and copy number effect, act in the same direction regarding increasing and decreasing oriC-proximal/distal genes during exponential and stationary phase. Therefore, it is only possible to determine the impact of each factor on gene expression by isolating a single factor.

To determine the influence of promoter regulation and the copy number effect on the gene expression pattern of exponential phase compared to stationary phase, one factor has to be altered and the other has to be kept constant. As changes in the genetic regulatory network would be difficult due to the diversity of regulatory proteins and DNA topology (Sobetzko, 2016; Martis B et al., 2019; Klein et al., 2021), we decided to significantly alter the copy number effect. By relocation of the oriC into the terminus region, one would expect an opposite copy number profile compared to wild type. In exponential phase, genes initially located in the ter domain in wild type would have a higher copy number than genes initially located in the oriC-domain in wild type. This inversion of gene copy number could result in an unaltered, a disturbed or an inverted expression profile when compared to stationary phase depending on the impact of each of the regulatory factors. However, relocating the oriC into the terminus can cause massive biological problems. In the terminus region, the replication forks are trapped at specific DNA sites called ter sites, which are bound by the terminus utilization substance protein (Tus) (Hill, 1992). This protein-DNA complex unidirectionally arrests DNA replication, allowing replication forks to pass ter sites only in the origin-to-terminus direction. An oriC in the terminus region would, therefore, lead to replication fork stalling shortly after initiation and prevent replication of the remaining chromosome. To circumvent this problem, we generated an E. coli MG1655 Δtus strain to abolish replication stalling at ter sites. This allowed the replication forks to pass freely from the former wild type terminus to the former oriC region.

Another problem is head-on replication-transcription conflicts of the highly transcribed ribosomal RNA (rrn) operons and the replication forks, as the rrn operons are transcribed in oriC-ter direction. These head-on collisions seem to significantly delay fork progression and especially the rrnCABE cluster and the rrnH operon causes substantial problems to replication progression (Ivanova et al., 2015). Therefore, we needed to alter the transcription direction of the rrn operons. The Cre-lox and FLP/FRT systems, which are based on site-specific recombinases, allow excision and inversions of chromosomal DNA flanked by two identical target sites depending on its relative orientation. By flanking rrn operons with facing FRT or loxP site pairs, it is possible to invert the transcription direction and circumvent head-on replication-transcription machinery collisions, when relocating the oriC to the terminus region. To minimize crosstalk between FRT/loxP sites of different inversion sites, different FRT/loxP variants were used for each pair (Missirlis et al., 2006; Turan et al., 2010). As the rrnCABE cluster consists of four nearby ribosomal operons in oriC-proximity, we only used one pair of FRT sites to invert the whole region instead of inverting every single operon on its own. All insertions and deletions were carried out using the CRISPR SWAPnDROP system, which allows consecutive chromosomal changes based on CRISPR/Cas9 counter-selection (Teufel et al., 2022).

For the relocation of the oriC into the terminus, we first replaced the native oriC with the F-plasmid origin of replication oriS flanked by a pair of tandem FRT sites to allow excision of the oriS at a later stage using the FLP/FRT system. This also allows for a parallel inversion of the ribosomal RNA operons together with the oriS deletion in a single step to avoid head-on collisions in intermediate strains in case of a sequential genome editing approach. After the replacement of oriC with oriS, we inserted the native oriC into the terminus region of the chromosome. The strain was then transformed with a plasmid harboring the Cre recombinase and Flippase under the control of the pBAD promoter. Additionally, we used CRISPR/Cas9 to actively remove oriS DNA looped-out during excision by FLP/FRT-recombination to prevent reintegration of generated mini replicon. In summary, the final strain was able to freely invert its rrn operons and remove the placeholder oriS upon induction of the CRISPR/Cas9 and recombination systems to generate a strain with an inverted copy number gradient.

After induction, streaked colonies appeared in different sizes ranging from very small to wild type-like size. OriS elimination could only be found in small and middle-size colonies. Surprisingly, rrn operon inversions occurred rarely and could not be found in combination with the oriS elimination. Colonies of different sizes and with oriS knockout were re-streaked for further investigation. Re-streak of the small-size colonies resulted in a mix of small and middle-sized colonies indicating instability of the strain due to frequent suppressor mutation. For stability reasons, we used one of the middle-sized colonies that originated from a re-streaked small-size colony for further investigations. These colonies remained stable and homogeneous with respect to colony size and morphology for several restreaks. MFA-analysis of the clone during exponential phase revealed an inversion spanning half of the chromosome, mainly the left replichore (see Figures 2A–D). This inversion resulted in a relocation of the oriC from the terminus back into the wild type oriC region with a distance of about 381 kb from the native oriC site. Thereby, for most right Replichore genes, a wild type oriC distance configuration was restored whereas most genes of the left replichore remained inverted with respect to oriC distance. Furthermore, for the rrnCABE cluster and the rrnH operon the wild type leading strand arrangement was restored and therefore head-on collision with replication was prevented. This might have improved fitness and explains the increased colony size (Ivanova et al., 2015). Additionally, this strain revealed a significantly reduced maximal copy number at oriC during exponential phase compared to wild type (see Figures 1F, 2D). Consistently, the doubling time of this inversion strain (INV) at around 62 minutes is three times greater compared to the wild type (see Supplementary Figure S1A).

With its stability and the strongly reduced copy-number in exponential phase, this strain can be used for further investigation of copy number impact on the spatio-temporal expression pattern. If copy numbers have a strong impact, comparison of the exponential and stationary phase should show a reduced or abolished expression gradient along the oriC-ter axis. For the reference strain, we decided not to take the wild type E. coli MG1655 parental strain but a strain containing all the modifications leading to the INV strain to avoid potential effects of the modifications when comparing to the INV strain. The strain prior to induction of recombinases and CRISPR/Cas9 was selected. It shared all modifications except for the inversion and the presence of oriC in the terminus region as well as the oriS in the native oriC location. In order to obtain the reference strain, the terminus oriC was removed and the place-holder oriS at the native oriC site was replaced by oriC to recover the wild-type origin of replication. This reference strain finally resembled the wild type genome arrangement and at the same time shared the modifications required to generate the INV strain (see Supplementary Figures S1A, B). The fitness and integrity of the reference strain was then tested and compared to the E. coli MG1655 parental strain. Concerning fitness and growth speed (19.6 min doubling time) , the reference strain and E. coli MG1655 showed no differences. MFA-analysis showed a replication profile similar to E. coli MG1655 regarding spatial copy number distribution (see Figures 1F, 2E). Comparison of RNA-seq data of E. coli MG1655 and the reference strain showed a similarity in gene expression that equals the similarity of replicates of a single strain (see Supplementary Figures S7C, D). Hence, from the perspective of mRNA levels, the reference strain cannot be distinguished from its wild type precursor E. coli MG1655. In analogy to the RNA-seq data analysis of E. coli CSH50 wild type (Sobetzko et al., 2013), a sliding window approach was used to determine spatial biases of up/down regulation (see Section Material and methods). RNA-seq analysis of the exponential and stationary phase revealed the same spatio-temporal gene expression pattern seen in the reference study (see Figures 1A, 2A). Due to the high degree of similarity of the reference strain and the E. coli MG1655 parental strain, the reference strain is referred to as wild type (WT) for the rest of the manuscript.

For comparison of WT and INV transcription patterns, RNA-seq of the INV and WT strains was carried out in triplicates. In analogy to WT, we analyzed the expression profile of the INV strain in the exponential (EXP) and stationary (STAT) phase. For better comparison, the expression profile was mapped against the WT chromosome. Mapping against the INV chromosome would alter the coordinate system of the chromosome with respect to WT comprising altered replichores and inverted regions. Comparison and interpretation of the data would therefore be difficult. As a consequence of the sliding window approach, windows spanning the inversion break points are not present in both data sets and are therefore excluded from the analysis (see Section Material and methods and Figure 3A). As seen in Figure 3A, the spatial expression gradient of the INV strain along the oriC-ter axis is mostly attenuated compared to WT. Nevertheless, local characteristic peaks of the spatial pattern are still consistent with the wild type pattern suggesting a state of the promoter regulatory system similar to WT (see Supplementary Figure S2). This indicates that copy number effects may play a major role in the formation of the gradual expression pattern while promoter-specific regulation remains unchanged. More compellingly, for the left replichore in Supplementary Figure S3A, genes closer to the terminus show a higher expression in exponential phase compared to the stationary phase. Due to the inversion, these genes are situated close to oriC_INV in the INV strain. Although general copy number is reduced significantly due to the slower growth rate, the effect of inverted copy numbers in the inversion region is still reflected in the expression profile in the INV strain.

To verify the effects of copy number and study them isolated from other regulatory effects, exponential phase of WT and INV strains were compared in the same growth phase. While both strains differ in strength and orientation of the copy number, in the same growth phase, differences in promoter regulation are expected to be minimal. The comparison revealed a very strong expression bias gradient along the oriC-ter axis (see Figure 3B). The vast majority of oriC-proximal genes show a higher expression in the WT strain, while genes closer to the terminus are predominantly higher expressed in the INV strain. Interestingly, in the putative absence of difference in promoter regulation, the gradient is even more pronounced than between the exponential and stationary phase indicating that the copy number effect is the source of the gradual spatial expression pattern. Additionally, the characteristic local peaks seen in the analysis of WT EXP/STAT and INV EXP/STAT cannot be observed. This underpins the promoter regulatory origin of the local peaks between the exponential and stationary phase. Furthermore, the expression biases of both WT EXP/STAT (see Figure 2F) as well as the WT/INV EXP (see Figure 3B) follow their corresponding copy number differences (see Figures 3C, D). The expression bias even reflects the small steps in copy number differences at the inversion break points (see Figures 3B, D). Due to the inversion of the major part of the left replichore in the INV strain, copy number differences between WT and INV strain are not symmetrical on both replichores. The impact of the inverted copy number on the left replichore of the INV strain is rather small due to the relatively flat overall copy number gradient of the INV strain (see Figure 2C). However, the impact of the reverted copy number of loci on the left replichore is causing in a steeper copy number difference gradient of the left replichore compared to the right replichore in the comparison of WT/INV during exponential phase (see Figure 3D). The steeper gradient is also reflected in the gene expression bias gradient (see Figure 2F), indicating a close relation of gene copy numbers and the gene expression bias.

Even though the previous data suggest a major role of the copy number effect on the expression profile, the exact impact is not yet quantified. Whether other regulatory factors systematically contribute positively or negatively to the spatio-temporal expression pattern is still an open question. Multiple copies of a gene cause an increase in gene expression proportional to the number of copies. Consequently, the average expression fold change should match the corresponding copy number differences, if the copy number effect is the dominant factor. If other systemic regulatory factors influence the spatio-temporal pattern, the average expression fold change should deviate significantly from the copy number difference. As seen in Figure 3E, the average expression fold change of WT EXP/STAT corresponds well to the copy number difference (see Figure 3C). For the case of WT/INV EXP (see Figures 3D, F), where copy numbers were systematically reduced in the INV strain, the fold changes also matched the copy number differences supporting the role of copy numbers in forming spatial expression patterns. In this case, the characteristic local peaks (see Supplementary Figure S2) observed between exponential and stationary phase is flattened out, indicating a promoter regulatory source between these phases.

In certain cases, it is important to remove copy number effects from expression data. Such cases are mutant studies in which regulatory effects of the mutant are investigated. A growth defect often observed in regulator mutants would introduce a bias caused by copy number differences between wild type and mutant (see Supplementary Figure S5) (Beber et al., 2016). Consequently, gene expression data and deduced regulatory interactions might be biased. Therefore, we tested this scenario by first normalizing the WT exponential phase expression data by the copy number difference between WT and INV in exponential phase. If copy numbers are a major player, the removal of the copy number effect in WT with respect to INV should generate an expression pattern similar to INV in exponential phase. Consequently, EXP_WTcor/STAT_WT should resemble EXP_INV/STAT_INV. The comparison revealed a remarkable similarity (see Supplementary Figure S4). This shows that copy number data can be used to compensate for copy number differences between samples and underpins the impact of copy numbers on forming spatial expression patterns.

We have shown that the copy number effect plays a major role in forming a spatio-temporal gene expression pattern between exponential and stationary phase. This may also indicate a central role in the regulation of individual genes and pathways. However, a spatial bias induced by copy numbers does not necessarily imply a major role in single gene regulation. Promoter regulation may alter gene expression several hundred-fold (Lu, 1997). Regarding total fold change, the fold change of copy number can be neglected in such cases. To estimate the impact of copy number relative to other regulatory factors, we analyzed the single gene expression fold change data of WT EXP/STAT. The expression fold change of a gene is determined by its difference in copy number and promoter regulation. We have shown that on a large scale (300 genes average), the expression fold change follows the copy number between exponential and stationary phase. Therefore, we can remove the copy number effect (f_copy) of single genes from its expression fold change (f_total) to determine the influence of the remaining promoter regulation (f_reg).

We compared the copy number influence factor and the promoter influence factor of each individual gene and divided them into two categories: predominantly copy number regulated and predominantly promoter regulated, depending on the proportion of each factor on the total fold change of gene expression. About 40% of the genes are predominantly regulated by copy number when comparing exponential with stationary phase (see Figure 4A). For more than three quarter (78%) of genes, its copy number still covers more than one fourth of its total fold change. This underlines the broad relevance of the copy number effect in gene expression. However, there are also genes which have a significantly higher fraction of promoter regulation. For those genes, the copy number effect presumably plays a minor role in regulation.

When comparing the influence factors of WT/INV EXP, the majority of the genes (67%) are predominantly regulated by the copy number effect with a reduced influence of other regulators (see Figure 4B). In this case, about 92% of the genes have a copy number influence factor, which is greater or equal to one fourth of the influence factor of other regulators. This reflects the mild influence of promoter regulation and the copy number dominance in this experimental design. The remaining fraction of altered promoter regulation could originate from the altered expression of regulators including RpoS, Fis, Crp, or StpA located in the inversion region of the INV strain and its secondary effects.

The large set of genes dominantly regulated by copy numbers may indicate a concerted regulatory mechanism. In the case of WT/INV EXP, the experimental setup was rather artificial and did not follow a process in the life cycle of E. coli. Therefore, we focused on the WT EXP/STAT experiment to see if there is a link between regulation by copy number and specific cell functions. For genes predominantly regulated by the copy number, frequencies of functional categories were investigated (see Figure 4C). We found a significant over representation of genes in the “Coenzyme metabolism”, “Amino acid metabolism”, and “Trafficking/secretion” categories, while in “Energy production and conversion” genes regulated by copy number were under represented. For coenzyme metabolism, it was shown that in E. coli coenzyme synthesis is directly correlated to growth (Hartl et al., 2017). E. coli is capable to effectively adjusting de novo coenzyme synthesis to counteract varying dilution rates during growth, but the regulation is still unknown. Hence, a direct linking to replication rounds and therefore to copy number appears plausible. Also, amino acid metabolism is involved in biomass formation and is, therefore, a plausible candidate for coupling to copy numbers. Analysis of individual metabolic pathways further supports a coherent regulation by copy numbers. Here, specific pathways were strongly enriched in copy number dominated genes e.g., the aspartate pathway (see Supplementary Figure S3). Interestingly, in this pathway, the copy number regulation of intermediate pathway steps is complemented by promoter regulation at key steps at the entry and exit points. Lists with all metabolic pathways from E. coli and the corresponding copy number influence, as well as raw data for fold change analysis, can be found in the Supplementary material.

As coenzyme and amino acid metabolism are essential in bacteria, their genes are evolutionary conserved. With respect to their regulation by copy number, those genes may also exhibit conservation regarding the location on the chromosome. Genes, which are coupled to growth and copy number should be located close to the oriC or at least maintain their relative position to oriC. In order to investigate the evolutionary conservation of genes predominantly regulated by copy number, we first divided all genes into different sets depending on the extent of copy number regulation (see Figure 4D). Three opposing sets were generated with increased stringency for either copy number or promoter regulation dominance. We then estimated the variation of those genes in two species (Dickeya dadantii and Vibrio natriegens) with respect to their oriC distance in E. coli. D. dadantii is part of the Enterobacterales and is the causative agent of bacterial stem and root rot affecting potatoes and other crops, while V. natriegens is part of the Vibrionales and of increasing biotechnological relevance. The stronger the gene regulation is dominated by the copy number, the less variation in oriC distance is observed in these two species compared to E. coli. In contrast, the stronger their regulation is dominated by promoter regulation, the more variation is detected. This indicates an oriC distance conservation of genes regulated by copy number and high spatial flexibility of genes governed by promoter regulation. This result refines the general observation of conserved positioning of genes with respect to the oriC distance in bacteria (Sobetzko et al., 2012).

Selective pressure for copy number regulation could depend on the extent of the available copy number effect. Consequently, for faster growing species, a higher portion of genes active during exponential growth could exploit copy number effects for regulatory purposes and migrate toward oriC. The two selected species flank E. coli with respect to doubling time during exponential growth. D. dadantii (~100 min) has a longer doubling time than E.coli (approx. 20 min) whereas V. natriegens (approx. 10 min) exhibits a far shorter doubling time. As DNA polymerase speed is a limiting factor for fast growing bacteria, a reduced doubling time is reflected in intensified overlapping replication increasing copy numbers (see Figure 4E). Consequently, in these species, three different levels of copy numbers are realized allowing us to investigate the impact of copy numbers on chromosome evolution with respect to gene location. For comparison of the three species, genes predominantly expressed during exponential phase in E. coli (p-value <0.05) were selected. Due to the difference in expression between exponential and stationary phase, copy number can potentially positively regulate these genes also in other organisms. The distribution of orthologs in the three species reveals a link between growth rate and stringency of gene positioning (see Figure 4F). Orthologs in the slow growing bacterium D. dadantii were less focused on the oriC proximal region than orthologs in E. coli. Orthologs in V. natriegens, the fastest growing bacterium, were even more focused in the oriC proximal region than in E. coli. We further investigated this observation using a larger set of species. For most species, the oriC position is not determined (Gao and Zhang, 2007) or hidden in countless publications. Therefore, we devised a method that is based on the oriC-ter symmetry found across the bacterial kingdom (Sobetzko et al., 2012). Hence, the chromosomes of two species match best with respect to oriC distance if oriCs of both species are superimposed (see Supplementary Figure 8 and Section Material and methods “constellation analysis”). Using this approach, oriC distance correlation was determined for species pairs of various phylogenetic classes including gram-positive and gram-negative bacteria. The species within these classes were selected to be different in their family membership to ensure a standardized evolutionary distance. To approximate growth rates, we used the correlation of growth rate and the number of ribosomal operons of a species. This correlation was first verified using growth rates of Couturier and Rocha (2006) and NCBI 16S rRNA annotations (see Figure 5A). Species in each class were split into slow growing (16S rRNA ≤ 3) and fast growing (16S rRNA ≥ 6). For all species pairs in these sets, oriC distance correlation was determined. Consistent with the initial analysis in Figure 4F, fast growing species showed a stronger correlation of oriC distance between its orthologs than slow growing species in all investigated classes (see Figure 5B). This indicates that copy number regulation is also involved in the evolutionary shaping of bacterial chromosomes proportional to its regulatory potential.

The INV and WT strains were cultivated in LB medium (10g/L tryptone, 5g/L yeast extract, 10g/L NaCl) at 37°C under aerobic conditions in flasks with shaking at 200 rpm. For the analyzes of the exponential phase, both strains were harvested at an OD₆₀₀ of 0.3 pelleted and immediately suspended in RNA later (Thermo Fisher) according to the manufacturer's instructions. The stationary phases of the WT and INV strains were harvested when no significant changes in OD₆₀₀ could be observed for about 20 minutes. Subsequently, cells were also pelleted and suspended in RNAlater. Samples were then split for DNA- and RNA-sequencing. Isolation of bacterial genomic DNA was performed according to Bruhn et al. (2018). For RNA-sequencing, lysis of cells and subsequent isolation of total RNA were carried out using the lysing matrix B/FastPrep sample preparation system (MP Biomedicals) and the miRNeasy Mini Kit (Qiagen), respectively. Ribosomal RNA depletion (RNA) and library preparation (RNA/DNA) was conducted by Eurofins Genomics using the Illumina Technology (strand-specific; paired-end; 2x150bp read length). All samples (INV/WT exponential and stationary phase) were carried out in biological triplicates. For the comparative genomics analysis of three species (Dickeya dadantii 3937; Escherichia coli MG1655; Vibrio natriegens ATCC 14048), all species were grown in rich medium. For D. dadantii and E. coli LB medium was used. For the halophilic V. natriegens LBV2 (LB + 204mM NaCl, 4.2mM KCl, 20.14mM MgCl₂) was used. All three species were grown under aerobic conditions in baffled flasks with orbital shaking at 200 rpm. For optimal growth, E. coli and V. natriegens were grown at 37°C and D. dadantii at 30vC. Cells were harvested at OD₆₀₀=0.3 in mid-exponential phase. DNA was also extracted according to Bruhn et al. (2018) and Illumina-sequenced by Eurofins Genomics yielding 5M 150bp paired-end reads.

Data concerning sigma factor and transcription factor regulation was obtained from regulonDB (v10.9). For chromosomal sigma factor distribution, genes solely regulated by the respective sigma factor were selected. Data concerning functional identity were derived from the NCBI COG database (https://www.ncbi.nlm.nih.gov/research/cog-project/).

For the construction of the inversion strain (INV) to dissect the impact of promoter regulation and copy number effect as well as for reference strain (WT) we used CRISPR SWAPnDROP to make all relevant chromosomal changes in E. coli MG1655 (Teufel et al., 2022). For each deletion and insertion, a different (pSwap) plasmid was constructed harboring homology regions, sgRNAs and inserts. All primers used for the amplification of the homology regions, for the sgRNA construction as well as for each insert are available in the Supplementary material. An overview of each chromosomal edit done in each of the strains is given in Table 1.

DNA read mapping was done with the R QuasR package. Marker frequency analysis (MFA) was performed to measure copy number (Maduike et al., 2014; Ivanova et al., 2015; Kemter et al., 2018). Genome coverage of exponential phase samples was first averaged over 5kb sliding windows relative to the corresponding stationary phase to get robust estimates of local copy numbers (see data points in MFA plots). A log₂ linear regression of local copy numbers was performed for each replichore separately. The intersection ordinate of the two replichore regression curves was used as oriC and terminus (ter) copy number estimates. The data were normalized to a terminus copy number of 1 to simplify illustrations. For copy number difference, the fold change between the regression curves of the investigated growth phases and strains was calculated at the corresponding locus. For copy numbers of individual genes, the ordinate of the respective replichore regression curve at the gene locus was used.

RNA-sequencing reads of each gene were first normalized for gene length and the total number of reads in each sample. All samples of this study were quantile-normalized in one batch to harmonize differences in overall gene expression distributions caused by technical variation. For the expression bias analyzes, the differences in the number of up- and down-regulated genes between the growth phases were determined for sliding windows of 300 genes. The chromosomal location of each window was set to the average location of all genes in the window. For the fold change analyzes, the average gene expression fold changes between the growth phases were determined for sliding windows of 300 genes. As the spatial expression pattern represents a systematic bias in expression data, the average fold changes were further corrected for relative frequency biases. This systemic bias originates from the imbalance of relative frequency when one component is enriched leading to a depletion of all other components. In the concrete case, copy number causes an increase of oriC-proximal gene expression levels, which in turn reduces expression levels of the terminus-proximal genes. This results in negative average fold change values in the terminus region. However, we assume the average fold change to be 1 at the chromosomal location where no copy number difference is present between samples. Therefore, we corrected the fold changes accordingly. The location of no copy number difference was first extracted from the copy number difference curve. The average fold change bias at that location was determined by taking the ordinate value of a regression curve of the fold change data on both replichores. All fold changes were corrected for that ordinate value. In the spatial analyzes of WT and INV, reference chromosome coordinates were set to wild type. The inversion in the INV strain causes new neighborhoods of genes at the break points of the inversion. Therefore, windows spanning these breakpoints contain gene sets that are not present in WT (e.g., WT oriC-proximal genes paired with terminus-proximal genes). These windows were omitted in the analysis as no counterpart was present in WT.

Orthologs of genes in the three species were determined using proteinortho v6 (Lechner et al., 2011). Only orthologs with a single copy (no paralogs) in all three species were considered. Gene positions were transformed to relative oriC proximities and normalized by half of the chromosome size. Consequently, values range between 1 (oriC) and 0 (ter).

For the analysis of variability of gene position, promoter and copy number regulated genes were split into 3 subsets with increasing stringency of the regulatory type (promoter or copy number). For each gene in a set the difference of oriC proximity was determined for E. coli vs D. dadantii and E. coli vs V. natriegens. The two resulting differences were averaged to reduce species-specific biases. Then, the average and standard errors of these averaged differences were determined for the different sets.

For the comparison of ortholog positions in the three species, genes significantly up-regulated (p-value <0.05) in exponential phase relative to stationary phase of E. coli WT expression data were chosen.

The significance of functional groups was determined by generating 1,000 random sets of genes of the same size as the set of predominantly copy number regulated genes. For these sets mean frequencies m and the corresponding standard deviations s for functional groups were determined to compute a z-score z

for each functional group, where x is the number of genes in the respective function group of the experimental data.

The full NCBI set of completely assembled genomes was first screened for NCBI taxonomy information to cluster species according to phylogenetic categories. The remaining set was split into distinct phylogenetic classes that were analyzed separately. To cover the diversity of a class and avoid a representation bias, one species was selected out of every family of the class. Each species was a randomly selected representative species of its family, according to the NCBI database. The family category was chosen to select species with a defined range of evolutionary distance for later comparison. Furthermore, categories below a family with closer evolutionary distance yielded little chromosomal diversity for analysis. Species within a class were split into the set of fast or slow growing species. This was done by the number of 16S rRNAs that correlate with growth speed. The number of 16S rRNA was extracted from the species annotation files (GFF3) provided by NCBI. Data about the doubling times were taken from Couturier and Rocha (2006). For all species pairs within a set, orthologs were determined using proteinortho v6 (Lechner et al., 2011). Only orthologs with a single copy (no paralogs) in the two species of a pair were considered to determine reliable chromosome positions of orthologs. The oriC-ter axis was determined by finding the putative oriC position in both species that give rise to the best Pearson correlation coefficient of distances to oriC (see Section Constellation analysis). These correlation coefficients were used as an indicator of the strength of oriC-ter axis symmetry.

For most species, oriC position is not determined. However, the oriC position causes a chromosomal symmetry due to positional conservation of genes relative to oriC (Sobetzko et al., 2012). This can be used to determine the oriC-ter axis. To determine the axis, the best overlay of chromosomes of two related species is determined. First, an arbitrary oriC position is assigned individually to both chromosomes (e.g., position 100,000 for species 1 and position 500,000 for species 2). The relative distances to all genes on both chromosomes to the assigned oriCs are determined (see Figure 5A top row). For all orthologous gene pairs of the two species, the distance to the respective oriC is compared yielding a correlation coefficient (Pearson correlation coefficient). This approach is repeated with other arbitrary oriC positions until all combinations were tested (see Figure 5B). The positions yielding the highest correlation are taken to be oriC. In (see Figure 5B) several optima are present due to intrinsic symmetries. For example, taking the correct oriC position of E. coli and D. datantii yields the same maximal correlation as taking the correct terminus position of E. coli and the terminus position of D. datantii as this directly infers that oriC is also superimposed. Which of the two optima are chosen has no impact on further analysis in this study as the level of the correlation is central not the related positions that yielded the maximal correlation. To distinguish which optimum is associated with oriC and ter, additional information such as dnaA or ribosomal operon location can be used. However, this is not relevant for the study. The strength of the correlation is an indicator of the conservation of gene position relative to oriC. For the analysis, it is important to consider the total evolutionary distance, as closely related species, in general, show a higher degree of conservation. This can be accomplished by choosing species within the same phylogenetic ranges (e.g., class, family, phylum etc.). Furthermore, there is a lower limit for this method, when little chromosomal rearrangements occurred between species. In this case, constellation analysis indicates this by a diagonal line of similarly high correlations. This is caused by the fact that correlation remains stable when rotating putative oriC positions on both chromosomes simultaneously (e.g., moving in 5000 bp steps in a clockwise direction). If gene positions are similar on both chromosomes, also distances to the oriCs change accordingly during synchronous rotation and retain the same correlation coefficient. The more genes have moved to other locations e.g., the opposite replichore, the more correlation differs during synchronous rotation. In Figure 5B a faint diagonal line can still be seen, indicating a small percentage of orthologous genes with a similar chromosome position in both species. As long as the two optima can still be clearly determined, the method is applicable.