Comparative transcriptomics across the prokaryotic tree of life.

Whole-transcriptome sequencing studies from recent years revealed an unexpected complexity in transcriptomes of bacteria and archaea, including abundant non-coding RNAs, cis-antisense transcription and regulatory untranslated regions (UTRs). Understanding the functional relevance of the plethora of non-coding RNAs in a given organism is challenging, especially since some of these RNAs were attributed to 'transcriptional noise'. To allow the search for conserved transcriptomic elements we produced comparative transcriptome maps for multiple species across the microbial tree of life. These transcriptome maps are detailed in annotations, comparable by gene families, and BLAST-searchable by user provided sequences. Our transcriptome collection includes 18 model organisms spanning 10 phyla/subphyla of bacteria and archaea that were sequenced using standardized RNA-seq methods. The utility of the comparative approach, as implemented in our web server, is demonstrated by highlighting genes with exceptionally long 5'UTRs across species, which correspond to many known riboswitches and further suggest novel putative regulatory elements. Our study provides a standardized reference transcriptome to major clinically and environmentally important microbial phyla. The viewer is available at http://exploration.weizmann.ac.il/TCOL, setting a framework for comparative studies of the microbial non-coding genome.

INTRODUCTION

The growing availability of next-generation sequencing technologies has led to a burst of studies characterizing the transcriptomes of many prokaryotic species, often charting transcription start sites (TSSs) at single-nucleotide resolution (1–9). These studies unraveled unexpected complexity in the transcriptomes of bacteria and archaea including numerous long regulatory 5′UTRs, non-coding RNAs, alternative operon structures, internal promoters and abundant cis-antisense transcription (10–12).

Comparative transcriptome studies of various prokaryotes had shown that transcriptional maps may differ even among closely related species, challenging functional interpretations of ncRNA transcription with suspected transcriptional noise (13–17). While the functional importance of non-conserved transcription remains debated (18), cases in which evolutionary conservation is observed were shown to have the potential to highlight functionally important features in the organisms’ transcriptional dynamics (19–21). To date, comparative transcriptomic studies were limited to few organisms, usually comparing closely related species.

A major challenge in comparative transcriptomics studies is the construction of reliable and comparable transcriptome maps across multiple diverse species. The first published transcriptome maps were based on manual curation of multiple putative TSSs. With the dramatic reduction of sequencing costs, however, using manual curation for construction of transcriptome maps becomes increasingly prohibitive because it is not scalable to large number of organisms and multiple conditions. Thus, in the last couple of years several methods had been developed to automatically infer TSSs. Some of these methods are tailored analyses for differential 5′-RNA sequencing data (dRNA-seq, (11)) as in the primary transcriptome of H. pylori (1). These include methods that are primarily based on analyzing two RNA-seq libraries, one untreated and one enzymatically treated to enrich for primary transcripts (22–24). Additional methods that analyze TSSs integrate comparative signals from closely related species (25), while other methods attempt to delineate transcript boundaries (both TSS and termination) based on statistically significant local differences in RNA-seq coverage (26,27). However, methods that can infer primary TSSs based on multiple transcriptomic and genomic features, taken together, are still lacking.

In this study, we present the TCOL web-server (Transcriptomes Compared across the tree Of Life) that includes several important novel components. First, TCOL encompass the largest compendium of comparable transcriptomes to date, spanning representatives of the major phyla of the prokaryotic tree of life. Second, all transcriptomes were sequenced using standardized approach, combining RNA-seq and 5'-end sequencing, generating comparable data. Third, we reconstructed the transcriptome maps using a computational pipeline that utilized machine learning to accurately infer primary TSSs based on multiple genomic and transcriptomic features. Fourth, TCOL includes a comparative transcriptome browser, allowing for simultaneous viewing of transcriptome maps for homologous genes across divergent organisms. Finally, the TCOL transcriptomes are BLAST-searchable with query sequences provided by users, broadening the utility of this web-server to the microbiology community.

MATERIALS AND METHODS

Culture conditions

