Strains

B. subtilis strains were generated from 168 (bearing the trpC2 mutation), and E. coli genomic DNA was extracted from MG1655. Linearized plasmids and genomic DNA were transformed using either standard protocols relying on natural competence60 or supercompetent B. subtilis strains bJD086 (xylose inducible, wild-type) and bJD087 (xylose inducible, Δrho::Kan) based on ref. 61, or bCE017 (anhydrotetracycline (aTc) inducible, wildtype) and bCE019 (aTc inducible, Δrho::Kan) based on ref. 62. bJD086 and bJD087 contain the xylose-inducible comK induction system from SCK6 (ref. 61) inserted at amyE. bCE017 and bCE019 contain the aTc-inducible comK system from (ref. 62) inserted at amyE.

Sequences of plasmids and integrations were confirmed by Sanger or nanopore sequencing. All strains are listed in Supplementary Table 2. All plasmids were generated in E. coli DH5-α cells using standard protocols and are listed with additional details in Supplementary Table 2.

Genomes

The RefSeq63 entries for B. subtilis (NC_000964.3) and E. coli (NC_000913.3) were used to map the sequencing reads and generate genome fragments in silico, and the corresponding annotations were used to classify fragments overlapping with sequences sense or antisense to genes.

List of oligonucleotides

Supplementary Table 2 contains the list of oligonucleotides used for strain construction and sequencing library generation.

Supercompetent B. subtilis transformation

To clone individual B. subtilis strains and generate the B. subtilis reporter library, linearized plasmid was introduced into B. subtilis using an inducible competence system mediated by the xylose-inducible comK allele from ref. 61 or the aTc-inducible comK allele from ref. 62 integrated at the amyE locus. Briefly, a colony of bJD086 or bCE017 (wild-type background), or bJD087 or bCE019 (Δrho background) was picked into LB and grown to OD600 0.8–1.1 at 37 °C with vigorous shaking. At this point, 1% xylose (w/v) (for bJD086 and bJD087) or 10 ng ml−1 aTc (for bCE017 and bCE019) was added to induce expression of ComK. Induced cultures were shaken at 37 °C for an additional 1.25–2 h, at which point 200 ng linearized plasmid was added per 100 μl induced cells. The cells were cultured with DNA for an additional 1–1.5 h before plating on selective plates for overnight growth at 37 °C.

Reporter plasmid library construction

The plasmid library was constructed through isothermal assembly of tagmented genomic DNA fragments into the reporter plasmid backbone, pJD16. To generate the inserts, a 1:1 mixture of B. subtilis and E. coli genomic DNA was fragmented by Nextera XT tagmentation (Illumina). Then, 6 ng of the tagmented genomic DNA was PCR amplified with oJD173/oJD174 for 10 cycles using Q5 DNA polymerase (New England Biolabs). PCR products 200–400 bp in length were size selected using an 8% TBE polyacrylamide gel (Thermo Fisher). This size-selected DNA was then cleaned up with a DNA Clean & Concentrator-5 column (Zymo Research). To linearize the pJD16 backbone, the pJD16 plasmid was PCR amplified with oJD091/oJD092 using Q5 DNA polymerase. The PCR product was extracted from an agarose gel using a Zymoclean Gel DNA Recovery Kit (Zymo Research) and DpnI (New England Biolabs) digested for 60 min at 37 °C. The DpnI-treated DNA was cleaned up with a DNA Clean & Concentrator-5 column.

To construct the plasmid pool, 200 ng of the pJD16 backbone PCR was mixed with insert at a 1:4 molar ratio and assembled at 50 °C for 15 min in a Gibson Assembly reaction (New England Biolabs). This reaction was then cleaned up with a DNA Clean & Concentrator-5 column and transformed into electrocompetent E. coli (New England Biolabs 10-β cells). Transformants were plated on LB with 100 μg ml−1 carbenicillin in 245-mm square bioassay dishes (Corning) and grown at 37 °C overnight. Cells were then scraped off the plates for ZymoPure II Maxiprep plasmid extraction (Zymo Research).

Transformation and integration of reporter library into the B. subtilis genome

