Mycobacterium smegmatis strain mc²155 was grown in Middlebrook 7H9 supplemented with ADC (Albumin Dextrose Catalase, final concentrations 5 g/L bovine serum albumin fraction V, 2 g/L dextrose, 0.85 g/L sodium chloride, and 3 mg/L catalase), 0.2% glycerol and 0.05% Tween 80. For the exponential phase experiment (Dataset 1), 50 ml conical tubes containing 5 ml of 7H9 were inoculated with M. smegmatis to have an initial OD = 0.01. Cultures were grown at 37°C and 200 rpm. Once cultures reached an OD of 0.7–0.8, they were frozen in liquid nitrogen and stored at -80°C until RNA purification. For hypoxia experiments (Dataset 2), a protocol similar to the Wayne model (Wayne and Hayes, 1996) was implemented. Briefly, 60 ml serum bottles (Wheaton, product number 223746, actual volume to top of rim 73 ml) were inoculated with 36.5 ml of M. smegmatis culture with an initial OD = 0.01. The bottles were sealed with rubber caps (Wheaton, W224100-181 Stopper, 20 mm) and aluminum caps (Wheaton, 20 mm aluminum seal) and cultures were grown at 37°C and 125 rpm to generate hypoxic conditions. Samples were taken at an early stage of oxygen depletion when growth had slowed but not completely stopped (15 h) and at a later stage when a methylene blue indicator dye was fully decolorized and growth had ceased (24 h). These time points were experimentally determined according to growth curve experiments (see Supplementary Figure S1). 15 ml of each culture were sampled and frozen immediately in liquid nitrogen until RNA extraction.

RNA was extracted as follows: frozen cultures stored at -80°C were thawed on ice and centrifuged at 4,000 rpm for 5 min at 4°C. The pellets were resuspended in 1 ml Trizol (Life Technologies) and placed in tubes containing Lysing Matrix B (MP Bio). Cells were lysed by bead-beating twice for 40 s at 9 m/sec in a FastPrep 5G instrument (MP Bio). 300 μl chloroform was added and samples were centrifuged for 15 min at 4,000 rpm at 4°C. The aqueous phase was collected and RNA was purified using Direct-Zol RNA miniprep kit (Zymo) according to the manufacturer’s instructions. Samples were then treated with DNase Turbo (Ambion) for 1 h and purified with an RNA Clean & Concentrator-25 kit (Zymo) according to the manufacturer’s instructions. RNA integrity was checked on 1% agarose gels and concentrations were determined using a Nanodrop instrument. Prior to library construction, 5 μg RNA was used for rRNA depletion using Ribo-Zero rRNA Removal Kit (Illumina) according to the manufacturer’s instructions.

After rRNA depletion, RNA samples from each biological replicate were split in three, in order to generate two 5′-end differentially treated libraries and one RNAseq expression library (next section). RNA for library 1 (“converted” library) was treated either with RNA 5′ pyrophosphohydrolase RPPH (NEB) (exponential phase experiment, Dataset 1), or with 5′ polyphosphatase (Epicenter) (hypoxia experiment, Dataset 2), in order to remove the native 5′ triphosphates of primary transcripts, whereas RNA for Library 2 (“non-converted” library) was subject to mock treatment. Thus, the converted libraries capture both 5′ triphosphates (converted to monophosphates) and native 5′ monophosphate transcripts, while non-converted libraries capture only native 5′ monophosphates (see scheme in Supplementary Figure S2A). Library construction was performed as described by Shell et al. (2015a). A detailed scheme showing the workflow of 5′-end libraries construction, the primers and adapters used in each step, and modifications to the protocol are shown in Supplementary Figure S2B.

One third of each rRNA-depleted RNA sample was used to construct RNAseq expression libraries. KAPA stranded RNA-Seq library preparation kit and NEBNext Ultra RNA library prep kit for Illumina (NEB) were used for Dataset 1 and Dataset 2, respectively, according to manufacturer’s instructions. The following major modifications were introduced into the protocols: (i) For RNA fragmentation, in order to obtain fragments around 300 nt, RNA was mixed with the corresponding buffer and placed at 85°C for 6 min (Dataset 1), or at 94°C for 12 min (Dataset 2). (ii) For library amplification, 10 or 19–23 PCR cycles were used for Dataset 1 and Dataset 2, respectively. The number of cycles was chosen according to the amount of cDNA obtained for each sample. After purification, DNA concentration was measured in a Qubit instrument before sequencing.

For 5′-end-mapping libraries from Dataset 1, Illumina MiSeq paired-end sequencing producing 100 nt reads was used. For 5′ end directed libraries from Dataset 2 as well as for all expression libraries, Illumina HiSeq 2000 paired-end sequencing producing 50 nt reads was used. Sequencing was performed at the UMass Medical School Deep Sequencing Core Facility. Quality of the generated fastq files was checked using FastQC.