Most data were derived from previously published studies, with the exception of the following organisms for which data were obtained within this study: E. coli was grown in LB medium overnight, diluted 1:200 into fresh LB and incubated at 37°C until OD600nm of ∼0.3 (early-mid log). B. subtilis was grown in TB, similarly diluted into fresh TB until reaching mid log OD600nm of ∼0.6. G. oxydans 621H was grown in pH 6 in constant temperature of 30°C in two media: (i) mannitol fermentation (50 mM glucose) (ii) glycerol fermentation (50 mM glycerol) until reaching OD600nm of 1.12 and OD600nm of 0.9, for mannitol and glycerol, respectively. C. acetobutylicum was cultured in phosphate limited conditions at (i) pH 4.5 (ii) pH 5.7 until OD600nm of 3.7. T. thermophilus was cultured in TB medium at 70°C until OD600nm of 3.0. S. acidocaldarius was cultured in ‘Brock medium’ (28) with 0.1% tryptone added, until early stationary phase. For all other organisms culture conditions details were published (Supplementary Table S1).

RNA-seq protocols and reads mapping

Total RNA-seq (strand insensitive) and 5′ end transcriptome (strand sensitive) libraries were produced as described previously (19,29) (Supplementary Materials). List of transcriptomes in this study is depicted in Table 1, and sequencing depth and growth conditions of all organisms are detailed in Supplementary Table S1. Briefly, the 5′ differential RNA-seq protocol included the construction of 2 libraries: A treated library in which RNA is incubated with Tobacco Acid Pyrophosphatase (TAP, Epicentre), which we term TAP(+) and an untreated library termed TAP(−). The relative number of TAP(+) and TAP(−) reads that are mapped to the same site is indicative of the likelihood of that site to represent primary TSS. Higher TAP(+)/TAP(−) ratios typically represent primary TSSs while lower values represent processing sites (Supplementary Materials).

Table 1.

The transcriptome maps in the study. A quantitative summary of major transcriptional features including various types of transcriptional start sites (TSSs), inferred small RNAs (sRNAs) and inferred operons

Organism	Accession	Growth conditions	genes	gTSS	iTSS	aTSS	nTSS	sRNAs	Operonsa	RNA-seq and TSS data
Bacillus subtilis	NC_000964	Terrific Broth, mid-log phase	4420	1292	20	29	128	75	1242	This study
Desulfovibrio vulgaris	NC_002937	DSMZ medium	3466	1365	204	407	215	99	1193	(42)
Escherichia coli	NC_000913	mid-log phase	4495	1118	10	35	83	36	1059	This study
Listeria monocytogenes	NC_003210	log phase 37°C; log phase 30°C; hypoxia; stationary phase 37°C; ΔsigB; ΔprfA	2936	1388	26	19	109	36	1282	(19)
Bdellovibrio bacteriovorus	NC_005363	attack phase; growth phase	3627	407	8	16	46	23	395	(43)
Catenulispora acidiphila	NC_013131	DSMZ medium	8984	1996	0	96	197	118	1908	(42)
Clostridium acetobutylicum	NC_003030	pH 4.5; pH 5.7	3777	866	5	1	112	26	839	This study
Gluconobacter oxydans	NC_006677	mannitol; glycerol	2499	1268	135	642	168	87	1113	This study
Kangiella koreensis	NC_013166	DSMZ medium	2681	627	16	16	49	16	601	(42)
Lactobacillus brevis	NC_008497	ATCC medium	2265	777	6	7	71	29	764	(42)
Lactococcus lactis	NC_002662	M17 medium	2400	746	5	7	82	31	714	(42)
Pseudomonas aeruginosa	NC_008463	28°C; 37°C (LB medium)	5975	2224	5	195	321	158	1728	(29)
Synechococcus WH7803	NC_009481	Artificial seawater medium	2585	1176	0	45	31	9	1125	(21)
Synechococcus WH8102	NC_005070	Artificial seawater medium	2581	1005	0	52	44	5	957	(21)
Spirochaeta aurantia	Saur_Contig1177	DSMZ medium	4124	1284	8	212	84	34	1229	(42)
Sulfolobus acidocaldarius	NC_007181	yeast extract, stationary phase	2330	1050	0	32	133	45	1033	This study
Sulfolobus solfataricus	NC_002754	glucose; cellobiose; minimal	3034	1094	0	124	202	109	1040	(4)
Thermus thermophilus	NC_005835	mid-log phase	2035	731	0	109	41	15	688	This study

gTSS = A TSS found upstream to coding gene.