The supercompetent transformation protocol described above was used to integrate the reporter library into B. subtilis strain bJD086, replacing the comK induction cassette at the amyE locus. Briefly, a colony of bJD086 was picked into 22 ml LB and induced with 1% xylose (w/v) once the culture reached OD600 1.02. After 1.5 h of growth in xylose, 15 ml of induced culture was combined with 1.5 ml (30 μg) of plasmid pool linearized with ScaI-HF (New England Biolabs). This digest was prepared according to the manufacturer’s protocol with 1 μg plasmid per 50 μl reaction. After 1.5 h of incubation with DNA, cells were pelleted, resuspended in residual media and plated on 245-mm bioassay dishes containing 100 μg ml−1 spectinomycin. After growth overnight, cells were scraped into LB with 100 μg ml−1 spectinomycin and approximately 500 million cells were back-diluted into LB with 100 μg ml−1 spectinomycin. This culture was grown for 5 h, and 1-ml aliquots of culture were mixed 1:1 with 40% glycerol and frozen at −80 °C.

Library growth and collection

To perform chloramphenicol selection on the B. subtilis library to enrich for genomic fragments with transcription termination activity, the library was grown in LB for five generations, split into +BCM-Bz and DMSO control (−BCM-Bz) conditions for pre-treatment and then back-diluted into selective media containing chloramphenicol and BCM-Bz or chloramphenicol and DMSO. Additional details on the timing and OD600 measurements for this process are described in the Supplementary Methods. In total, cell pellets for four conditions were collected: +BCM-Bz culture before and after chloramphenicol selection, and −BCM-Bz culture before and after chloramphenicol selection.

Genomic DNA sequencing

Sequencing libraries were prepared to quantify the frequency of each genomic fragment variant in the four experimental conditions. Genomic DNA was extracted from the cell pellets using the Wizard Genomic DNA Purification kit (Promega) following the manufacturer’s instructions. To reduce RNA contamination, the genomic DNA was then cleaned up with the Genomic DNA Clean & Concentrator-10 (Zymo Research). Libraries were generated using a two-step Q5 PCR protocol. To attach unique molecular identifiers, a two-cycle PCR was performed using 2 μg of genomic DNA per 150 μl PCR (scaled as needed for library complexity) with primers oJD180/oJD182. This PCR was cleaned up with two rounds of magnetic bead clean-up at a 1:1 ratio of beads to sample (PCRClean DX bead, Aline Biosciences). Half (preselection libraries) or one-quarter (post-selection libraries) of this cleaned-up reaction was then amplified in a second PCR with primers oJD181/oJD183. The second PCR was gel purified using an 8% TBE gel. These libraries were sequenced with 50-bp paired-end reads on a G4 sequencer (Singular Genomics).

Fragment quantification

Quantification of genomic fragments in each sample was handled using custom Python scripts. First, the paired-end sequencing reads were aligned to the B. subtilis (NC_000964.3) and E. coli (NC_000913.3) genome using the bowtie package (version 1.2.3)64. Then, 40 nt of each read was used for alignment, and alignments with more than 1 mismatch or an insert size of more than 5,000 bp were discarded (bowtie arguments: -trim3 10 -v 1 -X 5000). Any two reads with the same aligned start and end positions and unique molecular identifier were collapsed to a single read. The number of unique reads mapping to each genomic fragment was calculated and normalized to the total number of reads in the sample. The enrichment of each fragment in each condition (+BCM-Bz or −BCM-Bz) is the ratio of its frequency in the post-selection sample to the corresponding preselection sample. Reads mapping to the E. coli genome were combined and used for normalization to total reads in the enrichment calculations, but were not included in further analysis. Raw read counts for all fragments and samples are available in Supplementary Table 4.

Thresholding and pseudocounting enrichments

To reduce noise in the enrichment values, a threshold of ≥50 reads per fragment was applied to the preselection samples and a threshold of ≥25 reads per fragment was applied to the post-selection samples. Fragments that fell below the post-selection read threshold owing to depletion during selection were assigned a pseudo-enrichment value of 0.015 (labelled as below detection in Figs. 1d and 2a, and Extended Data Figs. 1 and 3d) provided that their raw enrichment was below 0.5 (−BCM-Bz) or 1 (+BCM-Bz). A total of 161,256 fragments passed these thresholds or were assigned a pseudo-enrichment in both conditions (+BCM-Bz and −BCM-Bz), and their raw read counts, normalized read counts and reported enrichments are available in Supplementary Table 4. In Figs. 1d and 2a, fragments with a pseudo-counted enrichment were assigned a random enrichment value to jitter points for visualization.

Fragment classification