Paired-end reads generated from 5′-end-directed libraries were mapped to the M. smegmatis mc²155 NC_008596 reference genome. In order to reduce noise from the imprecision of transcriptional initiation, only the coordinate with the highest coverage in each 5 nt window was used for downstream analyses. For read filtering, different criteria were used for the 2 datasets according to the library depth and quality (see Supplementary Figure S3). In order to discriminate between TSSs and CSs, the ratio of the coverage in converted/non-converted libraries for each detected 5′ end was calculated. To focus our analyses on the 5′ends that are relatively abundant in their local genomic context, we employed a filter based on the ratio of 5′ end coverage to expression library coverage in the preceding 100 nt. 5′ ends for which this ratio was ≤0.05 were removed. After this filter, 15,720 5′ ends remained and were further analyzed using a Gaussian mixture modeling to differentiate TSSs from CSs with a high confidence in Dataset 1 (Figure 1A). For this analysis, we used the iterative expectation maximization (EM) algorithm in the mixtools package (Benaglia et al., 2009) for R (version 1.1.0) to fit the mixture distributions.

Reads were aligned to the Mycobacterium smegmatis str. mc²155 reference genome (Accession number NC_008596) using Burrows-Wheeler Aligner (Li and Durbin, 2009). For comparison of gene expression levels according to presence or absence of Shine-Dalgarno sequences, RPKMs were calculated for all genes. The DEseq2 pipeline was used to evaluate the changes in gene expression in hypoxia (Love et al., 2014).

For analysis in Figure 1D, TSSs were classified as follows: those coordinates located ≤500 bp upstream from an annotated gene were considered to be primary TSSs (pTSS). Coordinates located within an annotated gene were classified as internal (iTSS) or N-associated internal TSSs (N-iTSSs) if they were located within the first 25% of the annotated coding sequence. N-iTSSs were considered for reannotation as a pTSSs only if their associated gene lacked a pTSS. TSSs located on the antisense strand of a coding sequence, 5′ UTR, or 3′ UTR were considered as antisense (aTSS). 5′ UTRs boundaries were assigned after assignment of pTSSs to genes annotated in the mc²155 reference genome (accession number NC_008596). When a gene had more than one pTSS, the longest of the possible 5′ UTRs was used for assignment of aTSSs. In the case of genes for which we did not identify a pTSS, we considered a hypothetical leader sequence of 50 bp for assignment of aTSSs. For assignment of aTSSs in 3′ UTRs, we arbitrarily considered a sequence of 50 bp downstream the stop codon of the gene to be the 3′ UTR. Finally, TSSs not belonging to any of the above-mentioned categories were classified as orphan (oTSSs).

Adjacent genes with the same orientation were considered to be co-transcribed if there were at least 5 spanning reads between the upstream and the downstream gene in at least one of the replicates in the expression libraries from Dataset 1. After this filtering, a downstream gene was excluded from the operon if: (1) it had a TSS ≤ 500 bp upstream the annotated start codon on the same strand, and/or (2) had a TSS within the first 25% of the gene on the same strand, and/or (3) the upstream gene had a TSS within the last 50–100% of the coding sequence. Finally, the operon was assigned only if the first gene had a primary TSS with a confidence ≥95% according to the Gaussian mixture modeling.

For CS categorization in Figure 4D, we established stringent criteria in order to determine the frequency of CSs in each location category relative to the amount of the genome comprising that location category. For 3′ UTR regions, we considered only CSs that were located between 2 convergent genes. To assess frequency relative to the whole genome, we considered the sum of all regions located between two convergent genes. For 5′ UTRs we considered all CSs located between 2 divergent genes, and the sum of all leader lengths for genes having a pTSS whose upstream gene is in the opposite strand (divergent) determined in this study was used for assessing relative frequency. For 5′ ends corresponding to cleavages between co-transcribed genes we used the operon structures determined in this study, and the sum of all their intergenic regions was used for assessing relative frequency. Finally, for CSs located within coding sequences all genes were considered, as all of them produced reads in the expression libraries. The sum of all coding sequences in NC_008596 genome was used for assessing relative frequency, after subtracting overlapping regions to avoid redundancy.

For validation of TSSs and CSs, RNA samples from M. smegmatis were split in two and treated with or without RPPH (NEB) in order to remove the native 5′ triphosphates of primary transcripts or not, respectively. Then, an adapter oligo SSS1016 (CTGGAGCACGAGGACACTGACATGGACTGAAGGAGTrArGrArArA, where nts preceded by “r” are ribonucleotides and the rest of the oligo is composed of deoxyribonucleotides) was ligated to the RNA 5′ ends using T4 RNA ligase (NEB). Prior to ligation, 8 μl of RNA sample were combined with 1 μl of 1 μg/μl adapter oligo and incubated at 65°C for 10 min. For ligation, the 9 μl of RNA-oligo mix were combined with: 10 μl 50% PEG8000, 3 μl 10X ligase buffer, 3 μl 10 mM ATP, 3 μl DMSO, 1 μl Murine RNase inhibitor (NEB), and 1 μl T4 ligase (NEB). Ligation reactions were incubated at 20°C overnight and then cleaned using RNA Clean and Concentrator 25 kit (Zymo). Both RPPH-treated and mock-treated samples were used for cDNA synthesis. Reactions in absence of reverse transcriptase were performed to control for genomic DNA contamination. For amplification of specific 5′ ends, PCR was done using a forward primer SSS1017 binding to the adapter oligo (CTGGAGCACGAGGACACTGA) and a reverse (specific) primer binding near the predicted 5′ end (see Supplementary Table S1). For PCRs, a touchdown protocol in which the annealing temperature was reduced 1°C every cycle was performed as follows: (i) initial step of DNA denaturation at 95°C for 5 min, (ii) 17 cycles of 95°C for 30 s, 72–55°C (touchdown) for 20 s and 68°C for 25 s, (iii) 20 cycles of 95°C for 30 s, 55°C for 20 s and 68°C for 25 s and (iv) a final elongation step at 68°C for 5 min. Each amplified fragment was sequenced using the specific primer. A TSS or CS was validated if (i) the 5′ end position coincided with that mapped the 5′ end libraries and (ii) the PCR product was more abundant in the RPPH than in the no RPPH treatment (TSS) or the PCR product was equally abundant in the RPPH and in the no RPPH treatment (CS).

