Strains used in this study are listed in Table 1. Escherichia coli strain DH5a was grown in lysogeny broth (LB) medium. mc²155 wild type strain and its derivatives were grown at 37°C in Middlebrook 7H9 medium (Difco Becton Dickinson, USA) supplemented with 0.5% (v/v) glycerol and 0.05% (v/v) Tween 80, or on Middlebrook 7H10 agar (Difco Becton Dickinson, USA) supplemented with 0.5% (v/v) glycerol (Tang et al., 2014). When necessary, antibiotics were added to the culture containing either 50 μg/mL of kanamycin or 100 μg/mL of ampicillin.

mc²155 cells equivalent to 30 OD₆₀₀ (e.g., 30 mL of 1 OD₆₀₀ of one culture) were harvested from each sample (three growth phases each in two biological replicates) and then quickly-frozen in liquid nitrogen. The bacteria pellet was transferred into a mortar containing liquid nitrogen, and was lysed by constant grinding (in the presence of liquid nitrogen), into a fine powder. The powder was transferred to a 2 mL tube containing 1 mL TRIzol Reagent (Invitrogen, CA, USA), and shaken vigorously, until the powder dissolved uniformly. Then, 200 μL of chloroform was added to the mixture and mixed vigorously for 15 s. The mixture was incubated on ice for 5 min, followed by centrifugation at 10,000 × g for 15 min at 4°C. Equal volumes (400 μL) of upper clear supernatant and isopropanol were mixed in a fresh 1.5 mL tube to precipitate RNA. It was mixed immediately by gently inverting 8~10 times, and stored on ice for 15 min, followed by centrifugation at 10,000 × g for 30 min at 4°C. The upper clear phase was discarded, and 1 mL of 75% ethanol was added to wash the RNA pellet. The sample was gently washed by inverting 4–6 times, followed by centrifugation at 5,000 × g for 2 min at 4°C. Ethanol washing was repeated twice. Following the wash steps, ethanol was discarded. RNA was eventually dried by exposing to air for 5 min (air-dry) at room temperature. Finally, the RNA pellet was dissolved in 50 μL of RNase-free water, and stored at −80°C. Note that all reagents and materials used in this experiment were RNase-free.

For RNA-seq, 10 μg of total RNA from each sample was first treated with RNase-free DNase I (Takara, Japan) to prevent the potential contamination of genomic DNA. Ribosomal RNAs were removed using the RiboZero rRNA removal kit (Epicentre, USA) for gram-positive organisms prior to sequencing analysis. Hundered nanograms of rRNA-depleted RNA from each sample was fragmented into 200–300 nts and used as a template for randomly primed PCR. Strand-specific cDNA libraries were prepared by standard techniques for subsequent Illumina sequencing using the mRNA-seq Sample Prep kit (Illumina, USA). The resulting cDNAs were sequenced on an Illumina HiSeq 2500 (Shao et al., 2015).

The quality of the raw sequence data was first assessed by FastQC, the sequences further filtered and trimmed by Trimmomatic (Bolger et al., 2014). Clean reads were then used to map to the reference genome of mc²155 (NC_008596.1) by Bowtie2 (Langmead and Salzberg, 2012). BEDTools (Quinlan and Hall, 2010) analysis was performed on the alignments generated by Bowtie2 in order to quantify the mapped reads. To facilitate the comparison of expression levels of different samples, reads for CDSs were normalized with the total read count per sample and presented in the form of Reads Per Kilobase Per Million Reads (RPKM). The quantity of differential expressions of all transcripts was obtained using DEGseq (Wang et al., 2010) to include fold change values. In context, DEGs are those with two-fold changes or more, with FDR (false discovery rate) <0.001. As for single-nucleotide resolution transcriptome map, cleaned reads were mapped to mc²155 genome by BLASTN to get their genome location. Then Perl script was used to count the expression level of each nucleotide. Data were visualized using Artemis software.