iTSS = internal TSS within a gene.

aTSS = cis-antisense TSS overlapping a gene.

nTSS = non-coding, intergenic TSS.

aIncluding overlapping operons.

Determination and assessment of reliability of TSSs

A machine learning approach was utilized to infer TSSs from all 5′ mapped sites (the set of putative TSSs). The inference algorithm consists of two main steps: (i) Heuristic determination of the training set and (ii) supervised Random Forest learning (pipeline steps in Supplementary Figure S1). The learning is based on a set of 15 features which capture relevant TSS characteristics including genomic, whole-transcriptome RNA-seq coverage patterns, and 5' differential RNA-seq (TAP(+)/(−), see list in Supplementary Table S2). A Random Forest machine learning procedure was used to infer and quantify reliability of the TSS determination. Finally, the classification of TSS types (e.g. gTSS as gene-related TSS or nTSS as non-coding TSS) was performed automatically based on genomic and transcriptomic features (Supplementary Extended Methods).

Inference of operons, sRNAs and putative ORFs

The general approach to infer operons (transcriptional units, TUs), sRNAs and putative ORFs was as previously described (19,29). Briefly, operon beginnings are defined by location of primary gTSSs while operon ends are defined by the genomic coordinates of the last gene which is likely to be expressed within the same TU. Inclusion of downstream genes in the same TU is dependent on both genomic and transcriptomic features (Supplementary Extended Methods). sRNAs inference is based on the presence of a TSS and expression in intergenic regions (Supplementary Extended Methods). Putative ORFs are inferred based on bioinformatic detection of putative start and stop codons in intergenic regions. Putative ORFs of at least 30 amino acids were subjected to BLASTX search against the non-redundant protein database.

The genomic browser and comparative analysis of gene families

The genomic browser is based on GenBank annotation for each organism in this study. The genomic browser annotation is augmented by sRNA and cis-regulators if found in Rfam repositories (30) or found with covariance models (31) homology search against the Rfam database (E-value < 10−3).

The comparative approach is focused on protein families. We constructed the set of homologous proteins by first producing the network of all homologies with blastp E-value < 10−3 and delineation of best bi-directional hits (BBH) between homologs. For readability we will refer to sets of BBH homologs genes as orthologs throughout the manuscript. For each locus across all organisms in our study all BBH loci are listed and a comparative viewer with aligned transcriptomes for the entire set is available. Further comparative analysis of gene families is facilitated by depicting, for each locus, the entire list of homologs (orthologs and paralogs) found in our compendium with the same COG (Clusters of Orthologous Groups) affiliation (32,33). The COG affiliation is shown separately for three phylogenetic levels: ‘TOL’, ‘Domain’ and ‘Phylum’ for orthologous groups across the entire Tree Of Life, Domain-specific, or Phylum-specific, respectively.

RESULTS

General approach and the TCOL transcriptome browser

To facilitate comparison of transcriptomes of diverged organisms we use standardized sequencing protocols and a fully automated computational pipeline for the reconstruction of transcriptome maps, including accurate inference of TSSs, and transcription-dependent prediction of sRNAs and operons. We produced a user-friendly webserver that visualizes the transcriptome maps, available at http://exploration.weizmann.ac.il/TCOL. The transcriptome browser provides a multi-track view of homologous genes across the tree of life, providing a unique comparative view into the transcriptional activity of gene families across diverged organisms.

The compendium of transcriptomes in this study encompasses the two prokaryotic domains of life (bacteria and archaea) with 18 representative species across 10 different phyla/subphyla (Figure 1). The transcriptomes of multiple organisms, typically grown in one or two conditions, were sequenced at a high depth both by standard Illumina RNA-seq, which captures the full length of the transcript, and using a 5′-end-specific sequencing protocol, which we previously described (4,19) (Table 1; Supplementary Table S1). For organisms analyzed in more than one growth condition, the transcriptome map for the comparative view is derived from merging the activity at the different conditions, resulting with greater sequencing depth and more comprehensive transcriptome maps. However, data tracks for the individual conditions are also presented in the TCOL browser, allowing the study of condition-specific transcription.

Figure 1.