For validation of the MSMEG_0063 promoter, an M. smegmatis mutant strain lacking the region comprising the genes MSMEG_0062-MSMEG_0066 was transformed with either of the 3 following constructs: (i) Wt promoter, which has the gene MSMEG_0063 with the native predicted promoter region and the downstream genes MSMEG_0064-MSMEG_0066, (ii) Δpromoter, in which the predicted promoter region for MSMEG_0063 was deleted, and (iii) Mutated promoter, in which two point mutations were introduced in the predicted -10 region of the MSMEG_0063 promoter. These constructs were inserted in the L5 site of the M. smegmatis genome.

In order to study the transcriptome structure of M. smegmatis, RNAs from triplicate cultures in exponential phase were used to construct 5′ end mapping libraries (Dataset 1) according to our previously published methodologies (Shell et al., 2015a,b) with minor modifications. Briefly, our approach relies on comparison of adapter ligation frequency in a dephosphorylated (converted) library and an untreated (non-converted) library for each sample. The converted libraries capture both 5′ triphosphate and native 5′ monophosphate-bearing transcripts, while the non-converted libraries capture only native 5′ monophosphate-bearing transcripts (Supplementary Figure S2). Thus, assessing the ratios of read counts in the converted/non-converted libraries permits discrimination between 5′ triphosphate ends (primary transcripts from transcription start sites) and 5′ monophosphate ends (cleavage sites). By employing a Gaussian mixture modeling analysis (Figure 1A) we were able to identify 5,552 TSSs in M. smegmatis with an observed probability of being a TSS ≥ 0.95 (high confidence TSSs, Supplementary Table S2). A second filtering method allowed us to obtain 222 additional TSSs from Dataset 1 (Supplementary Figure S3). A total of 5,774 TSSs were therefore obtained from Dataset 1. In addition, data from separate libraries constructed as controls for the hypoxia experiment (Dataset 2) in “The Transcriptional Landscape Changes in Response to Oxygen Limitation’ were also included in this analysis to obtain TSSs. After noise filtering (Supplementary Figure S3), 4,736 TSSs from Dataset 2 were identified. The union of the two datasets yielded a total of 6,090 non-redundant high confidence TSSs, of which 4,420 were detected in both datasets (Supplementary Figure S4 and Supplementary Table S2).

Although not all 5′ ends could be classified with the Gaussian mixture modeling, we were able to assign 57% of the 5′ ends in Dataset 1 to one of the two 5′ end populations with high confidence (5,552 TSSs and 3,344 CSs). To validate the reliability of the Gaussian mixture modeling used to classify 5′ ends, we performed two additional analyses. First, according to previous findings in Mtb (Cortes et al., 2013) and other well studied bacteria (Sass et al., 2015; Berger et al., 2016; Čuklina et al., 2016; D’arrigo et al., 2016), we anticipated that TSSs should be enriched for the presence of the ANNNT -10 promoter consensus motif in the region upstream. Evaluation of the presence of appropriately-spaced ANNNT sequences revealed that 5′ ends with higher probabilities of being TSSs are enriched for this motif, whereas for those 5′ ends with low probabilities of being TSSs (and thus high probabilities of being CSs) have ANNNT frequencies similar to that of the M. smegmatis genome as a whole (Figure 1B). Secondly, we predicted that TSSs should show enrichment for A and G nts at the +1 position, given the reported preference for bacterial RNA polymerases to initiate transcription with these nts (Lewis and Adhya, 2004; Mendoza-Vargas et al., 2009; Mitschke et al., 2011; Cortes et al., 2013; Shell et al., 2015b; Thomason et al., 2015; Berger et al., 2016). Thus, we analyzed the base enrichment in the +1 position for the 5′ ends according to the p-value in the Gaussian mixture modeling (Figure 1C). These results show a clear increase in the percentage of G and A bases in the position +1 as the probability of being a TSS increases, while the percentage of sequences having a C at +1 increases as the probability of being a TSS decreases. These two analyses show marked differences in the sequence contexts of TSSs and CSs and further validate the method used for categorization of 5′ ends.