To verify the TSSs identified by RNA-seq, 5′-RACE was performed using a tobacco acid pyrophosphatase (TAP)-based 5′-RACE kit according to manufacturer's instructions as previously described (Ali et al., 2017). The selected regions were amplified from 5′-RACE cDNA library using adapter primers (outer and inner primers) and gene specific primers (Table S1). After multiplex amplification, PCR products were purified and ligated into pMD19-T vector, and then transformed into E. coli DH5α (Table 1). Clones were sequenced and the first nucleotides next to 5′-RACE adaptor were considered to be the TSSs.

N-terminal sequencing was performed on proteins translated from leaderless mRNAs to identify their translation initiation amino acid residues. To get the test proteins, leaderless genes fused with 6 × His tag were over-expressed under the control of their own promoters based on the expression vector pMV261, and then, respectively, transformed into mc²155 strain to obtain a series of derivative strains. The native promoter could hardly be induced, so two highly expressed leaderless genes MSMEG_4921 and MSMEG_6422 were selected. For protein purification, the two mc²155 derivatives were harvested at mid-exponential phase, followed by cell lysis using French press. The clear lysate was loaded onto a Ni²⁺-NTA affinity column for protein purification. Purified proteins were applied on a 15% SDS-PAGE and then transferred to PVDF blotting membrane for Edman sequencing using PPSQ-31A protein sequencer (Shimadzu, Japan). The strains and primers used are listed in Table 1 and Table S1, respectively.

Operons prediction was based on the transcriptome map, visualized in Artemis. The prediction was carried out according to previously described method (Wang et al., 2013b). Briefly, genes in an operon share the same orientation, and transcripts should cover the intergenic region of the two adjacent genes. In addition, the ratio of the two gene expression levels should be <2. According to these two principles, we manually inspected the transcription of all 6,947 genes across the whole genome.

β-galactosidase activity was performed in mc²155 by translational/transcriptional fusion of promoter to lacZ gene based on the expression vector pMV261. For translational fusion, the promoter regions, 5′-UTRs and partial N-terminal amino acids were fused with lacZ gene, while promoter and 5′-UTRs only were included in transcriptional fusion. Phsp60 promoter was used as control when assessing the relative β-galactosidase activity of selected promoters (Yang et al., 2012). Reporter plasmids were transformed into wild type strain to obtain corresponding recombinant strains. All strains were grown in 7H9 medium at 37°C. β-galactosidase activity was determined as previously described (Miller, 1972). The strains and primers used are listed in Table 1 and Table S1, respectively.

The growth curve of mc²155 in 7H9 broth is a prerequisite for its transcriptome profile which comprises a series of growth-phase-dependent RNA-seq experiments. Results indicate that mc²155 reached mid-exponential phase at 16 h; early-stationary phase at 26 h and mid-stationary phase at 39 h (Figure S1). Two biological replicates were harvested from each of these three time points, and their total RNAs were extracted for subsequent experiments. All of the six RNA samples obtained from three different growth phases (two biological replicates) were found to exhibit high integrity with RNA Integrity Numbers (RINs) more than 9.5. Strand-specific cDNA libraries were then constructed and sequenced on an Illumina HiSeq 2000 platform. After removing low quality reads, 10 million clean reads were obtained for each library. The complete genome and trancriptome information of mc²155 are shown in Figure 1.

In order to temporally investigate distinct biological processes, we analyzed gene expression profiles of the growth-phase-dependent RNA-seq data (Table S2). Comparison of RPKM values for each growth phase revealed significantly lower transcriptional levels for mid-stationary phase genes than those in the other two growth phases (Wilcoxon test, p < 0.001; Figure 2A). The differences in gene expression levels among the three growth phases were also analyzed with respect to different Cluster of Orthologous Groups of proteins (COG). DEGs were analyzed using mid-exponential phase RNA-seq data as control. In the mid-stationary phase, DEGs (especially down-regulated genes) are obviously increased in some clusters when compared to early-stationary phase (Figures 2B,C). We propose that in a cluster when the number of down regulated genes is far more than up-regulated genes, the cluster would be considered as down-regulated and vice versa. Thus, in the early-stationary phase (26 h), “carbohydrate transport and metabolism [G]” cluster was up-regulated, whereas “coenzyme transport and metabolism [H]”, “lipid transport and metabolism [I]” and “translation, ribosomal structure, and biogenesis [J]” clusters were down-regulated (Figure 2B). In the mid-stationary phase, “nucleotide transport and metabolism [F]”, “coenzyme transport and metabolism [H]”, “translation, ribosomal structure and biogenesis [J]”, and “replication, recombination and repair [L]” clusters were down-regulated; however, no cluster was up-regulated (Figure 2C).

