Bacterial culture conditions

Bacteroides strains were routinely cultured in an anaerobic chamber (Coy Laboratory Products) with an anaerobic gas mix (85% N2, 10% CO2 and 5% H2) at 37 °C. Routine cultivation of all strains was performed in TYG medium and on Brain Heart Infusion Supplemented (BHIS) plates. For a detailed description of media composition and culture conditions for the RNA-seq analysis, refer to Supplementary Table 1.

Growth assays in the presence of diverse carbon sources were carried out in minimal medium supplemented with 0.5% of a suitable carbon source, as follows. A single colony of wild-type B. thetaiotaomicron VPI-5482 (AWS-001) was inoculated into 5 ml minimal medium–glucose and incubated anaerobically for 24 h. Then, 1 ml of this culture was centrifuged (2,000g for 3 min) to pellet bacterial cells that were resuspended in an equal volume of minimal medium (without a carbon source). This was subsequently used to inoculate (1:100 dilution) minimal medium containing an appropriate carbon source and incubated for the indicated time, following which aliquots (optical density equivalents of ~4) were collected for RNA extraction.

Stress response assays were performed in TYG medium as indicated below. A single colony of AWS-001 was inoculated into 5 ml TYG medium and incubated anaerobically overnight. The next day, it was sub-cultured into TYG medium and grown to the mid-exponential phase (~7 h; OD600 = 2.0). This culture was sub-divided into 5 ml fractions corresponding to each stress condition and centrifuged to pellet the bacterial cells, as before. The pellet was resuspended in an equal volume of TYG medium containing the indicated concentration of a stressor and incubated for a further 2 h, following which samples were collected for RNA extraction.

Bacterial genetics

A detailed list of the strains, plasmids and oligonucleotides used in this study can be found in Supplementary Table 10. To create ∆masB, we employed a previously established method48. To this end, we assembled 1-kilobase sequences flanking the deletion site into the pExchange-tdk suicide vector and introduced this construct into E. coli S17-1 λpir. The resulting transformants were mated with B. thetaiotaomicron Δtdk (AWS-003) and the resulting conjugants were selected on 5-fluoro-2′-deoxyuridine plates. Single recombinants were isolated on BHIS agar with 200 μg ml−1 gentamicin and 25 μg ml−1 erythromycin (BHISgent/erm). Double recombinants, leading to either scarless deletion mutants or wild-type revertants, were identified by their ability to grow on BHIS agar with 200 μg ml−1 5-fluoro-2′-deoxyuridine while being unable to grow on BHIS agar with 25 μg ml−1 erythromycin. The masB complementation strain (masB+) was assembled using a version of the pNBU2 vector system, as previously described48. The complete masB gene sequence was integrated into pNBU2 using Gibson Assembly (New England Biolabs (NEB)) to ensure transcription from the native TSS. This construct was conjugated into the ∆masB strain via E. coli S17-1 λpir, as described above.

The other deletion mutants (∆cur, ∆1675 and ∆cur∆masB) were generated using the pSIE1 plasmid system49. Briefly, 750-nucleotide flanking regions around the deletion site were Gibson assembled (NEB) into the linearized pSIE1 plasmid (SpeI and BamHI digested). The assembled plasmid was subsequently introduced into B. thetaiotaomicron via conjugation with E. coli S17-1 λpir and the conjugants were streaked onto BHISgent/erm plates. Resistant colonies were cultured overnight in TYG medium without antibiotics, and dilutions (10−1 to −3) were plated onto BHIS agar with 100 ng ml−1 anhydrotetracycline (aTC). Colony PCR and Sanger sequencing were used to confirm the intended deletions. The cur+ complementation strain resulted from Gibson Assembly (NEB) of full-length cur with pWW3452 (AWP-015) such that transcription was under the control of the phage promoter on the plasmid. It was ensured that the 3′ end of the transcript maintained a reading frame with the downstream FLAG and His tags. The construct was conjugated into the ∆cur background as described above.

Total RNA purification and removal of genomic DNA