To study the genome architecture of M. smegmatis, the 6,090 TSSs were categorized according to their genomic context (Figure 1D,E and Supplementary Table S2). TSSs located ≤500 nt upstream of an annotated gene start codon in the M. smegmatis str. mc²155 (accession NC_008596) reference genome were classified as primary TSSs (pTSS). TSSs within annotated genes on the sense strand were denoted as internal (iTSS). When an iTSS was located in the first quarter of an annotated gene, it was sub-classified as N-terminal associated TSS (N-iTSS), and was further examined to determine if it should be considered a primary TSS (see below). TSSs located on the antisense strand either within a gene or within a 5′ UTR or 3′ UTR were grouped as antisense TSSs (aTSSs). Finally, TSSs located in non-coding regions that did not meet the criteria for any of the above categories were classified as orphan (oTSSs). When a pTSS also met the criteria for classification in another category, it was considered to be pTSS for the purposes of downstream analyses. A total of 4,054 distinct TSSs met the criteria to be classified as pTSSs for genes transcribed in exponential phase. These pTSSs were assigned to 3,043 downstream genes, representing 44% of the total annotated genes (Supplementary Table S3). This number is lower than the total number of genes expressed in exponential phase, in large part due to the existence of polycistronic transcripts (see operon prediction below). Interestingly, 706 (23%) of these genes have at least two pTSSs and 209 (7%) have three or more, indicating that transcription initiation from multiple promoters is common in M. smegmatis. We used 5′ RACE to confirm seven selected pTSSs (Supplementary Table S1), all of which mapped to the same position by both methods. Four of these were novel TSSs not reported by Li et al. (2017).

A total of 995 iTSSs (excluding the iTSSs that were also classified as a pTSS of a downstream gene, see Supplementary Figure S5 for classification workflow) were identified in 804 (12%) of the annotated genes, indicating that transcription initiation within coding sequences is common in M. smegmatis. iTSSs are often considered to be pTSSs of downstream genes, to be spurious events yielding truncated transcripts, or to be consequences of incorrect gene start annotations. However, there is evidence supporting the hypothesis that iTSSs are functional and highly conserved among closely related bacteria (Shao et al., 2014), highlighting their potential importance in gene expression.

We were also able to detect antisense transcription in 12.5% of the M. smegmatis genes. Antisense transcription plays a role in modulation of gene expression by controlling transcription, RNA stability, and translation (Morita et al., 2005; Kawano et al., 2007; Andre et al., 2008; Fozo et al., 2008; Giangrossi et al., 2010) and has been found to occur at different rates across bacterial genera, ranging from 1.3% of genes in Staphylococcus aureus to up to 46% of genes in Helicobacter pylori (Beaume et al., 2010; Sharma et al., 2010). Of the 1,006 aTSSs identified here (excluding those that were primarily classified as pTSSs), 881 are within coding sequences, 120 are within 5′ UTRs and 72 are located within 3′ UTRs (note that some aTSS are simultaneously classified in more than one of these three subcategories, Supplementary Figure S6). While we expect that many of the detected antisense transcripts have biological functions, it is difficult to differentiate antisense RNAs with regulatory functions from transcriptional noise. In this regard, Lloréns-Rico et al. (2016) reported that most of the antisense transcripts detected using transcriptomic approaches are a consequence of transcriptional noise, arising at spurious promoters throughout the genome. To investigate the potential significance of the M. smegmatis aTSSs, we assessed the relative impact of each aTSS on local antisense expression levels by comparing the read depth upstream and downstream of each aTSS in our RNAseq expression libraries. We found 318 aTSSs for which expression coverage was ≥10-fold higher in the 100 nt window downstream of the TSS compared to the 100 nt window upstream (Supplementary Table S4). Based on the magnitude of the expression occurring at these aTSS, we postulate that they could represent the 5′ ends of candidate functional antisense transcripts rather than simply products of spurious transcription. However, further work is needed to test this hypothesis. Finally, 78 oTSSs were detected across the M. smegmatis genome. These TSSs may be the 5′ ends of non-coding RNAs or mRNAs encoding previously unannotated ORFs.

Out of the 995 iTSSs identified, 457 were located within the first quarter of an annotated gene (N-iTSSs). In cases where we could not predict a pTSS with high confidence, we considered the possibility that the start codon of the gene was misannotated and the N-iTSS was in fact the primary TSS. Although we do not discount the possibility that functional proteins can be produced when internal transcription initiation occurs far downstream of the annotated start codon, we only considered N-iTSSs candidates for gene start reannotation when there was a start codon (ATG, GTG, or TTG) in-frame with the annotated gene in the first 30% of the annotated sequence. In this way, we suggest re-annotations of the start codons of 213 coding sequences (see Supplementary Figure S5 and Supplementary Table S5). These N-iTSSs were considered to be pTSSs (N-iTSSs → pTSSs) for all further analyses described in this work.

