RiboRid: A low cost, advanced, and ultra-efficient method to remove ribosomal RNA for bacterial transcriptomics

RNA sequencing techniques have enabled the systematic elucidation of gene expression (RNA-Seq), transcription start sites (differential RNA-Seq), transcript 3′ ends (Term-Seq), and post-transcriptional processes (ribosome profiling). The main challenge of transcriptomic studies is to remove ribosomal RNAs (rRNAs), which comprise more than 90% of the total RNA in a cell. Here, we report a low-cost and robust bacterial rRNA depletion method, RiboRid, based on the enzymatic degradation of rRNA by thermostable RNase H. This method implemented experimental considerations to minimize nonspecific degradation of mRNA and is capable of depleting pre-rRNAs that often comprise a large portion of RNA, even after rRNA depletion. We demonstrated the highly efficient removal of rRNA up to a removal efficiency of 99.99% for various transcriptome studies, including RNA-Seq, Term-Seq, and ribosome profiling, with a cost of approximately $10 per sample. This method is expected to be a robust method for large-scale high-throughput bacterial transcriptomic studies.

Introduction Genetic information encoded in the genome is transferred to proteins via messenger RNAs (mRNAs). Thus, investigating mRNAs is a central approach to elucidating the fundamentals of cellular functions. Multiple techniques such as quantitative polymerase chain reaction (qPCR), microarray, and RNA sequencing (RNA-Seq) have been developed to quantitatively measure mRNAs inside a cell or their changes in response to a variety of environmental and genetic perturbations [1,2]. Since ribosomal RNAs (rRNAs) comprise more than 90% of the total RNA in a cell, their efficient removal is one of the most important tasks for reliable genome-wide transcriptomics studies [3]. Unlike eukaryotic mRNAs that can be selectively enriched by virtue of their poly A tail structure [4], bacterial and archaeal mRNAs do not, or rarely, possess such features; thus, removal of the predominant rRNA species from total RNA is critical for downstream applications.
Several methods for removing bacterial rRNAs have been developed and commercialized. Subtractive hybridization with nucleic acid probes that are reverse complementary to rRNAs is the most popular method that has been commercialized in multiple systems, such as Ribo-Zero, MICROBExpress, RiboErase, and RiboMinus [5][6][7]. However, as the listed methods depend on short conserved regions on rRNAs, they tend to be inconsistent in partially degraded RNA samples and species that have divergent rRNA sequences compared to the probe sequence. Exonucleolytic digestion of processed RNA with monophosphate at the 5 0 end has also been devised for transcriptomics [8]. However, this method has relatively low efficiency and is limited by the fact that primary transcripts protected by triphosphate from 5 0phosphate-dependent terminator exonuclease (TEX) may not be a precise representation of mRNA levels in a cell, because a considerable amount of mRNA exists as a processed form. Several methods have been proposed for rRNA removal based on duplex-specific nucleasebased digestion, electrophoretic size selection, and sequence-specific blockage of reverse transcription, although they are not as efficient as commercial systems [9][10][11][12].
Ribonuclease H (RNase H) is an endoribonuclease that specifically digests RNA strands of RNA:DNA hybrids. The RNase H-based selective digestion of rRNAs has been proposed as a cost-effective alternative method for depleting cellular rRNAs because of its high removal efficiency [3,13]. An RNA-Seq study of formalin-fixed paraffin-embedded (FFPE) archival tissue using RNase H was reported previously [3]. However, this method was not revisited in bacterial samples until the recent temporal discontinuation of the Ribo-Zero. Although an rRNA depletion method was assessed for bacterial RNA-Seq as a low-cost alternative [13], additional optimization is required for efficient rRNA depletion. These optimization methods include the following: (i) optimization of reaction conditions to prevent nonspecific hybridization; (ii) design of the oligonucleotide probes complementary to the rRNAs; (iii) addition of probe specific for premature rRNAs (pre-rRNAs) that contain extended sequences from the operonic structure of rRNA operons; and (iv) application of the method to other transcriptomic studies such as sequencing the 3 0 ends of bacterial transcripts (Term-Seq) and ribosome profiling. To this end, we revisited the rRNA depletion method to improve and streamline the experimental procedure by adjusting the reaction temperature, ionic strength, removal of pre-rRNAs, and redesigning probe sequences. The advanced method presented here, called RiboRid, could remove bacterial rRNAs up to 99.99% and prevent nonspecific binding of oligonucleotide probes to other cellular RNA species. Furthermore, RiboRid has rRNA removal efficiency comparable to the most efficient commercial method (Ribo-Zero) when performing technically challenging transcriptomic studies, such as Term-Seq and ribosome profiling.