All of the RNA-seq samples were collected as biological duplicates. Total RNA was isolated by the hot-phenol method, as follows. Briefly, bacterial cultures containing a total of ~4 OD equivalents of cells were collected and a one-fifth volume of stop mix was added (5% vol vol−1 water-saturated phenol; pH > 7.0; 95% vol vol−1 ethanol)50. Cell lysis was achieved by incubation with lysozyme (600 µl; 0.5 mg ml−1) and sodium dodecyl sulfate (60 µl of a 10% solution) for 2 min at 64 °C with the subsequent addition of NaOAc (66 µl of a 3 M solution). Extraction with 750 μl phenol (ROTIAqua-Phenol) was carried out at 64 °C for 6 min, followed by the addition of 750 μl chloroform. Precipitation of RNA from the aqueous phase was achieved with twice the volume of ethanol and 3 M NaOAc (30:1) mix and incubated at −80 °C overnight. The samples were then centrifuged and the pellets washed with ethanol (75% vol vol−1), followed by resuspension in 50 µl RNase-free water. Traces of genomic DNA were removed by treating ~40 µg total RNA with 5 U DNase I (Fermentas) and 0.5 µl SUPERase·In RNase Inhibitor (Ambion) in a reaction volume of 50 µl. Samples for dRNA-seq were prepared by pooling equimolar amounts (each 100 ng) of total RNA from each condition.

cDNA library preparation and sequencing

For dRNA-seq, samples were treated according to a previous protocol16. Before synthesizing cDNA, pooled total RNA was fragmented by ultrasound (four pulses of 30 s each at 4 °C) and then treated with T4 Polynucleotide Kinase (NEB). The RNA sample was then split evenly, one half of which was treated with Terminator Exonuclease to enrich for primary transcripts, whereas the other half remained untreated. The samples were then poly(A)-tailed using poly(A) polymerase and the 5′-PPP removed using 5′ polyphosphatase (Epicentre Biotechnologies). RNA adaptors were ligated and the synthesis of first-strand cDNA was performed using M-MLV reverse-transcriptase and oligo(dT) primers. The cDNA was subsequently amplified to a concentration of ~10–20 ng µl−1, purified (Agencourt AMPure XP kit; Beckman Coulter Genomics) and fractionated in a size range of 200–600 bp. The libraries were deep-sequenced on an Illumina NextSeq 500 system using 75 bp read length at Vertis Biotechnologie.

For conventional RNA-seq, samples were first depleted of ribosomal RNA (rRNA) using the Pan-Prokaryote riboPOOL kit (siTOOLs Biotech). This involved incubation of 1 µg total RNA with 100 pmol rRNA-specific biotinylated DNA probes at 68 °C for 10 min, followed by a shift to 37 °C for 30 min in 0.25 mM ethylenediaminetetraacetic acid (EDTA), 2.5 mM Tris-HCl (pH 7.5) and 500 mM NaCl. Depletion of rRNA–DNA hybrids was achieved by two 15-min incubation periods with streptavidin-coated magnetic Dynabeads MyOne Streptavidin C1 beads (0.45 mg; Thermo Fisher Scientific) in 0.25 mM EDTA, 2.5 mM Tris-HCl (pH 7.5) and 1 M NaCl at 37 °C. The samples were then purified using the Zymo RNA Clean & Concentrator kit along with DNase I treatment (Zymo Research). Libraries were prepared with the NEBNext Multiplex Small RNA Library Prep kit for Illumina (NEB) according to the manufacturer’s instructions and the following modifications. Samples were fragmented at 94 °C for 2.75 min, per the NEBNext Magnesium RNA Fragmentation Module (NEB) with subsequent RNA purification using the Zymo RNA Clean & Concentrator kit. The fragmented RNA was then 3′ dephosphorylated, 5′ phosphorylated and decapped with 10 U T4 Polynucleotide Kinase ± 40 nmol ATP and 5 U RNA 5′ Pyrophosphohydrolase (NEB). After each step, RNA was purified as mentioned above. The fragmented RNA was then ligated to adapters (3′ SR and 5′ SR, pre-diluted 1:3 in nuclease-free water) and the cDNA was amplified for 14 cycles. These barcoded libraries were purified using MagSi-NGSPREP Plus beads (AMSBIO) at a 1.8:1 ratio of beads to sample volume. Libraries were checked for quantity and quality using a Qubit 3.0 Fluometer (Thermo Fisher Scientific) and a 2100 Bioanalyzer with the High Sensitivity DNA kit (Agilent). Pooled libraries were sequenced on the NextSeq 500 platform (Illumina) at the Core Unit SysMed of the University of Würzburg.