Furthermore, DEGs were subjected to Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways analysis to find out temporally associated biological functions. As shown in Figure 3, pathways with member genes significantly regulated (Q-value < 0.05) were marked in bold and underlined. As mentioned above, transcriptional levels of most genes were substantially reduced in mid-stationary phase; however, most of the genes in “inositol phosphate metabolism”, “degradation of aromatic compounds” and “oxidative phosphorylation” pathways were up-regulated (Figure 3A and Figure S2). In contrast, the majority of genes were significantly down-regulated, especially genes in “translation” (including “ribosome” and “aminoacyl-tRNA biosynthesis” pathways) process (Figure 3B and Figure S2).

Transcriptional start sites (TSSs) typically refer to the 5′-ends of primary transcripts, but not 5′-ends of the processed transcripts. In our transcriptome data, most of the 5′-ends of transcripts (except for those of processed transcripts) are indeed their TSSs, thus we used the term “TSSs” in the following text. To locate TSSs in genome-wide scale, clean reads were mapped to mc²155 genome to generate single-nucleotide resolution transcriptome maps. However, the uneven transcriptome map caused by sequencing makes it hard to identify all the TSSs, especially internal TSSs which could typically imply the existence of small RNAs. Therefore, to improve the accuracy and reliability of our study, we only focused on the TSSs found at the beginning of mRNA transcripts. We manually inspected the transcriptome maps separately and successively. As shown in Figure 4A, the whole view of strand-specific transcriptional maps of MSMEG_3068-3078 gene cluster was uneven, which is common in RNA-seq caused by sequencing bias. The first nucleotide significantly enriched at the beginning of transcript is considered to be the TSS. Single-nucleotide resolution of transcriptional maps showed that the transcriptions of MSMEG_3071 and MSMEG_3070 genes were clearly enriched at TSSs (Figures 4B,C). In this way, we totally identified 2139 TSSs throughout the mc²155 genome (Table S3). To further confirm TSSs identified by RNA-seq, eight transcripts were randomly selected and verified by the canonical 5′-RACE method (Table S4). The results showed that TSSs identified by RNA-seq and 5′-RACE experiments overlapped at an identical nucleotides, which validates the accuracy of RNA-seq based TSSs identification. In general, promoters are located immediately upstream of TSSs. Thus, sequences 50 nts upstream of identified TSSs were scanned for potential conserved promoter motif by MEME. As a result, a conserved −10 motif (TANNNT) was found ranging 7–12 bp upstream of the identified TSSs in 69.5% of the genes, but no conserved −35 motif could be identified. These findings are similar to a previous study in M. tuberculosis (Cortes et al., 2013).

MSMEG_2196 is a c-di-GMP synthetase encoding gene, the TSS of which was first determined in a previous study (Bharati et al., 2013). In our study, its TSS was detected by RNA-seq and was validated by 5′-RACE (Figure 5). Interestingly, the reported TSS is 88 nts downstream of our identified location. Moreover, a conserved −10 motif could be found upstream of our identified TSS but not in the previous one, indicating that our study could provide more credible data for gene organization. Taken together, our systematic identification of TSSs in genome scale could facilitates the research of gene expression regulation in mc²155.