Two criteria were used to identify Rho-terminated fragments from their enrichment scores. First, fragments terminated by Rho should be enriched in the absence of BCM-Bz, when Rho is active. Rho-terminated fragments had to therefore exceed an enrichment of 1 in the −BCM-Bz condition, which roughly corresponds to ≥65% termination efficiency (Extended Data Fig. 1). Second, the enrichment observed in the absence of BCM-Bz should be lost when Rho is inactivated by BCM-Bz treatment. Thus, the enrichment of Rho-terminated fragments in the +BCM-Bz condition had to be more than twofold lower (0.45×) than the enrichment in −BCM-Bz. This cut-off was chosen to minimize misclassification of fragments with known intrinsic terminators as Rho-dependent termination sites (less than 1.6% misclassified). A total of 9,995 fragments were classified as driving Rho-dependent termination activity, the positions of which are available in Supplementary Table 4.

Fragments referred to as ‘non-terminated’ throughout the main text and figures were defined as those with an enrichment less than 1 in the −BCM condition. The positions of these fragments are identified in Supplementary Table 4. (Note that these fragments may include some cases of weak termination and are further separated in the classifications described below.)

Criteria for the additional classifications shown in Supplementary Fig. 1 and a discussion of potential limitations in fragment classification are provided in the Supplementary Methods.

Fragment overlap with sense and antisense regions and intrinsic terminators

Fragments overlapping with sense and antisense regions of the B. subtilis genome were identified using the annotated RefSeq coding sequence (CDS) features for NC_000964.3. Sense and antisense fragments were required to fully overlap with the same or opposite strand of an annotated CDS, and any fragments that overlapped with multiple coding regions (on either strand) were discarded. In total, 47,550 sense and 54,692 antisense fragments of B. subtilis are shown in Fig. 2a. Fragments were annotated as containing an intrinsic terminator if any of the wild-type positions of intrinsic termination identified in ref. 34 fell within the fragment (n = 3,037 fragments). The termination efficiencies used in Extended Data Fig. 1 were determined from the readthrough fraction in Rend-seq of ΔpnpA B. subtilis34.

Excess C and maximum %T metrics, linear model and Rho target score

To identify sequences with a skewed C content relative to G, we wanted to capture information about the magnitude of the difference between C and G counts in a defined region (representing a single Rho termination site) and also capture the cumulative effect of multiple regions where C content exceeds G content, which we hypothesized would increase the probability of Rho termination in the fragment. Both features are captured by the excess C score, which is described in detail in the Supplementary Methods.

Maximum %T was determined by counting the number of threonines (Ts) in every 150-nt window tiling a sequence fragment and then taking the maximum T count as a fraction of the window length. We found that the performance of metrics related to T content in distinguishing the Rho-terminated fragments was not improved by incorporating information about relative adenine content. Fragments shorter than 150 bp (2.1% of the input library) were excluded from sequence feature analysis.

To evaluate the explanatory power of the excess C and maximum %T metrics and predict enrichments in silico, a linear model was trained to predict fragment enrichments from these two scores. The ‘Rho target score’ describes the enrichment predicted by this model. The training process and coefficients for this model and corresponding boundary line are described in the Supplementary Methods.

LacZ reporter strains

DNA fragments with potential Rho-dependent transcription termination sites were cloned upstream of lacZ, under the control of the IPTG-inducible pSpankHy promoter in the vector pJD19. Detailed methods describing the process to generate the LacZ strains are available in the Supplementary Methods.

For qualitative assays of β-galactosidase activity, blue–white screening was performed on X-gal (5-bromo-4-chloro-3-indolyl β-D-galactopyranoside) (GoldBio) plates. For each strain, a colony was picked into 5 ml LB and grown at 37 °C with vigorous shaking to OD 0.50–3.5. These cultures were then diluted to OD600 0.005 in LB, and 3 μl of the dilutions were spotted onto LB-agar plates with 1 mM IPTG and 200 μg ml−1 X-gal. Plates were incubated at 37 °C overnight and imaged 16–17 h after plating.

The process for generating recoded brnQ and HGH sequences, and additional details on the brnQ homologue sequences, are provided in the Supplementary Methods.

Rho homologue classification

For Lentilactobacillus buchneri and the genomes shown in Fig. 6b, classifications of genomes with and without a rho homologue were taken from ref. 34. For all other Bacilli species (Fig. 3c–e and Extended Data Fig. 7), the Rho homologue classifications in ref. 42 were used.

In silico generation of B. subtilis and E. coli genome fragments