Read processing and mapping

Generated reads were quality checked using FastQC (version 0.11.8) and adapters were trimmed using Cutadapt (version 1.16) with Python (version 3.6.6), using the following parameters: -j 6 -a Illumina Read 1 adapter=AAGATCGGAAGAGCACACGTCTGAACTCCAGTCA -a Poly A=AAAAAAAAAAA --output=out1.fq.gz --error-rate=0.1 --times=1 --overlap=3 --minimum-length=20 --nextseq-trim=20 3_1. For both sequencing data types (dRNA-seq and conventional RNA-seq), READemption51 (version 0.4.5) was used to map reads to the B. thetaiotaomicron VPI-5482 reference genome (NC_004663.1) and plasmid (NC_004703.1). Details of the alignment statistics can be found in Supplementary Table 11.

Transcriptome annotation

We used the ANNOgesic pipeline27 (version 0.7.33) to update the annotations of TSSs, terminators, operons and noncoding RNAs, as previously described16. In short, TSSs were identified using the TSSpredator function, which compares the relative enrichment of reads between Terminator Exonuclease-treated and untreated libraries to identify enriched peaks that are characteristic of the protected 5′ ends of primary transcripts. TSSs were classified on the basis of this enrichment and distance, relative to a coding gene. Primary TSSs were identified as having the highest enrichment within a 300-bp region upstream of an open reading frame. All other TSSs within this region were classified as secondary TSSs. Internal TSSs were defined as originating on the sense strand within a coding sequence, whereas antisense TSSs were those that originated on the antisense strand and overlapped with or within 100-bp flanks of a sense gene. All remaining TSSs were classified as orphan TSSs. To predict terminators, ANNOgesic utilizes two heuristic algorithms, one of which scans the genome for Rho-independent terminators (TransTermHP) and the other predicts terminators by detecting a decrease in read coverage between two adjacent genes. Leveraging the wealth of our diverse conditional datasets, we additionally utilized the operon detection function of ANNOgesic (default settings) to predict both operons and sub-operons based on the other detected features; namely, TSSs, transcripts and genes supplied as general feature format (GFF) files.