Based on TSS identification, we systematically explored the 5′-UTR length throughout the transcriptome (Figure 6A and Table S3). As shown in Figure 6B, 40.5% of the 5′-UTRs were <10 nts in length, and were considered to be leaderless mRNA; 28.6% of transcripts have 5′-UTR length from 11 to 50 nts, and 14.8% from 51 to 100 nts. In addition, 9.1% of transcripts exhibited long 5′-UTRs (with length more than 100 nts). Previous study indicated that long 5′-UTRs might hide regulatory elements such as riboswitch and transcription-attenuation region, which would regulate the transcription of downstream genes in a more sophisticated way (Breaker, 2011; Wang et al., 2013a; Rosinski-Chupin et al., 2014; Dersch et al., 2017). It is worth mentioning that about 7% TSSs were identified downstream of the annotated start codon, which were rather unlikely and were therefore considered to be mis-annotated.

We have observed that ~40.5% (866/2139) of the transcripts with annotated TSSs had a 5′-UTR <10 nts in length, likely generating leaderless mRNAs. Moreover, the TSSs of 35.48% (759/2139) mRNAs overlapped with the first nucleotide of annotated start codon (AUG/GUG). In bacteria, the Shine-Dalgarno (SD) sequence helps to recruit the ribosome to the mRNA to initiate translation. One would ponder how ribosome binds to mRNA to start the translation procedure without an SD sequence. Can these mRNA be efficiently translated? Do these widespread leaderless mRNAs use other specific translation mechanisms or are they simply mis-annotated? Therefore, it is important to confirm the translation of leaderless mRNA and identify the initial amino acid of the translated protein. Shell et al. (2015) tried to reveal the widespread translation of leaderless mRNA in mc²155 by ribosome profiling; however, the result could only indicate that ribosome did bind to the leaderless mRNAs but could not directly verify the translation of leaderless mRNAs. Moreover, they had also performed a proteomics study (Orbitrap LC-MS) on M. tuberculosis H37Rv to confirm this feature. As a result, a large number of proteins translated from leaderless mRNAs were identified. However, the key point in the validation of leaderless feature is to determine whether the first N-terminal amino acid of the translated protein is in accordance with the initial nucleotides of leaderless mRNA. LC-MS cannot be used to efficiently identify the initial amino acid though some peptides detected in that study were identical to the annotated start codons. Thus, we performed Edman degradation, which is the most credible method in initial amino acid identification, on the MSMEG_4921 and MSMEG_6422 proteins both translated from leaderless mRNAs (Figures 7A,B). To get the test protein, leaderless genes fused with 6 × His tag were over-expressed under the control of their own promoters based on the expression vector pMV261, and were respectively transformed into mc²155 strain to obtain a series of derivative strains. Purified proteins were applied on a 15% SDS-PAGE and then transferred to PVDF blotting membrane for Edman sequencing. The results indicated that their five N-terminal amino acids were identical to those annotated in the genome (Figures 7C,D). This strongly confirmed the leaderless feature of these leaderless mRNAs. Noteworthy, the N-terminal translation initiator Met was removed by methionine amino peptidase (MetAP), which is often crucial for the function and stability of proteins. Collectively, leaderless mRNA seems to be the common feature of mycobacteria. Further studies are required to understand how these leaderless mRNAs could be successfully translated into proteins. Several studies have shown that the unique translation can be initiated with high efficiency at leaderless transcripts via the undissociated 70S/80S ribosome and the initiator tRNA (Moll et al., 2004; Udagawa et al., 2004; Andreev et al., 2006; Giliberti et al., 2012). However, it is unclear whether mycobacterial strains employ this mechanism to initiate the translation of the high-proportion leaderless mRNA.

The majority of sequenced bacterial genomes were annotated by algorithm-based software, which can be error-prone. As shown in Figure 6, the lengths of 5′-UTRs of 150 transcripts were negative, indicating that these corresponding genes were very likely mis-annotated in these assay conditions (three different growth phases in 7H9 medium). Further explorations found that 124 out of these 150 transcripts used AUG/GUG as the first three nucleotides at TSS (Table S5), which were similar to the leaderless mRNAs. Given the fact that leaderless mRNAs are widespread in mc²155, and an AUG/GUG at TSS is sufficient for the initial translation of leaderless mRNA in mc²155 (Shell et al., 2015), we therefore propose that in these assay conditions, the 124 mis-annotated genes were transcribed into leaderless mRNA. Furthermore, distances between the identified TSSs and the first nucleotide of annotated start codons were in multiples of three in length, indicating that the corrected gene structural annotation does not shift the original open reading frame. Figure 8 schematically shows the transcriptional maps of two mis-annotated proteins MSMEG_1874 and MSMEG_6901, as well as the corrected versions. Together, we propose that RNA-seq could be highly useful to correct gene mis-annotation.