To predict operon structure, we combined 5′ end libraries and RNAseq expression data. We considered two or more genes to be co-transcribed if (1) they had spanning reads that overlapped both the upstream and downstream gene in the expression libraries, (2) at least one TSS was detected in the 5′ end-directed libraries for the first gene of the operon, and (3) the downstream gene(s) lacked pTSSs and iTSSs (for more detail, see section “Materials and Methods”). Thus, we were able to identify and annotate 294 operons with high confidence across the M. smegmatis genome (Supplementary Table S6). These operons are between 2 and 4 genes in length and comprise a total of 638 genes. Our operon prediction methodology has some limitations. For example, operons not expressed in exponential growth phase could not be detected in our study. Furthermore, internal promoters within operons can exist, leading to either monocistronic transcripts or suboperons (Guell et al., 2009; Paletta and Ohman, 2012; Skliarova et al., 2012). We limited our operon predictions to genes that appear to be exclusively co-transcribed, excluding those cases in which an internal gene in an operon can be alternatively transcribed from an assigned pTSS. Finally, our analysis did not capture operons in which the first gene lacked a high-confidence pTSS. Despite these limitations, our approach allowed us to successfully identify new operons as well as previously described operons. Previously reported operons that were captured by our predictions included the furA-katG (MSMEG_6383-MSMEG_6384) operon involved in oxidative stress response (Milano et al., 2001), the vapB-vapC (MSMEG_1283-MSMEG_1284) Toxin–Antitoxin module (Robson et al., 2009) operon, and the ClpP1-ClpP2 (MSMEG_4672-MSMEG_4673) operon involved in protein degradation (Raju et al., 2012).

Most bacterial promoters have two highly conserved regions, the -10 and the -35, that interact with RNA polymerase via sigma factors. However, it was reported that the -10 region is necessary and sufficient for transcription initiation by the housekeeping sigma factor SigA in mycobacteria, and no SigA -35 consensus motifs were identified in previous studies (Cortes et al., 2013; Newton-Foot and Gey van Pittius, 2013; Li et al., 2017; Zhu et al., 2017). To characterize the core promoter motifs in M. smegmatis on a global scale we analyzed the 50 bp upstream of the TSSs. We found that 4,833 of 6,090 promoters analyzed (79%) have an ANNNT motif located between positions -6 to -13 upstream the TSSs (Figure 2A). In addition, 63% of the promoters with ANNNT motifs have a thymidine preceding this sequence (TANNNT). This motif is similar to that previously described in a transcriptome–wide analysis for Mtb (Cortes et al., 2013) and for most bacterial promoters that are recognized by the σ⁷⁰ housekeeping sigma factor (Ramachandran et al., 2014; Sass et al., 2015; Berger et al., 2016; Čuklina et al., 2016; D’arrigo et al., 2016). However, no apparent bias toward specific bases in the NNN region was detected in our study or in Mtb, while in other bacteria such as Escherichia coli, Salmonella enterica, Burkholderia cenocepacia, Pseudomonas putida, and Bacillus subtilis an A/T preference was observed in this region (Jarmer et al., 2001; Ramachandran et al., 2014; Sass et al., 2015; Berger et al., 2016; D’arrigo et al., 2016). We were unable to detect a consensus motif in the -35 region either using MEME server (Bailey et al., 2015) or manually assessing the possible base-enrichment in the -35 region. Analysis of the sequences in the immediate vicinity of TSSs revealed that G and A are the most frequent bases at the +1 position, and C is considerably more abundant at -1 (Figure 2B).

Interestingly, we identified several alternative motifs in the -10 promoter regions of transcripts lacking the ANNNT motif (Figure 2A). One of these, (G/C)NN(G/C)NN(G/C), is likely the signature of M. smegmatis’ codon bias in the regions upstream of iTSSs. The other three sequences are candidate binding sites for alternative sigma factors, which are known to be important in regulation of transcription under diverse environmental conditions. However, the identified consensus sequences differ substantially from those previously described in mycobacteria (Raman et al., 2001, 2004; Sun et al., 2004; Lee et al., 2008a,b; Song et al., 2008; Veyrier et al., 2008; Humpel et al., 2010; Gaudion et al., 2013). The TSSs having these sigma factor motifs and the associated genes are listed in Supplementary Table S7. We next examined the relationship between promoter sequence and promoter strength, as estimated by the read depths in the 5′ end converted libraries. As shown in Figure 2C, the expression levels of transcripts with ANNNT -10 motifs are on average substantially higher than those lacking this sequence. In addition, promoters with the full TANNNT motif are associated with more highly abundant transcripts compared to those having a VANNNT sequence, where V is G, A or C. These results implicate TANNNT as the preferred -10 sequence for the housekeeping sigma factor, SigA, in M. smegmatis. As shown in Figure 2C, expression levels of transcripts having the motif 2 in Figure 2A were significantly increased when compared to the total pool of transcripts lacking the ANNNT motif.