Representative transcriptomes presented on the tree of life. The compendium of transcriptomes in this study is denoted as colored lineages across the prokaryotic tree of life consisting of 18 representative species across different phyla including Gamma, Alpha and Delta Proteobacteria, Cyanobacteria, Firmicutes, Deinococcus-Thermus, Actinobacteira, Spirochaetes and the Archael Crenarchaeotes phylum. Tree was taken from (44,45), reconstructed with (45).

Reconstruction of transcriptome maps

A quantitative summary of all transcriptome maps in this study is presented in Table 1. Primary TSSs were inferred using a machine learning approach, integrating relevant genomic and transcriptomic information in an unbiased manner using multiple features (Supplementary Table S2). In order to estimate the accuracy of our automated inference of TSSs, we compared these to a set of manually curated TSSs in three organisms for which published TSS maps are available - Listeria monocytogenes (19), Pseudomonas aeruginosa (29) and Sulfolobus solfataricus (4). In comparison with these benchmarks, both the sensitivity and precision were high, with sensitivity (True Positive, TP) of 79.3%, 74.5% and 64.1%, and precision (TP/TP+FP) of 78.3%, 70.9% and 64.1% for Listeria, Pseudomonas and Sulfolobus, respectively (Supplementary Materials Figure S2). Further inspection of a subset of TSSs in which the automatic annotation and the benchmark, manually curated data, disagreed showed that in the vast majority of cases the automatic inference is more likely to report the correct TSS (Supplementary Materials).

Using the TCOL transcriptome viewer

Comprehensive documentation of the functions available in TCOL is detailed in the ‘help’ page of the web server.

Figure 2 illustrates some of the main features of TCOL. Users can start by choosing a representative organism (Figure 2A). The next step is choosing a desired sample (e.g. growth condition), and pressing ‘Show’ to open the transcriptome viewer (Figure 2C). A user can navigate the browser by genome position, gene locus name or annotated Rfam genes. Genomic and transcriptomic information appears as ‘tooltips’ upon hovering with the cursor, including TSS information (e.g. hover over plotted TSS arrow to reveal position, number of supporting reads and reliability score of TSS), coverage information (number of RNA-seq reads per position), information on genes and regulatory elements and operon inferences. The user may click on a gene to open a ‘Gene page’ (Figure 2D) containing specific gene and homologs information, including links to the homologs transcriptional maps, and a link to the comparative viewer that visualizes all homologs together (Figure 2E; Figure 3). Furthermore, in each ‘Gene page’ the ‘COG-based homology’ section contains links to the transcriptomes of all the genes in our data set with the same COG (Clusters of Orthologous Groups) (32,33).

Figure 2.

The Transcriptomes Compared across the Tree of Life (TCOL) online web-server. Interactive browser with detailed transcriptional maps of 18 representative species across the microbial tree of life. Transcriptome maps can be accessed either by (A) choosing an organism or (B) by using BLAST with a user-provided query sequence. (C) Transcriptome maps can be navigated with multiple functions and are information-rich. Some of the information appears as tooltips when hovering above the features with the cursor. (D) ‘Gene page’ (when clicking on a gene in the browser) containing information about the gene including its homologs. (E) The comparative transcriptome viewer shows aligned transcriptomes of homologous genes, centered on the queried gene. This view is accessible from the ‘Gene page’ and from the BLAST results page.

Figure 3.

Comparative transcriptome views of homologous genes with long 5'UTRs. (A) Threonyl-tRNA synthetase (thrS) has a known T-box (RF00230) in Listeria, Clostridium and Bacillus. (B) The transcription termination factor rho gene (COG1158) in which regulatory structure was previously inferred in Pseudomonas (Pseudomon-Rho). (C) The cytochrome d ubiquinol oxidase subunit I (COG1271) gene (D) The cytochrome o ubiquinol oxidase subunit II (COG1622) gene. Brown horizontal arrows denote RFAM-identified cis-regulators.

Alternatively, browsing can start by BLAST searching a query sequence against the database of the TCOL compendium (Figure 2B). Each of the returning hits includes links to both the comparative transcriptome map as well as the maps for individual species. There are three different BLAST options: blastp, blastn and blastx.