Some genes in prokaryotes are organized into operons, which are then transcribed into polycistronic transcripts. It is typically considered that genes in an operon are functionally related. Zhang and He (2013) reported that overexpression of MSMEG_6080 in mc²155 led to cell expansion and aggregation; however, co-overexpression of MSMEG_6080 and MSMEG_6078, located in a similar operon, alleviated the phenotypes described above instead. Furthermore, a gene in the middle of an operon, which can be independently transcribed in response to certain stimuli, is called sub-operon. Therefore, identification of genes that are grouped together into operons/sub-operons is important to studies of gene function and complex regulatory networks.

To investigate the structures of operons and sub-operons in genome scale, we manually inspected the transcription of all 6,947 genes across the genome, and finally identified 2,233 operons/sub-operons (Table S6). The 2,233 operons/sub-operons were further classified into five groups, according to the prediction of operons by DOOR (Mao et al., 2014). These five groups are: (1) confirmed, which means that our results were consistent with those predicted by DOOR (Figure 9A); (2) extended, in which we found that more genes could be transcribed into the polycistronic mRNA (Figure 9B); (3) dismissed, where less genes were found to be transcribed together (Figure 9C); (4) new, where some genes were transcribed together but were not annotated by DOOR (Figure 9D); and (5) alternative, indicating the existence of sub-operon (Figure 9E). To further illustrate the reliability of our results, several operons of each group were randomly selected for reverse transcription PCR (RT-PCR), and the results were consistent with RNA-seq data. The selected operons/sub-operons in each group is shown in Table 2.

As mentioned above, mc²155 is widely accepted as a model organism in mycobacterial researches and has usually been used to overexpress homologous proteins from M. tuberculosis and M. bovis due to its non-pathogenic and fast-growing features (Soares et al., 2014; Deng et al., 2016). The temperature-sensitive Phsp60 promoter could induce the expression of downstream genes on heat shock conditions, typically via a temperature shift from 37 to 45°C (Batoni et al., 1998). The inducibility makes Phsp60 the most widely used promoter in mycobacteria. However, there are also some problems existing in the inducible expression system. Therefore, more promoters need to be discovered for driving gene expression in mycobacteria.

It is widely known that many non-coding RNAs such as rRNA, tRNA, 6S RNA, possess a highly active promoter. Here, we intended to screen alternate highly active promoters for gene expression. According to RNA-seq data, 22 gene candidates with high RPKM values were selected (Table S7). By fusing these promoter sequences to lacZ gene, we successfully certified that 8 out of 22 candidates exhibiting high β-galactosidase activities that are 3–5 times higher than the control promoter Phsp60 (Figure 10). Additionally, some promoters showed temporal characteristics. As shown in Figures 10C,D, the promoter activities of P_{MSMEG_3050} and P_{MSMEG_3084} increased with time. In contrast, P_{MSMEG_4891} only exhibited high promoter activities in the early growth phase (Figure 10E). However, P_{MSMEG_5081} showed a relatively constant expression level during the whole detection period (Figure 10F). It is noteworthy that the β-galactosidase activity assays do not correlate well with RNA-seq data. This may be due to the fact that RNA-seq data only denote the transcriptional level. However, β-galactosidase activity assays reflect the expression levels of both transcription and translation. The disagreement between mRNA and protein levels usually happens, since many modifications may occur either at post-transcriptional or translational levels (McCarthy and Gualerzi, 1990; Alifano et al., 1994; Kozak, 2005).

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.