5′ UTRs play important roles in post-transcriptional regulation and translation, as they may contain regulatory sequences that can affect mRNA stability and/or translation efficiency. Whereas in most bacteria 5′ UTR-bearing (“leadered”) transcripts predominate, this is not the case for Mtb, in which near one quarter of the transcripts have been reported to be leaderless (Cortes et al., 2013; Shell et al., 2015b). To investigate this feature in M. smegmatis, we analyzed the 5′ UTR lengths of all genes that had at least one pTSS. We found that for 24% of the transcripts the TSS coincides with the translation start site or produces a leader length ≤5 nt, resulting in leaderless transcripts (Figure 3A). This is less than the 40% reported for M. smegmatis in a smaller TSS-mapping study (Li et al., 2017), and suggests that the proportions of leaderless transcripts are in fact similar for M. smegmatis and Mtb. A total of 1,099 genes (including those re-annotated in section “Mapping, Annotation, and Categorization of Transcription Start Sites”) have leaderless transcripts, and 155 of those (14%) are also transcribed as leadered mRNAs from separate promoters. Two of the pTSSs we validated by 5′ RACE (Supplementary Table S1) belong to leaderless transcripts. For leadered transcripts, the median 5′ UTR length was 69 nt. Interestingly, 15% of the leaders are >200 nt, suggesting that these sequences may contain potential regulatory elements. We then sought to compare the leader lengths of M. smegmatis genes with the leader lengths of their homologs in Mtb. For this analysis we used two independent pTSS mapping Mtb datasets obtained from Cortes et al. (2013) and Shell et al. (2015b) (Figure 3B). To avoid ambiguities, we used only genes that had a single pTSS in both species. Our results show a statistically significant correlation of leader lengths between species, suggesting that similar genes conserve their transcript features and consequently may have related regulatory mechanisms. Additionally, comparison of leaderless transcription in M. smegmatis and Mtb revealed that 62% or 73% of the genes that are only transcribed as leaderless in M. smegmatis also lack a 5′ UTR in MTB, according to Cortes et al. (2013) or Shell et al. (2015b), respectively (Supplementary Table S8). We next assessed if leaderless transcripts are associated with particular gene categories, and found the distribution across categories was uneven (Figure 3C). The three categories “DNA metabolism,” “Amino acid biosynthesis,” and “Biosynthesis of cofactors, prosthetic groups and carriers” were significantly enriched in leaderless transcripts (p-value < 0.05, hypergeometric test), while “Signal transduction,” “Transcription,” and “Transport and binding proteins” appear to have fewer leaderless transcripts.

We next evaluated the presence of the Shine-Dalgarno ribosome-binding site (SD) upstream of leadered coding sequences. For this analysis, we considered those leaders containing at least one of the three tetramers AGGA, GGAG or GAGG (core sequence AGGAGG) in the region -6 to -17 relative to the start codon to possess canonical SD motifs. We found that only 47% of leadered coding sequences had these canonical SD sequences. Thus, considering also the leaderless RNAs, a large number of transcripts lack canonical SD sequences, suggesting that translation initiation can occur through multiple mechanisms in M. smegmatis. We further compared the relative expression levels of leaderless and leadered coding sequences subdivided by SD status. Genes expressed as both leadered and leaderless transcripts were excluded from this analysis. We found that on average, expression levels were significantly higher for those genes with canonical SD sequences than for those with leaders but lacking this motif and for those that were leaderless (Figure 3D). Together, these data suggest that genes that are more efficiently translated have also higher transcript levels. Similar findings were made in Mtb, where proteomic analyses showed increased protein levels for genes with SD sequences compared to those lacking this motif (Cortes et al., 2013).

As GTG or ATG codons are sufficient to initiate leaderless translation in mycobacteria (Shell et al., 2015b; Potgieter et al., 2016), we used this feature to look for unannotated ORFs in the M. smegmatis NC_008596 reference genome. Using 1,579 TSSs that remained after pTSS assignment and gene reannotation using N-iTSSs (see Supplementary Figure S5) we identified a total of 66 leaderless ORFs encoding putative proteins longer than 30 amino acids, 5 of which were previously identified (Shell et al., 2015b). 83% of these ORFs were predicted in other annotations of the M. smegmatis mc²155 or MKD8 genome [NC_018289.1, (Gray et al., 2013)], while 10 of the remaining ORFs showed homology to genes annotated in other mycobacterial species and Helobdella robusta and two ORFs did not show homology to any known protein. The TSS of ORF15 was validated by 5′RACE. These results show that automatic annotation of genomes can be incomplete and highlight the utility of transcriptomic analysis for genome (re)annotation. Detailed information on these novel putative ORFs is provided in Supplementary Table S9.

As our methodology allows us to precisely map RNA cleavage sites in addition to TSSs, we sought to analyze the presence and distribution of cleavage sites in the M. smegmatis transcriptome. mRNA processing plays a crucial role in regulation of gene expression, as it is involved in mRNA maturation, stability and degradation (Arraiano et al., 2010). Mixture modeling identified 3,344 CSs with a posterior probability ≥0.9 (high confidence CSs) (Figure 1A and Supplementary Table S10). To determine the sequence context of the CSs, we used the regions flanking the 5′ ends to generate a sequence logo (Figure 4A). There was a strong preference for a cytosine in the +1 position (present in more than the 90% of the CSs) (Figure 4B), suggesting that it may be structurally important for RNase recognition and/or catalysis.

Cleaved 5′ ends can represent either degradation intermediates or transcripts that undergo functional processing/maturation. In an attempt to investigate CS function, we classified them according to their locations within mRNA transcripts (Figure 4C and Supplementary Table S10). We found that, after normalizing to the proportion of the expressed transcriptome that is comprised by each location category, cleaved 5′ ends are more abundant within 5′ UTRs and intergenic regions of operons than within coding sequences and 3′ UTRs (Figure 4D). Stringent criteria were used in these analyses to avoid undesired bias (Figure 4C and see section “Materials and Methods”). While one would expect the CSs associated with mRNA turnover to be evenly distributed throughout the transcript, enrichment of CSs within the 5′ UTRs as well as between two co-transcribed genes may be indicative of cleavages associated with processing and maturation. Alternatively, these regions may be more susceptible to RNases due to lack of associated ribosomes. Here we predicted with high confidence that at least 101 genes have one or more CSs in their 5′ UTRs (Supplementary Table S11).