Sequence specific digestion of ribosomal RNAs
Selective digestion of rRNAs using RNase H depends on hybridization of the short deoxynucleoside probes specific to rRNAs and digestion of the RNA:DNA hetero-duplex by RNase H. Thus, minimization of the nonspecific hybridization of the probes to other RNA species is critical for efficient and nonbiased transcriptome studies. Thus, the hybridization reaction of the Ribo-Zero method is carried out at a high temperature (68˚C) to reduce nonspecific hybridization. To this end, we aimed to optimize the experimental conditions in the RNase H-based rRNA depletion method to prevent unwanted RNA degradation by nonspecific binding (Fig  1). First, the RNA sample was denatured at 95˚C for only 1 s. Prolonged incubation at a denaturing temperature that may induce spontaneous hydrolysis of RNA was not necessary for all tested samples in this study. Next, we performed RNase H reaction at an elevated temperature (65˚C) to prevent unwanted hybridization of anti-rRNA oligonucleotide probes (ArOPs) to rRNA mRNA mRNA. To support this reaction, Hybridase thermostable RNase H, that has an optimum reaction temperature at 65˚C or higher, was used. Synthetic single-stranded deoxyoligonucleotides with melting temperatures higher than 68˚C were used as ArOPs to support efficient hybridization at a high reaction temperature (65˚C).
We used 32 nt-long single-stranded deoxyoligonucleotides as ArOPs that can bind to rRNA. When the melting temperature of the probe is lower than 68˚C, the position of the probe is moved to the rRNA until the melting temperature reaches 68˚C. However, when an appropriate probe could not be found or the adjusted probe overlapped with the adjacent probe, the length of the probe was shortened to 28 nt or extended up to 43 nt (see Methods). The probes were designed to bind approximately every 40 bp of rRNA and hybridized to rRNA in an RNA sample that had been denatured. Then, the RNA strand of RNA:DNA hetero-duplexes were subjected to RNase H digestion. Compared to the previous approaches, that use mild temperatures (22-45˚C) for the digestion reaction [3,13], RNA is relatively prone to spontaneous hydrolysis at the elevated reaction temperature in this study. The nonenzymatic RNA hydrolysis occurs at elevated temperatures due to the presence of magnesium cations (Mg 2+ ), which is required for RNase H activity [14]. Thus, the reaction buffer was modified to contain a lower Mg 2+ concentration, which is 6.4-fold lower than that in previously reported methods [3,13]. With the optimized protocol, we could detect removal of rRNA by electrophoretic analysis (S1 Fig), which can be a useful checkpoint before proceeding to downstream study. Subsequently, ArOPs were removed by DNase I digestion followed by column-based purification, to collect remaining mRNAs. This additional enzymatic digestion can be avoided when the 3 0 end of ArOP is blocked by C3-spacer or dideoxynucleotide modification, because they are not used as a substrate by ligases and polymerases; thus, they do not interfere with downstream library preparation.