Users interested in browsing all TSS or RNA-seq data globally can use the ‘Supplementary files’ tab that contains downloadable detailed information on the transcriptome annotation, including summary statistics (‘TSS statistics’), details of the TSS for each gene (‘TSS’), summary of total RNA-seq data statistics (‘RNA-seq summary’), number of reads mapped to each gene (‘Expression per gene’) and a summary file with additional transcriptome features including RFAM genes, operons and identified sRNAs (‘Additional transcriptome annotations’).

Comparative transcriptomics highlights cis-acting RNA regulators

One of the main applications of the comparative transcriptomics approach is to highlight cases in which recurring transcriptomic features are found across homologs, suggestive of regulatory importance and evolutionary conservation. We demonstrate this by focusing on orthologous genes showing long 5′ UTRs across the tree of life.

Long regulatory 5′ UTRs have the capacity to integrate environmental and inner-cellular cues into regulation of gene expression (34,35). Previous comparative genomics studies significantly enhanced the capacity to detect putative functional 5′ regulatory elements (36,37). However, it is likely that many additional regulatory 5′ UTRs elements remain to be found, including regulatory elements with conservation that is insufficient for reliable detection from genomic data alone (38,39).

Searching for gene families with long conserved 5′ UTR in multiple organisms across the tree of life, we found 41 gene families (COGs) with exceptional propensity for recurring long 5′ UTRs (Supplementary Table S3). Approximately three quarters of the gene families within this list recapitulated previously reported 5′ UTR regulatory elements (Supplementary Figure S4). These include 14 known riboswitches (including T-box, S-box, TPP, Cobalamin, FMN, glmS ribozyme and Glycine, Purine and Lysine riboswitches); 8 ribosomal leaders; 2 thermosensor RNA elements; and 6 cis-regulatory elements of various categories that had been reported including the rimP leader, Pseudomon-GroES, Pseudomon-Rho, Lacto-rpoB, mini-ykkC and leucine operon leader. We further found 8 novel putative candidates, inferred to harbor RNA-regulatory elements of unknown function in their 5′ UTRs, and additional 3 ribosomal leader candidates (Supplementary Table S3; Figure 3).

As examples, Figure 3 presents the comparative transcriptome viewer of four gene families with recurring long 5′ UTRs. Figure 3A illustrates the 5′ UTR of the threonyl-tRNA synthetase gene family (COG0441, thrS), where T-box 5′ UTR riboswitches are found in gram-positive bacteria (40). Among the 5 presented cases, the Rfam T-box leader (RF00230) is described only in Listeria, Clostridium and Bacillus (30). However, searching for structural homology (31), we found a similar structural RNA element in the 5′ UTRs of the gram-positive Lactobacillus and Lactococcus, suggesting that these Lactobacillales homologs are regulated by an element similar to the T-box riboswitch, but with structural homology that is too weak for detection based on sequence and/or structure prediction alone.

A regulatory element in the 5′UTR of the Rho transcription terminator was discovered bioinformatically in Pseudomonas aeruginosa, and termed Pseudomon-Rho (37). Among the homologs of the ‘Rho, Transcription termination factor’ gene family (COG1158), we found recurring long 5′UTRs across multiple genera including Pseudomonas, Gluconobacter, Desulfovibrio and Bacillus (Figure 3B). These results suggest that the Rho transcriptional terminator is being regulated by a leader 5′UTR element in multiple bacterial species in addition to Pseudomonas.

The gene families ‘cytochrome d ubiquinol oxidase subunit I’ (COG1271), and ‘cytochrome o ubiquinol oxidase subunit II’ (COG1622) are part of the CydABX and CyoABCD operons, respectively. These two different cytochrome oxidase complexes, responsible for cellular respiration, are known to be regulated by the oxygen state. We found, for both of these gene families, recurring long 5′UTRs across divergent bacteria, implying the existence of cis-acting RNA regulators within their 5′UTRs (Figure 3C and D).

CONCLUSION

The comparative transcriptomics approach was recently shown to be useful for studying closely related species – either strains within the same species (15) or species within the same genus (16,19,21). Broader evolutionary comparison across diverse species based on aggregating results from different studies (41) is limited due to methodological differences in RNA sequencing and computational analyses. Here we provide, for the first time, a comparative transcriptome viewer of representative species spanning the prokaryotic tree of life, generated by standardized methodologies in both data generation and data analysis. This constitutes an important step toward better understanding the similarities and differences in transcriptome regulation and dynamics in various organisms.