We detected cleaved 5′ ends within the coding sequences of 18% of M. smegmatis genes, ranging from 1 to over 40 sites per gene. We analyzed the distribution of CSs within coding sequences (Supplementary Figure S7), taking into consideration the genomic context of the genes. When analyzing the distribution of CSs within the coding sequences of genes whose downstream gene has the same orientation, we observed an increase in CS frequency in the region near the stop codon (Supplementary Figure S7A). However, when only coding sequences having a downstream gene on the opposite strand (convergent) were considered, the distribution of CSs through the coding sequences was significantly different (p-value < 0.0001, Kolmogorov-Smirnov D test) with the CSs more evenly distributed throughout the coding sequence (Supplementary Figure S7B). This suggests that the cleavage bias toward the end of the genes observed in Supplementary Figure S7A may be due to the fact that many of these CSs are actually occurring in the 5′ UTRs of the downstream genes. In cases where the TSS of a given gene occurs within the coding sequence of the preceding gene, a CS may map to both the coding sequence of the upstream gene and the 5′ UTR of the downstream gene. In these cases, we cannot determine in which of the two transcripts the cleavage occurred. However, cleavages may also occur in polycistronic transcripts. We therefore assessed the distributions of CSs in the operons predicted above. The distribution of CSs in genes co-transcribed with a downstream gene showed a slight increase toward the last part of the gene (Supplementary Figure S7C). This may reflect cases in which polycistronic transcripts are cleaved near the 3′ end of an upstream gene, as has been reported for the furA-katG operon, in which a cleavage near the stop codon of furA was described (Milano et al., 2001; Sala et al., 2008; Taverniti et al., 2011). The furA-katG cleavage was identified in our dataset, located 1 nt downstream of the previously reported position. A similar enrichment of CSs toward stop codons was also observed in a recent genome-wide RNA cleavage analysis in S. enterica (Chao et al., 2017), although in this case the high frequency of cleavage may be also attributed to the U preference of RNase E in this organism, which is highly abundant in these regions.

The sequence contexts of TSSs (Figure 2B) and CSs (Figure 4A) were markedly different, as G and A were highly preferred in the TSS +1 position whereas C was highly preferred in the CS +1 position, and TSSs were associated with a strong overrepresentation of ANNNT -10 sites while CSs were not. These sequence-context differences not only provide validation of our methodology for distinguishing TSSs from CSs, as discussed above, but also provide a means for making improved predictions of the nature of 5′ ends that could not be categorized with high confidence based on their converted/non-converted library coverage alone. Taking advantage of these differences, we sought to obtain a list of additional putative TSSs and CSs. Thus, of the 5′ ends that were not classified with high confidence by mixture modeling, we selected those that had an appropriately positioned ANNNT motif upstream and a G or an A in the +1 position and classified them as TSSs with medium confidence (Supplementary Table S12). In the same way, 5′ ends with a C in the +1 position and lacking the ANNNT motif in the region upstream were designated as medium confidence CSs (Supplementary Table S13). In this way, we were able to obtain 576 and 4,838 medium confidence TSSs and CSs, respectively. Additional validation of a medium confidence TSS was performed for gene MSMEG_0063 using 5′RACE. We were able to corroborate that, as predicted, transcription of this gene is initiated 139 bp upstream the coding sequence and that either deletion or mutation of the predicted -10 promoter region dramatically decreased transcription initiation (Supplementary Figure S8). These results support the value of TSS prediction based on -10 promoter region motif and base composition at +1 position, and highlight the importance of the -10 ANNNT promoter motif for mycobacterial transcription. Three medium confidence CSs (86927+, 87293+, and 5038902-) were also validated using 5′ RACE. Although we are aware of the limitations of these predictions, these lists of medium confidence 5′ ends provide a resource that may be useful for guiding further studies. 5′ ends that did not meet the criteria for high or medium confidence TSSs or CSs are reported in Supplementary Table S14.