Streamlining the rRNA removal protocol for RNA-Seq in Escherichia coli and closely related species
We first examined the capability of the RiboRid method for removing rRNAs from 500 ng total RNA samples from Escherichia coli and compared it with three commercial methods, MICROBExpress, Ribo-Zero (a legacy system that uses subtractive hybridization), and RiboErase. Without rRNA depletion, rRNA comprised over 97% of the total sequenced reads ( Table 1; untreated). Commercial kits showed successful removal of rRNA, except for the MICROBExpress method that showed low removal efficiency as reported elsewhere (Fig 2A) [15,16]. The RiboRid method using 108 ArOP sets (S1 Table) removed rRNAs from the E. coli total RNA samples at a level down to 1.63%, which is comparable to commercial methods at a fraction of the cost (Fig 2A and Table 1). Furthermore, the RiboRid method successfully depleted coverage from the entire 16S and 23S rRNA genes (Fig 2B). Technical mock samples that underwent the same RiboRid process with nuclease-free water instead of the Hybridase enzyme, showed virtually no difference to untreated control. In addition, we inspected RNA--Seq profile of highly expressed genes to examine enzymatic reaction introduced biases or RNA degradation as reported earlier [17,18], and did not detect any sign of RNA degradation or bias (S2 Fig). In addition, we did not observe non-specific enrichment of short rRNA reads that has been reported in the previous enzyme-based rRNA depletion approach (S3 Fig) [17]. Next, we analyzed gene expression by RNA-Seq that are performed with different rRNA removal methods. Except for the untreated and mock samples, all the rRNA-depleted samples and replicates were linearly correlated, even across the different methods (S4 Fig). The RNA--Seq dataset prepared from two different biological replicates using the RiboRid method showed high reproducibility, with a correlation of 0.977 (Pearson's R 2 ).
The number of genes observed by RNA-Seq was also significantly increased by rRNA removal. From the RNA-Seq result of Ribo-Zero, RiboErase, or RiboRid treated samples, expressions of 4,201, 4,166, and 4,122 genes were observed, respectively, while only 2,738 genes were detected in the control samples.
Furthermore, by the nature of the Poisson process of RNA-Seq, gene expression levels measured by multiple replicates showed heteroscedasticity; that is, the variance of gene expression measured between replicates increases as the expression level decreases [19,20]. The variation of biological replicates is a combination of technical and biological variation [19]. Given the experimental setup, the difference in observed variations across different methods were originated from different technical variations, since the libraries were prepared from aliquots of the same biologically replicated samples. Thus, we expected to observe different distribution of variation across expression level if rRNA removal method was unsuccessful or introducing bias. The variation of gene expressions from the samples prepared without rRNA depletion was much higher than that of the rRNA-depleted samples and was not heteroscedastic (S5 Fig). In contrast, samples prepared using Ribo-Zero, RiboErase, and RiboRid had identical distribution of experimental variations. Moreover, the average coefficient of variation (standard deviation divided by mean expression of a gene, which implies that the variance of measured expression adjusted by expression level) of the control samples was more than 3.6-fold higher than that of the rRNA-depleted sample, while the three methods had similar values. Considering the number of genes detected and the distribution of experimental variances, rRNA depletion greatly reduced the counting error of RNA-Seq, especially for genes with low abundance. Next, we tested RNA input up to 1 μg and different ratios of ArOPs to total RNA. There was virtually no difference in the rRNA content of sequenced reads, indicating that the rRNA depletion method works for various amounts of input RNA and ratios of ArOPs to RNA ( Fig  2C). In addition, correlations between these technically replicated samples were higher than 0.989 (S6A Fig). This shows that the method is sufficiently robust to accommodate practical variations in biological experiments. To examine the possible nonspecific binding of ArOPs to mRNA, we repeated RNA-Seq from the same RNA sample using the RiboRid method at 45˚C (Fig 2C). Although the RiboRid method was able to remove rRNA, the fraction of reads from rRNA increased to 27.4% because of the suboptimal reaction temperature of Hybridase. In addition, correlations of the samples with other technical replicates were relatively low (Pearson's R 2 = 0.744 ± 0.019) (S6B Fig). More importantly, we predicted 46 nonspecific bindings of ArOPs to mRNA (E-value < 1 using the Basic Local Alignment Search Tool (BLAST) (S2 Table). Eighteen of the probes had consecutive matches longer than 15 bp. Expression of the 18 genes measured by the RNA-Seq library prepared from the rRNA-subtracted RNA sample using the RiboRid method at 45˚C was significantly underestimated when compared to the normal RiboRid-treated sample (p = 0.027; Wilcoxon's rank-sum test) (Fig 2D). In contrast, we did not observe any significant difference in gene expression between untreated, mock, and RiboRid treated samples (S7 Fig). This illustrates the importance of a high reaction temperature in preventing nonspecific binding of probes to mRNA and demonstrates that high reaction temperature of RiboRid was able to avoid experimental bias.
Next, we applied the method and the same ArOPs to RNA-Seq of K. oxytoca KCTC1686, a bacterial species closely related to E. coli. Owing to the sequence similarity between the two bacteria, 76 out of 108 anti-E. coli ArOPs matched the rRNAs of K. oxytoca (Fig 3). Although gaps in K. oxytoca rRNAs, where no ArOP hybridizes, span a maximum length of 163 nt, 76 ArOPs were dense enough to remove rRNAs to a level less than 3% of mappable sequence reads in K. oxytoca RNA-Seq, with a correlation of 0.998 (Pearson's R 2 ) between biological replicates ( Table 1). rRNA removal of K. oxytoca using anti-E. coli ArOPs indicates that the number of ArOPs may be reduced to a much lower density. Thus, we designed a new probe set in which each probe was spaced at a longer distance. The new probe set consisted of 55 probes that could hybridize approximately every 80 nt on the rRNAs (S3 Table). We also adjusted the protocol further to reduce the cost and time of the RiboRid method. First, ArOPs were synthesized with a C3 spacer at the 3 0 end, which is not utilized as a substrate for ligases and polymerase in downstream library construction. Originally, oligonucleotide probes remained after the RiboRid reaction needed to be removed by DNase digestion, since they can interfere with the downstream PCR amplification step of sequencing library construction. However, C3 spacermodified oligonucleotides cannot be used as primers for DNA polymerase; thus, DNase I digestion of probes at the end of the RiboRid method can be avoided. Alternatively, the 3 0 end of pre-existing probes can be blocked by attaching dideoxynucleotide analogs using terminal deoxynucleotidyl transferases instead of C3 spacer modification. Then, the column-based clean-up step was replaced by SPRI paramagnetic bead purification. Together with these  modifications, the cost and processing time of the RiboRid method can be reduced considerably (S1 and S2 Protocol) without increasing the rRNA fraction or experimental bias ( Table 1 and S6C Fig).

Depletion of rRNA from the RNA samples from various bacterial species
To examine the versatility of the RiboRid method for diverse bacterial species, we applied RiboRid to two Gram-negative bacteria, Pseudomonas aeruginosa and Bacteroides thetaiotaomicron, and two Gram-positive bacteria, Eubacterium limosum and Staphylococcus aureus, with ArOPs specific for each bacterium ( Table 2). P. aeruginosa is a Gram-negative opportunistic pathogen that is a common cause of pneumonia and other infections in various parts of the body [21]. rRNAs were successfully removed from the laboratory derivative PAO1 isolate [22] using the standard RiboRid method to a level where rRNA comprised an average of 4% of mappable reads in RNA-Seq. However, we encountered an interesting profile of rRNA operons (Fig 4A). There were large amounts of RNA fragments in the intergenic region of the rRNA operons. Considering rRNA maturation and processing in bacteria, the fragments were predicted to be the pre-rRNAs that were not depleted due to the absence of complementary ArOPs. The pre-rRNAs comprised an average of 0.84% of the mapped reads in RNA-Seq. Although the amount of pre-rRNA reads was negligible when compared to mRNA reads, we tested two different approaches to further remove pre-rRNA species. First, a two-fold increase in Hybridase in the reaction did not affect the overall amount of non-mRNA or pre-rRNA reads (Fig 4A and 4B). In fact, supplementation of a few anti-pre-rRNA oligos (S4 Table) in the Hybridase reaction reduced the amount of pre-rRNA fragments from the RNA sample to 0.06% (Fig 4B).
We further applied the RiboRid method to two Gram-positive bacteria: S. aureus, a common pathogen that is one of the biggest threats to the health care system due to the emergence of strains with antibiotic resistance [23,24]; and E. limosum, a model strain of CO 2 -fixing acetogenic bacteria [25]. rRNAs from both strains were efficiently removed by the RiboRid method ( Table 2). The fraction of S. aureus rRNAs in RNA-Seq was lower than 1.72% regardless of the culture medium, with Pearson's correlation (R 2 ) between biological replicates higher than 0.955 (S8 Fig). In particular, rRNAs comprised only 0.02% of the total mapped reads when RNA was prepared from the cells grown in RPMI medium. The RiboRid method was also capable of reducing rRNAs from E. limosum down to 2.25% with high reproducibility (R 2 = 0.949) between biological replicates ( Table 2). B. thetaiotaomicron is one of the major constituents of the human gut commensal microbiome [26,27]. Unlike in other bacterial species, the RiboRid method left relatively high amount of rRNA from RNA sample of B. thetaiotaomicron. In detail, rRNA from the RiboRidtreated sample comprised 27.20 ± 5.03% ( Table 2), although the RNA expression measured from RiboRid-treated sample was highly correlated with that of Ribo-Zero-treated samples with a Pearson's correlation (R 2 ) of 0.915 ± 0.01 (S8 Fig). To improve the efficiency of rRNA removal, we designed a new ArOP set (ArOP Set 2) with different DNA sequences that have higher melting temperatures, so that they can hybridize to rRNA stronger at the elevated reaction temperature of RiboRid. Although the fraction of rRNA to the total mapped reads was reduced to 17.82 ± 0.05%, its correlation with the Ribo-Zero-or standard ArOP (ArOP Set 1)treated sample was decreased to 0.741 (R 2 ) (S8 Fig). Lowering the reaction temperature to 58˚C, which is at least 10˚C lower than the annealing temperature, did not improve the rRNA removal efficiency ( Table 2). A recent report also indicated that both Ribo-Zero and RNase Hbased removal of rRNAs Bacteroidetes dorei, were relatively inefficient, possibly due to its high genomic AT-content [13]. As illustrated by the duplex-specific nuclease-based rRNA removal method [16], further experimental and probe optimization may be required for organisms with high AT content, such as Bacteroidetes; however, the RiboRid method was able to remove rRNAs to a level less than 20% of total RNA.

Application of the RiboRid method to Term-Seq
To further validate the RiboRid method for other challenging transcriptomic analysis methods, we conducted Term-Seq [28].  cell at the nucleotide level [28]. Thus, it is sensitive to RNA degradation, which generates false signals. We successfully detected 1,308 transcript 3 0 ends across the transcriptome of the model actinomycetes, S. coelicolor, which is a Gram-positive soil bacterium well-known for producing various secondary metabolites [29]. Reads mapped on rRNAs comprised only 0.36% of the reads mapped to the genome ( Table 3). To address possible non-specific random RNA degradation during the RiboRid treatment that may decrease resolution of single nucleotide-sensitive study such as Term-Seq, we inspected profile of Term-Seq (S9 Fig). For example, transcript 3 0 end detected from two biological replicates coincides with each other in nucleotide-level (S9A and S9B Fig). Furthermore, the Term-Seq signal from two replicates were highly correlated (Pearson's R 2 of 0.9998) across the entire genome (S9C Fig). Thus, we concluded that RiboRid did not introduce any noise for the transcriptomic study. With successful depletion of rRNAs from total RNA of S. coelicolor with high GC content (72.4%), we were able to identify a conserved rho-independent transcriptional terminator motif at the 110 transcript 3 0 ends, which resembles the rho-independent transcriptional terminators of the closely related actinomycetes, S. lividans (Fig 5A) [30]. In addition, the location of transcription termination occurred on the T-rich tract of the conserved motif, which was consistent with the previous report (Fig 5B) [30]. This demonstrates that the RiboRid method can be applied to transcriptome sequencing other than RNA-Seq.

Application of the RiboRid method to ribosome profiling
Ribosome profiling is a transcriptome sequencing technique that surveys mRNAs that are actively translated by capturing ribosome-protected fragments (RPFs) [31]. It is one of the most challenging transcriptome analysis techniques because it generates highly fragmented short rRNA fragments during enzymatic digestion of RNA that are not protected by ribosomes. Depletion of highly fragmented short rRNA fragments is a major technical challenge in bacterial ribosome profiling. We performed ribosome profiling of E. coli using the RiboRid method, which was modified by adopting a previously developed streamlined protocol (see Methods) [32]. Initially, the method failed to remove rRNA fragments from the RPFs, such that 95.6% of sequencing reads were mapped to rRNA ( Table 3). Close inspection of the sequencing profile indicated that specific regions of rRNA were enriched (Fig 6A). Thus, we designed and supplemented six additional probes (additional set 1) targeting the enriched regions (S5 Table). The specific RNA fragments were clearly depleted after RiboRid treatment with the additional probe set 1 (Fig 6B). However, another rRNA fragment remained and became more noticeable relative to the fragments targeted by probe set 1. We constructed six additional probes (additional set 2) and RiboRid with the two additional probe sets reduced the fraction of reads mapped on rRNA down to the level shown by the Ribo-Zero-treated sample, which was approximately 70% (Fig 6C and Table 3). Meta-analysis of the ribosome profile on coding sequences showed a nucleotide-resolution profile of ribosome footprint and pausing on the start codon (Fig 6D), which is comparable to Ribo-Zero treated profile (S10 Fig), illustrating the capability of RiboRid in ribosome profiling [33]. Taken together, RiboRid effectively removed the fraction of rRNA even from highly fragmented polysomal RNA samples, indicating that it can be used for complicated transcriptomic studies other than RNA-Seq. The single base-pair resolution required for probing 3 0 end information and the precise polysomal location on the RNA were not compromised by the method developed in this study.

Discussion
Based on the transcriptomic sequencing results presented in this study, we demonstrated a strategy for effectively depleting rRNA from bacterial RNA samples. The RiboRid method is highly advanced and effective method, capable of removing rRNA from highly fragmented RNA samples without losing precise biological information of the mRNA. Owing to thorough experimental and protocol design considerations regarding hybridization temperature and probe design, the prevalent rRNA species in bacterial RNA samples could be efficiently removed without introducing any experimental bias or causing nonspecific mRNA removal. Furthermore, a short incubation time at denaturing temperature, low concentration of divalent magnesium cations, and mild pH prevented RNA degradation observed in the previous enzymatic approaches [17,18]. More importantly, the high hybridization and reaction temperature prevented nonspecific binding of probes to mRNA that may occur in the previous approaches using low hybridization temperature [13]. Optimization of the probe number could streamline the protocol such that the cost is dramatically reduced ($10 per reaction) when compared to the commercially available method (Ribo-Zero;~$80 per reaction), with comparable or better rRNA removal efficiency. Although this method is based on enzymatic RNA digestion, the total experiment time is only 1 h with a hands-on time of 30 min; this time is comparable to the hybridization-magnetic subtraction-based method. RiboRid is highly efficient and could remove rRNA with an efficiency of up to 99.99% in the S. aureus sample. However, in the case of B. thetaiotaomicron, rRNA comprised 17% of all mappable rRNA reads. This is possibly because of the high AT content of the organism, in  which hybridization of the probe at high temperature may be inefficient. However, this can be easily compensated for by using a slightly higher sequencing depth, as the sequencing cost is much lower than the difference between the costs of the RiboRid and commercial methods.
To facilitate the use of this method, a step-by-step protocol is available in the S1 Protocol. The protocol contains all the consumables, materials, and equipment required to perform the RiboRid reaction. In addition, an in-house Python script is freely available through a public repository (https://github.com/SBRG/RiboRid_Design); this script can design a probe set for any custom genome or rRNA sequence. In addition, we designed and deposited probe sets for representative bacterial species available in the RefSeq database (n = 5,467). This method provides a cost-effective, rapid, and powerful alternative means to deplete rRNA that outperforms previously developed and reported methods. In particular, this method is valuable for routine large-scale transcriptome studies and reduces the burden of high-cost commercial kits.

RNA extraction
Total RNA of E. coli, K. oxytoca, and B. thetaiotaomicron were extracted using the RNASnap method [34] with slight modifications. Briefly, cell pellets collected from 5 mL mid-exponentially grown culture (OD 600nm = 0.4) was resuspended in 100 μL RNASnap solution (18 mM EDTA, 0.025% SDS, 1% β-mercaptoethanol, and 95% formamide) and incubated at 95˚C for 7 min. The resuspension was centrifuged at 16,000 × g for 5 min and RNA with a size larger than 200 nt in the clear supernatant was purified using RNA Clean & Concentrator-5 Kit (Zymo, Irvine, CA, USA) according to the manufacturer's instructions. Total RNA of E. limosum was extracted from 100 mL mid-exponential growth culture (OD 600nm = 1.5). The collected cells were resuspended in 500 μL lysis buffer (20 mM Tris-HCl [pH 7.4], 140 mM NaCl, 5 mM MgCl 2 , and 1% Triton X-100) and ground using a mortar and pestle after flash freezing with liquid nitrogen. RNA was isolated from the supernatant of the ground sample using TRIzol Reagent (Thermo, Waltham, MA, USA) followed by an RNA Clean & Concentrator-5 Kit according to the manufacturer's instructions. Total RNA of P. aeruginosa and S. aureus was extracted using Quick-RNA Fungal/Bacterial Microprep Kit (Zymo) from 3 mL culture sampled at an OD 600nm of 0.4. Briefly, the cell pellet was resuspended in 800 μL RNA lysis buffer. The resuspension was then transferred into a ZR BashingBead Lysis Tube (0.1 and 0.5 mm) and homogenized in a bead beater. RNA in the cleared lysate (400 μL) was purified by column purification. Total RNA of S. coelicolor was extracted using the hot phenol method from 50 mL culture samples at the early exponential, transition, late exponential, and stationary growth phases (OD 600nm = 0.6, 2, 3.5, and 5, respectively). First, the collected cells were resuspended in Solution 1 (25 mM Tris-HCl [pH 8.0], 10 mM EDTA, 50 mM glucose, and 2 mg/mL lysozyme) and incubated for 10 min at 30˚C. The supernatant of the mixture was removed by centrifugation, and the cell pellet was resuspended in Solution 2 (50 mM sodium acetate [pH 5.2], 10 mM EDTA, and 1% sodium dodecyl sulfate). The suspension was mixed with an equal volume of phenol:chloroform (5:1) solution and incubated for 5 min at 65˚C. RNA was isolated by isopropanol precipitation from the aqueous phase and separated by centrifugation.

Anti-rRNA oligonucleotide probes (ArOPs) design
The probes used were 32 nt-long deoxynucleotides reverse complementary to the rRNA with melting temperatures at least 3˚C higher than the rRNA digestion reaction temperature (68˚C). The anti-rRNA oligonucleotide probes were designed using an in-house Python program (https://github.com/SBRG/RiboRid_Design). The program can be either imported into other pipelines or executed as a standalone program from the command line. As a basic input, the script takes in GenBank files (with annotated rRNA) or a FASTA file containing the rRNA sequences for the target organism. Given the rRNA sequences in either format, the script starts by building a consensus sequences for rRNAs. These consensus sequences were then used to build a DNA probe library. To build the library for each rRNA type, we started from the 50 th sequence position of the rDNA. The first probe was defined as the DNA sequence that started from this position and was the length of the user-defined probe. With a default probe length of 32, the probe starts at position 50 and ends at position 81. If the probe has a melting temperature above the defined melting threshold (68˚C used in this study), then the script looks for the next probe downstream with the user-defined maximum gap (default of 50 nt) between the probes. Typically, 12-25 probes per kb of rRNA are required for efficient removal of rRNA (one probe per 40-80 nt). The melting temperatures of the probes were calculated using the OligoAnalyzer Tool (IDT, Coralville, IA, USA). If the melting temperature of the probe is below the threshold, then the script will look for another probe up to the user-defined maximum search space (default of 10 bp) upstream of the starting search position and record the probe with the highest melting temperature. Note that the probe with the highest melting temperature was recorded even if the melting temperature of the probe was below the threshold. Following these steps, the script generates probes in a stepwise manner until all the rRNAs are covered. Once the design process is finished, the designed probe sequences are recorded in a FASTA file, and the metadata associated with each probe (containing information such as its melting temperature, start position, end position, etc.) is recorded in a.csv text file. Using the script, the melting temperatures of 32 nt-long probes were higher than 68˚C in 98.3% of the design attempts. In a small number of cases (1.7%), the melting temperature criteria were met by manually adjusting the lengths of probes from 28 to 43 nt (especially in the AT-rich pre-rRNA region).

RiboRid
Before performing rRNA depletion, RNA content was measured using a Qubit RNA HS Assay Kit (Thermo) with a Qubit 2.0 fluorometer (Thermo). RNA extract was subjected to DNase I treatment for removal of DNA contamination at 37˚C for 10 min in 15 μL reaction mixture containing 0.5-1 μg total RNA, 1.5 μL 10× DNase I Buffer, and 2 U RNase-free DNase I. Then, DNase I was inactivated at 75˚C for 10 min after addition of 15 μL Hybridase Complement Buffer comprising 90 mM Tris-HCl (pH 7.5) and 200 mM KCl. The ArOP mixture (1 μL) containing 5 pmol of each probe and 100 mM MgCl 2 was added to the DNase I-treated RNA sample. To hybridize the ArOP to RNA, the temperature of the RNA-ArOP mixture was heated to 90˚C for 1 s and cooled to 65˚C on a thermocycler. When the temperature of the mixture reached 65˚C, 10 U Hybridase Thermostable RNase H (Lucigen, Middleton, WI, USA) prewarmed at room temperature was added. The reaction was carried out at 65˚C for 20 min and 90˚C for 1 s to rehybridize the ArOP to the remaining rRNA, and then at 65˚C for an additional 10 min. Then, RNA larger than 200 nt was extracted using RNA Clean & Concentrator Kit. A nucleic acid-binding column was saved for later use. The remaining ArOPs were removed by DNase I treatment in a reaction composed of 10 U RNase-free DNase I and 5 μL 10× DNase I buffer in a total reaction volume of 50 μL. The DNase I reaction was carried out by consecutive 5 min incubations at 25,30,35,40, and 45˚C. The reaction was again purified with the RNA Clean & Concentrator Kit using a nucleic acid binding column saved from the previous clean-up. A step-by-step protocol is provided in S1 and S2 Protocol. The columnbased RNA purification step can be replaced by a solid-phase reversible immobilization (SPRI) bead-based purification method. In this study, RNA was purified using 1.8 volume of CleanNGS DNA & RNA SPRI Bead (Bulldog Bio, Portsmouth, NH, USA) according to the manufacturer's instructions. When using a 3 0 -C3 carbon spacer-modified ArOP that does not participate in downstream reactions (reverse transcription, adaptor ligation, and PCR), the DNase I treatment step for removing residual ArOPs after the hybridase reaction was not performed.

RNA-Seq
RNA-Seq libraries were constructed from approximately 100 ng rRNA-depleted RNA using a TruSeq Stranded mRNA LT Sample Prep Kit (Illumina, San Diego, CA, USA) or KAPA RNA HyperPrep Kit (Roche, Basel, Switzerland) according to the manufacturer's instructions. Constructed sequencing libraries were quantified using a Qubit dsDNA HS Assay Kit with a Qubit 2.0 fluorometer and a TapeStation 2200 (Agilent, Santa Clara, CA, USA) equipped with a High Sensitivity D1000 Screen Tape (Agilent). The library was sequenced using an Illumina platform. Information on NGS (run recipe and instrument) is summarized in S6 Table. Term-Seq Total RNA extract containing 5 μg total RNA was treated with DNase I at 37˚C for 15 min in a 50 μL reaction comprising 2 U RNase-free DNase I and 5 μL 10× DNase I buffer. The Term-Seq library was constructed as previously described [29]. Briefly, 5 0 -DNA adaptor was ligated to 1 μg of the pooled RNA in 25 μL reaction mixture containing 150 pmol amino-blocked 3 0 -DNA adaptor (5 0 -p-NNAGATCGGAAGAGCGTCGTGT-AmMO), 25 U T4 RNA Ligase 1 (NEB, Ipswich, MA, USA), 2.5 μL 10× T4 RNA Ligase 1 Buffer, 2.5 μL 10 mM ATP, 2 μL dimethyl sulfoxide, and 9.5 μL polyethylene glycol 8000. The reaction mixture was incubated for 2.5 h at 23˚C followed by RNA purification using 2.2× volumes of AMPure XP Beads (Agencourt, Beverly, MA, USA). The adaptor-ligated RNA was then subjected to the columnbased RiboRid method as described above. Then, rRNA-depleted RNA was fragmented by incubating at 72˚C for 90 s with 1 μL RNA Fragmentation Reagent in a total reaction volume of 10 μL, followed by RNA purification using 2.2× volumes of AMPure XP Beads. Then, cDNA was synthesized from rRNA-depleted RNA using 200 U SuperScript III Reverse Transcriptase and 10 pmol amino-blocked reverse transcription primer (5-TCTACACTCTTTCCC TACACGACGCTCTTC) in a total reaction volume of 30 μL. The synthesized cDNA was purified with 1.5× volumes of AMPure XP Beads. Then, cDNA 3-adaptor (5-p-NNAGATCGGAA GAGCACACGTCTGAACTCCAGTCAC-AmMO) was ligated into the 3-DNA adaptor ligation mixture for 8 h at 23˚C. The ligation product was purified using 1.8× volumes of AMPure XP Beads and amplified using Phusion High-Fidelity DNA Polymerase (Thermo). The library was sequenced using an Illumina platform. Information on NGS (run recipe and instrument) is summarized in S6 Table. Ribosome profiling using Ribo-Zero Ribosome profiling was conducted using a method described in a previous report [32]. Briefly, 50 mL E. coli culture was collected after 5 min of treatment with chloramphenicol (34 mg/mL) at an exponential growth phase. Cells were flash frozen with 0.5 mL lysis buffer (1% Triton X-100, 34 μg/mL chloramphenicol, 133 mM NaCl, 4.75 mM MgCl 2 , and 19 mM Tris-HCl, pH 7.5) and lysed using a mortar and pestle. Then, the supernatant containing 100 μg RNA was treated with 2,000 gel units of micrococcal nuclease (MNase; NEB). Polysomes were recovered from MNase-digested samples using Illustra MicroSpin S-400 HR columns (GE Healthcare, Chicago, IL, USA) followed by phenol:chloroform:isoamyl alcohol extraction. rRNA was removed from 1 μg polysome-protected RNA using a Ribo-Zero rRNA Removal Kit according to the manufacturer's instructions. The rRNA-subtracted RNA samples were phosphorylated by treating 10 U T4 polynucleotide kinase (NEB) at 37˚C for 1 h and purified with RNeasy MinElute columns (Qiagen, Hilden, Germany). Sequencing libraries were prepared from phosphorylated RNA samples using the NEBNext Small RNA Library Prep Set for Illumina (NEB) according to the manufacturer's instructions. The library was sequenced using an Illumina platform. Information on NGS (run recipe and instrument) is summarized in S6 Table. Ribosome profiling using RiboRid Polysomal RNA was prepared as described previously. The polysomal RNA samples were subjected to phosphorylation without rRNA removal by treating with 10 U T4 polynucleotide kinase at 37˚C for 1 h and purified with RNeasy MinElute columns. 5 0 sequencing adaptors were attached to phosphorylated RNA samples using the NEBNext Small RNA Library Prep Set for Illumina according to the manufacturer's instructions. After adaptor ligation, RNA samples were subjected to a column-based RiboRid reaction with dideoxy-modified ArOP and without DNase I treatment. RNA purification steps were performed using the RNA Clean & Concentrator Kit for RNA size > 17 nt. The dideoxy-modified ArOPs were prepared by incubating 5.4 nmol (50 pmol each probes) of ArOPs at 37˚C for 4 h with 50 U terminal deoxynucleotide transferase (TdT; Thermo), 20 μL 5× Reaction Buffer, 20 mM ddNTP (as a mixture of 5 mM each ddATP, ddTTP, ddGTP, and ddCTP) in 100 μL reaction mixture. The modified ArOP was purified using the Oligo Clean & Concentrator Kit (Zymo) as described by the manufacturer. Then, 3 0 RNA adaptor ligation, reverse transcription, and sequencing library amplification were performed using the NEBNext Small RNA Library Prep Set for Illumina. The library was sequenced using an Illumina platform. Information on NGS (run recipe and instrument) is summarized in S6 Table. Data processing Sequencing data were processed on a CLC Genomics Workbench (CLC Bio, Aarhus, Denmark). Raw reads were trimmed using the Trim Sequence Tool in NGS Core Tools with a quality limit of 0.05. Reads with more than two ambiguous nucleotides were discarded and the quality trimmed reads were mapped on either the reference genome or rRNAs using the following parameters: mismatch cost of 2, indel cost of 3, and similarity and length fraction of 0.9. Reference genome sequences were downloaded from NCBI under the accession numbers NC_000913.3 (E. coli), NC_016612 (K. oxytoca), NC_002516.2 (P. aeruginosa), NC_010079 (S. aureus DNA), NC_010063 (S. aureus plasmid pUSA300HOUMR), NC_012417 (S. aureus plasmid pUSA01-HOU), NC_004663 (B. thetaiotaomicron chromosomal DNA), NC_004703.1 (B. thetaiotaomicron plasmid DNA), NZ_CP019962.1 (E. limosum), and NC_003888.3 (S. coelicolor). Gene expression levels were calculated by counting strand-specific reads of genes (antisense) and normalized using DESeq2 [20]. Nonspecific ArOP binding was predicted by standalone BLASTN alignments of probes to genomic sequences lacking rRNA sequences. For motif analysis, DNA sequences 40 nt upstream to 20 nt downstream of the 3 0 ends of the transcript were analyzed using MEME software (v5.0.4) using the mode of zero or one site per sequence [35]. The found motif had an E-value of 4.2 × 10 −55 and each sequence aligned to the motif with a p-value lower than 1.32 × 10 −4 . The meta-analysis of ribosome profiles was performed using the STATR pipeline, as described previously [36]. Wilcoxon's rank-sum test was performed via statistical analysis using the SciPy package [37] Supporting information S1 Protocol. Column-based RiboRid. Genes with different degree of sequence similarity (E-value; BLASTN) were compared as a group. E<1; �15 nt: genes with E-value lower than 1 (BLASTN alignment with ArOP) and with 15 nt or more consecutive matches. ns: difference between fold-change distributions are not significant (Wilcoxon's rank-sum test). Box limits, whiskers, and center lines indicate 1 st and 3 rd quartiles, 10 th and 90 th percentiles, and the median of the distribution, respectively. Dots are individual genes. Number of subjected genes are indicated above the graph. ns: not significant (Wilcoxon's rank-sum test).