To analyse the nucleotide composition of sense and antisense sequences in the B. subtilis and E. coli genomes, 300-nt tiling fragments (step size: 50 nt) were generated from the sequence of all CDS features in the RefSeq annotation files. Features shorter than 300 nt were excluded. The antisense fragments were generated by taking the reverse complement of the sense sequence fragments. All fragments (n = 51,531 for B. subtilis and 57,477 for E. coli) were used to generate the Gaussian kernel density estimates shown in Fig. 3a, and a subset of 1,600 fragments, split evenly into sense and antisense strands, are shown in the scatter plot. For the unbiased leading and lagging strand fragments shown in Fig. 4, randomly selected start positions were used to generate 300-nt windows of genomic sequence.

For Fig. 4, the B. subtilis sense sequence fragments and randomly generated fragments were assigned to the leading or lagging strand based on their position relative to the origin (4,215,389 (ref. 65)), assuming replication arms of equal length.

Bacilli phylogenetic tree

The phylogenetic tree of Bacilli (Fig. 3b) was generated from the Genome Taxonomy Database (GTDB) tree66 (version 226, downloaded June 2025). A total of 155 Bacilli species analysed in ref. 42 and L. buchneri were included in the tree, with Rho homologue classification based on species name using refs. 34,42 as described above. Only a subset of the 309 Bacilli species shown in Fig. 3c was included in this tree for visualization purposes. See the Fig. 3 source data for the full list of species included in the tree.

Sense strand purine content in Bacilli and across Bacteria

To compare the purine content in sense sequences across Bacteria (Fig. 6b), we analysed a set of 1,551 RefSeq genomes with coupling predictions10 and an annotated origin of replication67. To correct for length variation across genes, the purine content was determined for a random 250-nt fragment of each annotated gene (genes shorter than 250 nt were discarded). Each gene was then assigned to the leading or lagging strand using the position of the origin of replication from ref. 67 and assuming replication arms of equal size. The average purine content was calculated separately for genes on the leading and lagging strands, and the leading strand averages for the 1,002 genomes that encode a rho homologue and were assigned to a named phylum in ref. 10 are shown in Fig. 6b. The designation of runaway versus coupled transcription was performed at the phylum level, based on ref. 10. Species in Bacillota, Fusobacteria, Thermotogota, Deferribacterota, Aquificota and Campylobacterota were classified as exhibiting runaway transcription. These data are available in the Fig. 6 source data.

The same process for analysis of sense sequence purine content was repeated for the 308 Bacilli genomes (GTDB classification of c__Bacilli) from ref. 42 and L. buchneri to generate the data for Fig. 3c.

Codon usage

For each of the Bacilli genomes (GTDB classification) with a rho homologue assignment from ref. 42, the nucleotide sequences of the annotated gene sequences (CDS features) were downloaded from the National Center for Biotechnology Information (NCBI) RefSeq database63. The frequency of each codon across the entire coding genome was then normalized to the amino acid frequency to determine the codon usage by amino acid. To determine the frequency of codons with a purine in the third position, the frequency of codons ending with A and G were summed for each amino acid. For each amino acid considered, the purine frequency in the third codon position for each Bacilli species is available in the Fig. 5 source data.

Extended data and supplementary figures

Methods related to all extended data and supplementary figures are provided in the Supplementary Methods.

Statistics

To generate the P values reported in Fig. 5 and Extended Data Fig. 7, a one-sided Mann–Whitney U test was performed. Assuming underlying distributions of the same shape, the alternative hypothesis for the test was that the median relative frequency (of codon or amino acid usage, respectively) was higher for Bacilli with rho than for species without rho. The values for the U statistic are as follows: Fig. 5: T (19,385), V (16,585), A (18,872), G (16,243) and P (13,820), and Extended Data Fig. 7: E to D (8,853), E to Q (8,507) and Y to F (6,056). Cohen’s d was used to report the effect size of the difference in means (absolute difference between the means divided by pooled standard deviation) and is reported in the corresponding figure. The number of each species is available in the figure legend, and source data are provided.

For Extended Data Fig. 3a,b, the Area Under the Curve statistic was computed for the distribution of the maximum CG skew or the excess C score, respectively, in non-terminated fragments (n = 136,728 fragments) versus Rho-terminated fragments (n = 9,983 fragments). Fragments shorter than 150 nt were excluded from the analysis. Source data are provided.

A description of the calculation of the R2 statistic reported in the text and Extended Data Fig. 4 is available in the Supplementary Methods in ‘Linear model and Rho target score’.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.