We sought to study the global changes occurring at the transcriptomic level in oxygen limitation employing a system similar to the Wayne model (Wayne and Hayes, 1996) (see section “Materials and Methods”). Two timepoints were experimentally determined in order to evaluate transcriptomic changes during the transition into hypoxia (Supplementary Figure S1). A different enzyme was used for conversion of 5′ triphosphates to 5′ monophosphates in these 5′-end libraries, and it appeared to be less effective than the enzyme used for the 5′ end libraries in Dataset 1. As a consequence, our ability to distinguish TSSs from CSs de novo in these datasets was limited. However, we were able to assess changes in abundance of the 5′ ends classified as high-confidence TSSs or CSs in Dataset 1, as well as identify a limited number of additional TSSs and CSs with high confidence (Supplementary Figure S4 and Supplementary Table S3). Corresponding RNAseq expression libraries revealed that, as expected, a large number of genes were up and downregulated in response to oxygen limitation (Supplementary Figure S9 and Supplementary Table S15). We next investigated the transcriptional changes in hypoxia by assessing the relative abundance of TSSs in these conditions. We found 318 high-confidence TSSs whose abundance varied substantially between exponential phase and hypoxia (Supplementary Table S16). A robust correlation was observed between the pTSS peak height in the 5′-end-directed libraries and RNA levels in the expression libraries for hypoxia (Supplementary Figure S10). In an attempt to identify promoter motifs induced in hypoxia, we analyzed the upstream regions of those TSSs whose abundance increased (fold change ≥2, adjusted p-value ≤ 0.05). Interestingly, we detected a conserved GGGTA motif in the -10 region of 56 promoters induced in hypoxia using MEME (Figure 5A and Supplementary Table S16). This motif was reported as the binding site for alternative sigma factor SigF (Rodrigue et al., 2007; Hartkoorn et al., 2010; Humpel et al., 2010). Additionally, the extended -35 and -10 SigF motif was found in 44 of the 56 promoter sequences (Figure 5A and Supplementary Table S16). SigF was shown to be induced in hypoxia at the transcript level in Mtb (Iona et al., 2016) and highly induced at the protein level under anaerobic conditions using the Wayne model in Mycobacterium bovis BCG strain and Mtb (Michele et al., 1999; Galagan et al., 2013). In M. smegmatis, SigF was shown to play a role under oxidative stress, heat shock, low pH and stationary phase (Gebhard et al., 2008; Humpel et al., 2010; Singh et al., 2015) and sigF RNA levels were detected in exponential phase at a nearly comparable level to sigA (Singh and Singh, 2008). Here, we did not detect significant changes in expression of the sigF gene in hypoxia at the transcript level. However, this is consistent with reported data showing that sigF transcript levels remain unchanged under stress conditions in M. smegmatis (Gebhard et al., 2008), as it was postulated that SigF is post-transcriptionally modulated via an anti-sigma factor rather than through sigF transcription activation (Beaucher et al., 2002). We noted that, in the case of TSSs whose abundance was reduced in hypoxia, almost the totality of the promoters contains the -10 ANNNT σ⁷⁰ binding motif. We then examined the presence of SigF motif in the regions upstream of 5′ ends that were not classified as high confidence TSSs. We speculate that 5′ ends associated with this motif may be potential TSSs triggered by hypoxia. We found 96 additional putative TSSs that were (1) overrepresented in hypoxia and (2) associated with appropriately-spaced SigF motifs (Supplementary Table S17). Three of the hypoxia-induced genes with SigF motifs (MSMEG_3460, MSMEG_4195 and MSMEG_5329) have homologous genes induced in hypoxia in Mtb (Park et al., 2003; Rustad et al., 2008).

It is well known that under anaerobic conditions mycobacteria induce the DosR regulon, a set of genes implicated in stress tolerance (Rosenkrands et al., 2002; O’Toole et al., 2003; Park et al., 2003; Roberts et al., 2004; Rustad et al., 2008; Honaker et al., 2009; Leistikow et al., 2010). The DosR transcriptional regulator was highly upregulated at both hypoxic timepoints in the expression libraries (13 and 18-fold at 15 and 24 h, respectively, Supplementary Figure S9) and 30 out of the 49 DosR-activated genes (Berney et al., 2014) were upregulated in our dataset. Thus, we hypothesized that the DosR binding motif should be present in a number of regions upstream the TSSs that were upregulated in hypoxia. Analysis of the 200 bp upstream the TSSs using the CentriMo tool for local motif enrichment analysis (Bailey and Machanick, 2012) allowed us to detect putative DosR motifs in 13 or 53 promoters, depending on whether a stringent (GGGACTTNNGNCCCT) or a weak (RRGNCYWNNGNMM) consensus sequence was used as input (Lun et al., 2009; Berney et al., 2014; Gomes et al., 2014) (Supplementary Table S16). At least two of the 13 genes downstream of these TSSs were previously reported to have DosR motifs by Berney et al. (2014) and RegPrecise Database (Novichkov et al., 2013) and two others are homologs of genes in the Mtb DosR regulon that were not previously described in M. smegmatis as regulated by DosR (Supplementary Table S16).

We then used CentriMo to search for DosR motifs in the regions upstream of 5′ ends that were not classified as high confidence TSSs, given that TSSs derived from hypoxia-specific promoters may have been absent from Dataset 1. We found 36 putative TSSs associated with 20 different genes (Supplementary Table S18), of which 11 have been shown to have DosR binding motifs (Berney et al., 2014). Five of these are homologs of genes in the Mtb DosR regulon.

There is evidence that mycobacterial mRNA is broadly stabilized under hypoxia and other stress conditions (Rustad et al., 2013; Ignatov et al., 2015). Thus, we anticipated that RNA cleavage should be reduced under hypoxia as a strategy to stabilize transcripts. We compared the relative abundance of each high confidence CS in stress and in exponential phase (Figure 5B) and found that RNA cleavage is significantly reduced in both hypoxia 15 and 24h on a global scale (Figure 5C). In contrast, relative abundance of TSSs did not decrease in these conditions, indicating that the reduction in CSs is not an artifact of improper normalization (Figure 5B). When the ratios of CSs abundance in hypoxia/normal growth of individual genes were analyzed, we observed the same behavior (Supplementary Figure S11). These results indicate that the number of cleavage events per gene decreases during adaptation to hypoxia, which could contribute to the reported increases in half-life (Rustad et al., 2013).