The automated prediction of noncoding RNAs was done using the srna function of ANNOgesic, which first compares predicted transcripts with known RNAs within the sRNA database (all sequences were downloaded from BSRD52) and the non-redundant protein database (ftp://ftp.ncbi.nih.gov/blast/db/FASTA/). Candidates that were contained in the sRNA database were retained, while those contained in the non-redundant protein database were excluded from further analysis. All remaining transcripts were classified as intergenic sRNA candidates if they possessed a defined TSS, stable secondary structure (RNAfold normalized folding energy < −0.05) and length within 30–500 nucleotides and did not overlap with any other genetic feature in either sense or antisense orientation. The prediction of cis-asRNAs was based on the same criteria along with the presence of an annotated gene on the opposing strand. Untranslated region (UTR)-derived sRNAs were classified as 5′ if they shared the TSS with an mRNA and were associated with a read coverage drop or processing site in front of the coding sequence. Similarly, 3′ UTR-derived sRNAs were predicted by either a TSS or processing site in the 3′ UTR and either a processing site or terminator shared with an mRNA. Intra-operonic sRNAs were associated with a TSS or processing site at the 5′ end and a read coverage drop or processing site at the 3′ end.

As with all computationally automated predictions, manual curation was necessary to ensure the accuracy of the global annotations. Manual confirmation of sRNA predictions was guided by the following criteria for excluding putative sRNA annotations: (1) lack of an identifiable promoter or processing site up to 50 bp upstream of a predicted sRNA’s 5′ end; (2) complete overlap with an annotated mRNA in the most recent genome update on the National Center for Biotechnology Information’s Nucleotide database (NC_004663.1; 21 August 2022); (3) complete overlap with an annotated terminator sequence; (4) overlap with cis-regulatory elements such as riboswitches or RNA thermometers; and (5) no evident change in read coverage relative to flanking regions. As a result, we report a total of 135 intergenic sRNAs in B. thetaiotaomicron. This includes an overlap of 91 sRNAs previously identified16, as well as 44 novel sRNA candidates. In total, 20 candidates previously annotated as intergenic sRNAs16 were reassigned here to different sub-classes (Supplementary Table 12): 11 were re-annotated as asRNAs (as a divergently encoded genetic feature was discovered in the present data), three were re-annotated as 3′ UTR-derived sRNAs (due to their overlap with the 3′ end of an mRNA) and six were re-annotated as 5′ UTR-derived sRNAs (as the present dataset allowed us to refine TSS annotations). Additionally, 14 previously annotated sRNA candidates16 were eliminated as they were very weakly expressed in the former dataset and their existence was not further supported by the present dataset, despite a greater sequencing depth.

Updated annotations of the TSSs, terminators, operons and noncoding RNAs can be accessed via Theta-Base 2.0 (http://micromix.helmholtz-hiri.de/bacteroides/). Coding genes and sRNAs can be interrogated using either their ID (BT_xxxx or BTncxxx) or—when applicable—their trivial name (for example, MasB).

Prediction of invertible DNA regions

Invertible DNA regions were predicted using the PhaseFinder -locate pipeline as described in ref. 28. Briefly, inverted repeats were determined by allowing no mismatches for repeats of a maximum of 11 bp, one mismatch for repeats up to 13 bp and two mismatches for repeats with lengths exceeding 19 bp. Homopolymeric inverted repeats were removed and the maximum GC content per inverted repeat was filtered to be between 15 and 85%.

Differential gene expression analysis

Differential gene expression analysis was performed using the R package edgeR (version 3.38.2)53,54. Using the filterByExpr function, the genes with a counts per million (CPM) value of >0.6635 (equivalent to around ten reads per sample) across all replicates under each growth condition (median library size of ~15 million reads) were retained for differential analysis. While calling for contrasts (using the makeContrast function), the analysis was sub-divided into two groups based on the respective control condition: all conditions with varying carbon sources (including starvation) were compared with glucose as a control, whereas all stress conditions were compared with the condition immediately before stress induction (that is, a mid-exponential phase culture in TYG medium). Differential gene expression data across all conditions relative to their respective control condition can be found in Supplementary Table 4.

Gene set annotation and enrichment analyses

We assembled a list of functionally annotated gene sets from the literature. We recovered annotations of PULs 01–96 from the Polysaccharide Utilization Loci Database4; CPSs and conjugative transposons from ref. 55; the genes transcribed from promoter motifs PM1 and PM2 from ref. 16; regulons from RegPrecise version 3.2 (ref. 56); the Cur regulon from ref. 13; annotated KEGG pathways and modules from the KEGG database (accessed on 1 December 2022); Gene Ontology terms from UniProt (accessed on 25 November 2021); and predicted KEGG modules and pathways and Gene Ontology terms from an eggNOG version 5.0 (ref. 57) annotation of the B. thetaiotaomicron genome. Gene set enrichment analysis was performed with the fgsea R package over all gene sets with more than nine genes, except for PULs, which were retained irrespective of their gene number. Genes were ranked based on the −log10[P value] × sign[FC] metric.

Northern blot

Northern blotting was performed as described previously16. In short, total RNA (2.5–10 µg) was electrophoretically resolved on a 6% (vol vol−1) polyacrylamide (PAA) gel containing 7 M urea and electro-blotted onto a membrane (Amersham Hybond-XL) at 50 V and 4 °C for 1 h. The blots were probed with gene-specific 32P-labelled oligonucleotides in Hybri-Quick buffer (Carl Roth) at 42 °C and subsequently exposed to a phosphor screen as required. Images were visualized using a phosphorimager (FLA-3000 Series; Fuji).

Reanalysis of TIS data

We reanalysed fitness data from a comprehensive B. thetaiotaomicron transposon mutant library that probed a suite of hundreds of different conditions, including 48 different carbon sources and 56 stress-inducing compounds19, in the context of our extensive noncoding RNA annotation (see above). This was done with the primary objective of identifying and possibly correlating gene expression from our transcriptomic dataset with mutant fitness data and thereby allowing us to draw biologically meaningful conclusions. To further streamline our analysis, we focused exclusively on independently encoded intergenic sRNAs since phenotypes pertaining to such mutants would probably not involve polar effects. Consequently, of the 135 intergenic sRNAs identified in B. thetaiotaomicron, we obtained fitness data for 81 sRNAs (~70%), of which 28 were associated with statistically significant effects (|t| > 4) in at least one successful experiment (Supplementary Table 8). A successful experiment requires that each gene is represented by a sufficient number of barcode counts58. The fitness of a gene is defined as the average log2[change in relative abundance of its mutants (|fit|)]. Negative and positive values mean that the sRNA mutants were less or more fit, respectively, than the average strain in the pool. Experiments with ‘jackpot’ effects, whereby the disruption of an sRNA resulted in a large competitive advantage versus the other mutants in the pool, were retained, but specifically labelled as strong phenotypes (|fit| > 2 and |t| > 5) (Supplementary Table 8). A third category, namely ‘combined’, comprised those phenotypes that were both strong and significant, per the above criteria.

Launch of Theta-Base 2.0

Theta-Base 2.0 (http://micromix.helmholtz-hiri.de/bacteroides/) was created using Micromix (https://github.com/BarquistLab/Micromix)59, which relies on Flask60 (back end) and Vue.js (front end), storing underlying visualization and expression data using MongoDB. The Clustergrammer plugin uses the API from the Ma’ayan laboratory61, while the heat map plugin follows the same front- and back-end architecture as the main site. Gene set annotations (Gene Ontology terms, KEGG pathways and modules, PULs, CPSs, conjugative transposons, promoter motifs and known regulons) were prepared as described in the section ‘Gene set annotation and enrichment analyses’ and can be found in Supplementary Table 3. The sRNA fitness dataset was adapted from Supplementary Table 8. Deployment of the back and front ends uses Gunicorn (https://readthedocs.org/projects/gunicorn-docs) and Nginx62.

B. thetaiotaomicron datasets can be manually selected by first clicking the Bacteroides Theta tab, followed by providing a title and then selecting an appropriate dataset using the dropdown menu. Users can select from a choice of expression data (that is, normalized in CPM or log2[FC] (compared with control conditions)) or between the entire dataset or specific sRNA fitness data. As an option, columns of interest can be further customized using the ‘Select columns’ box and by subsequently clicking the ‘Add’ tab. Users also have the option to add their own data, as outlined in additional tabs, such as by uploading a delimited file.

The resulting data tables are displayed in the browser and can be filtered or transformed. For example, the ‘Filter’ button allows data tables to be filtered using keywords with prompts to make the search process seamless. The ‘Functional annotation’ button permits the user to select from a large number of preset manually curated gene lists, such as ‘GO term’, ‘KEGG pathway’, ‘PUL’, ‘CPS’ and ‘CTn’, to name a few. A third button, ‘ncRNA’ permits selection of manually curated noncoding RNA gene sets (for example, ‘High-confidence intergenic sRNAs’, ‘Intergenic sRNAs’ and ‘Cis-antisense RNAs’). Once the desired genes have been filtered, they may be transformed by clicking the ‘Transform data’ button and performing operations such as rounding values, log conversion or calculating transcripts per million. Once datasets have been loaded by the user, they can be further examined using three visualization modes; namely, ‘Heatmap’, ‘Clustergrammer’ and ‘JBrowse’. Two- and three-dimensional heat maps can be generated using the Heatmap function. Note that the heat map defaults to a three-dimensional option, but users can manually switch to the two-dimensional option. Heat maps can be customized using the menu on the left that permits changes to the colour gradients and overall structure. Customized heat maps can be downloaded in SVG or PNG formats using the download tab. Alternatively, for clustering according to genes or conditions, ‘Clustergrammer’ is recommended. Selecting this tab generates a two-dimensional dynamic heat map of the data that can be further investigated using the menu on the left. Currently, the tool only permits a maximum of 200 rows to be loaded and users will be notified if more rows are selected. Customized heat maps and data tables can be downloaded using the ‘Take snapshot’ and ‘Download matrix’ buttons, respectively. For a detailed view of normalized coverage plots for the investigated conditions, in addition to those published in the first iteration of Theta-Base16, users can select the ‘JBrowse’ button63. Users are free to select from a range of updated annotations displaying high-resolution maps for noncoding RNAs (ncRNA), TSSs (TSSv3), terminators (term_v2), operons (Operon_structure), a transposon insertion map related to the fitness data (Tn_insertions) and invertons (Inverted_repeats).

On the top right of the website there are four buttons. The ‘Padlock’ button locks the current state of the site, allowing users to copy their URL and share with colleagues. The next (‘Download’) button allows users to download the currently selected dataset as an Excel or a delimited file (such as .csv). The ‘New document’ button will re-load the website so users can select another dataset. The ‘Help’ button—when clicked—will provide pop-over text explaining various features of the site.

Antibiotics growth curve analyses and agar strip assays

Bacterial growth curves were determined by inoculating a single colony each of AWS-003 (Δtdk; referred to as the wild-type in Fig. 4) and AWS-029 (∆masB) into 5 ml TYG medium and incubating overnight under anaerobic conditions. These cultures were sub-cultured (1:100 dilution) in 2 ml TYG medium containing the indicated final concentrations of doxycycline (Sigma-Aldrich), tetracycline (AppliChem) or a water control. The samples (200 µl volume) were incubated in a 96-well flat-bottom plate (Nunclon) at 37 °C (doxycycline) or 40 °C (tetracycline) with continuous shaking (double orbital) in a microplate spectrophotometer (BioTek Epoch 2). Optical densities were recorded every 20 min. The assay was performed in three biological replicates, each comprising technical duplicates.

Antibiotics strip assays were performed by dipping a sterile cotton swab into overnight TYG cultures of AWS-003 or AWS-029 and streaking on BHIS agar plates containing strips for doxycycline (EM103 (HiMedia; in Fig. 4c) or 92156 (Liofilchem; in Fig. 5g)) or tetracycline (EM056; HiMedia). The plates were incubated anaerobically for 48 h at 37 °C and images were taken. The minimal inhibitory concentrations were derived from the positions where the inhibition ellipses intersected the strips.

Gene co-expression analyses

Correlation of the expression of all B. thetaiotaomicron genes across all of the profiled carbon source and stress conditions was calculated by generating a correlation matrix (Pearson’s correlation score) of the z scores of the CPM values of each gene. To identify the correlation in expression between our gene sets (see ‘Gene set annotation and enrichment analyses’) and a given gene of interest (MasB in Fig. 5a and BT_1675 in Extended Data Fig. 10f), the median of the correlation values between all genes within a gene set and the gene of interest was calculated. Gene sets composed of fewer than ten operons were excluded from this analysis.

MS2 affinity purification and sequencing

A B. thetaiotaomicron ΔmasB mutant complemented with either MS2-MasB (AWS-062) or untagged MasB (AWS-036) was diluted 1:100 in TYG medium from an overnight culture grown anaerobically at 37 °C. At an OD600 of 2.0, expression of MS2-tagged MasB and untagged MasB was induced by the addition of 200 ng ml−1 aTC. After another 2 h of growth at 37 °C, 90 OD equivalents of the cultures were collected, centrifuged for 20 min at 2,000g and 4 °C and snap-frozen in liquid nitrogen. MS2 pulldown and RNA purification was performed as described in ref. 64, but with slight modifications to adapt the protocol to Bacteroides. Specifically, the column was washed only six (instead of eight) times with buffer A before elution. Elution itself was then induced with 600 µl (rather than 300 µl) elution buffer.

For library preparation (at Vertis Biotechnologie), the RNA samples were first fragmented using ultrasound (four pulses of 30 s each at 4 °C). Then, an oligonucleotide adapter was ligated to the 3′ ends of the RNA molecules. First-strand cDNA synthesis was performed using M-MLV reverse-transcriptase and the 3′ adapter as a primer. The first-strand cDNA was purified and the 5′ Illumina TruSeq sequencing adapter was ligated to the 3′ end of the antisense cDNA. The resulting cDNA was PCR-amplified to ~10–20 ng μl−1 using a high-fidelity DNA polymerase. The cDNA was purified using the Agencourt AMPure XP kit (Beckman Coulter Genomics) and analysed by capillary electrophoresis. For Illumina sequencing, the samples were pooled in approximately equimolar amounts. To deplete sequences derived from 5S rRNA, the cDNA pool was digested using probes specific for bacterial 5S and Cas9 endonuclease. Afterwards, the cDNA pool was fractionated in the size range of 200–600 bp using a preparative agarose gel. An aliquot of the size-fractionated pool was analysed by capillary electrophoresis. The cDNA pool was sequenced on an Illumina NextSeq 500 system using a read length of 2 × 150 bp.

Generated reads were quality-checked using FastQC (version 0.11.8) and adapters were trimmed using BBDuk with the following parameters: qtrim=r trimq=10 ktrim=r ref=bbmap/ressources/adapters.fa k=23 mink=11 hdist=1 tpe tbo. BBmap was used to map reads to the B. thetaiotaomicron VPI-5482 reference genome (NC_004663.1) and plasmid (NC_004703.1), as well as to the MS2-MasB sequence. Read quantification was performed using featureCounts (2.0.1). Differential abundance analysis between the MS2-MasB and untagged samples was conducted using edgeR65 in combination with RUVSeq66 to estimate the factor of unwanted variation using replicate sample with correction factor k=1.

IntaRNA prediction of sRNA–mRNA interactions

In silico interaction prediction between MasB and its putative mRNA targets fusA2 and BT_1675 was performed with the help of IntaRNA40,41 using the Vienna RNA package (2.4.14 and boost 1.7) at default settings along with the output flag (--out=pMinE:FILE.csv) to generate minimal energy values for intermolecular index pairs. For visualization, the resulting values were plotted in form of a heat map in R (version 4.2).

In vitro transcription and radiolabelling of RNA

DNA templates for in vitro transcription were amplified using genomic DNA and primer pairs carrying a T7 promoter (Supplementary Table 10). The in vitro transcription reaction was performed using the MEGAscript T7 kit (Thermo Fisher Scientific) followed by DNase I digestion (1 U; 37 °C; 15 min). RNA products were then excised from a 6% (vol vol−1) PAA-7M urea gel by comparison with a Low Range RNA ladder (Thermo Fisher Scientific) and eluted overnight in elution buffer (0.1 M NaOAc, 0.1% sodium dodecyl sulfate and 10 mM EDTA) on a thermoblock at 8 °C and 1,400 r.p.m. The next day, the RNA was precipitated in an ethanol:NaOAc (30:1) mix, washed with 75% ethanol and resuspended in 20 µl water (at 65 °C for 5 min).

Radioactive labelling of the in vitro-transcribed RNA was carried out by dephosphorylating 50 pmol RNA with 25 U calf intestine alkaline phosphatase (NEB) in a 50 µl reaction and incubating at 37 °C for 1 h. The dephosphorylated RNA was extracted using phenol:cholorform:isoamylalcohol (25:24:1) and precipitated as described above. Next, 20 pmol of this RNA was 5′ end-labelled (20 µCi 32P-γATP) using 1 U polynucleotide kinase (NEB) at 37 °C for 1 h in a 20 µl reaction. The labelled RNA was purified using a G-50 column (GE Healthcare) and extracted from a PAA gel as described above.

EMSA

EMSA was carried out in a reaction volume of 10 μl, containing 1× RNA structure buffer (Ambion), 1 μg yeast RNA (~4 μM final concentration), 5′ end-labelled MasB RNA (4 nM final concentration) and an mRNA segment of 137 nucleotides in length, spanning the predicted MasB target site within BT_1675 (see Fig. 5e) at final concentrations of 0, 8, 16, 32, 64, 128, 256, 512 and 1,024 nM. Following incubation at 37 °C for 1 h, 3 μl of 5× native loading dye (0.2% bromophenol blue, 0.5× TBE and 50% glycerol) was added to each tube. All of the samples were loaded on a native 6% (vol vol−1) PAA gel in 0.5× TBE buffer and run at 300 V and 4 °C for 3 h. The gel was dried, exposed and visualized using a phosphorimager (FLA‐3000 Series; Fuji). The experiment was repeated three times and quantified using ImageJ version 1.52s67 and GraphPad Prism version 9 for Windows (GraphPad Software; www.graphpad.com). The dissociation constant (Kd) was calculated via the one site-specific binding formula:

$$Y=B_}} \times X/(}}}+X)$$

where Y is the specific binding; X the concentration of radio ligand; Bmax the maximum binding in the same unit as Y; and Kd the dissociation constant in the same unit as X.

qRT-PCR analysis

For the qRT-PCR assays, Δtdk B. thetaiotaomicron (AWS-003; referred to as the wild-type in Fig. 5), ΔmasB (AWS-029) and masB+ (AWS-036) were grown anaerobically overnight at 37 °C in 5 ml TYG medium, then sub-cultured 1:100 in TYG and induced with 200 ng ml−1 aTC. Around four optical density equivalents of samples were collected at the early exponential phase (OD600 = 0.3), mid-exponential phase (OD600 = 2.0) and late exponential phase (OD600 = 3.7) for RNA extraction, as described above. For the starvation condition, the same strains were grown anaerobically for 24 h at 37 °C in 5 ml minimal medium supplemented with 0.5% glucose, and then sub-cultured 1:100 in 0.5% glucose-containing minimal medium supplemented with 200 ng ml−1 aTC. At OD600 = 2.0, the cultures were centrifuged, the supernatant was discarded and the pellet was resuspended in minimal media without a carbon source and incubated anaerobically at 37 °C for another 2 h. Around four optical density equivalents of the samples were collected for RNA extraction. qRT-PCR reactions were performed as described in ref. 16. A minimum of three biological replicates were pipetted and plates were analysed on a QuantStudio 5 instrument (Thermo Fisher Scientific).

Dual-plasmid fluorescence reporter assay

Strains of E. coli TOP10, which were engineered to carry translational fusions of superfolder GFP43 to different variants of the 5′ part of the BT_1675 coding sequence, were cultured in Lysogeny Broth medium supplemented with chloramphenicol (20 μg ml−1) and carbenicillin (100 μg ml−1) until an OD600 of 0.5 was reached. Subsequently, 100 μl of the cultures was collected and subjected to three washes with 1× phosphate-buffered saline before fixation with a 4% paraformaldehyde solution. The fluorescence intensity of GFP was measured in phosphate-buffered saline using flow cytometry (NovoCyte Quanteon; Agilent).

Statistics and reproducibility

Conventional RNA-seq of diverse growth conditions, dRNA-seq of pooled conditions and MAPS were performed in biological duplicates. Testing for differential expression or enrichment was performed using the generalized linear model likelihood ratio test implemented in edgeR53. Northern blots were performed in two biological replicates. EMSAs were performed in technical triplicates. qRT-PCR analysis was performed in a minimum of three biological replicates, each comprising technical duplicates, and a Mann–Whitney test was used to call significant comparisons. Growth curve experiments were performed in a minimum of three biological replicates, unless explicitly mentioned otherwise. Two-plasmid GFP reporter assays were performed in biological triplicates and Tukey’s multiple comparisons test was used to test for statistically significant differences. Minimal inhibitory concentration strip assays were performed in a minimum of five biological replicates for doxycycline and three biological replicates for tetracycline. No statistical method was used to predetermine the sample size. Instead, sample sizes were chosen based on previous experience and studies16. No data were excluded from the analyses. The experiments were not randomized. The investigators were not blinded to allocation during the experiments and outcome assessment.

Reporting summary

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