Temperature-dependent sRNA transcriptome of the Lyme disease spirochete

Background

Transmission of Borrelia burgdorferi from its tick vector to a vertebrate host requires extensive reprogramming of gene expression. Small regulatory RNAs (sRNA) have emerged in the last decade as important regulators of bacterial gene expression. Despite the widespread observation of sRNA-mediated gene regulation, only one sRNA has been characterized in the Lyme disease spirochete B. burgdorferi. We employed an sRNA-specific deep-sequencing approach to identify the small RNA transcriptome of B. burgdorferi at both 23 °C and 37 °C, which mimics in vitro the transmission from the tick vector to the mammalian host.

Results

We identified over 1000 sRNAs in B. burgdorferi revealing large amounts of antisense and intragenic sRNAs, as well as characteristic intergenic and 5′ UTR-associated sRNAs. A large fraction of the novel sRNAs (43%) are temperature-dependent and differentially expressed at the two temperatures, suggesting a role in gene regulation for adaptation during transmission. In addition, many genes important for maintenance of Borrelia during its enzootic cycle are associated with antisense RNAs or 5′ UTR sRNAs. RNA-seq data were validated for twenty-two of the sRNAs via Northern blot analyses.

Conclusions

Our study demonstrates that sRNAs are abundant and differentially expressed by environmental conditions suggesting that gene regulation via sRNAs is a common mechanism utilized in B. burgdorferi. In addition, the identification of antisense and intragenic sRNAs impacts the broadly used loss-of-function genetic approach used to study gene function and increases the coding potential of a small genome. To facilitate access to the analyzed RNA-seq data we have set-up a website at http://www.cibiv.at/~niko/bbdb/ that includes a UCSC browser track hub. By clicking on the respective link, researchers can interactively inspect the data in the UCSC genome browser (Kent et al., Genome Res 12:996-1006, 2002).

Lyme disease is the most common arthropod-borne disease in the United States of America and Europe, with 300,000 cases reported annually in the US [1]. Borrelia burgdorferi, a causative agent of Lyme disease, oscillates in nature between a tick vector and a vertebrate host [2–4]. The enzootic life cycle of B. burgdorferi is maintained by uninfected tick larvae feeding on an infected vertebrate, usually a small mammal. The infected larvae molt into nymphs and the spirochetes are transmitted to and infect a vertebrate at the next blood meal [4]. The infected nymphs are also the primary route of B. burgdorferi transmission to humans. B. burgdorferi must survive in and transition between two vastly different environments, the tick vector and the vertebrate host [4, 5]. Like many other pathogenic bacteria, B. burgdorferi senses and responds to environmental cues, such as a change in temperature [6–9], by regulating the gene expression of proteins necessary for survival [4, 5].

Bacterial gene expression is highly regulated at the level of transcription, which is catalyzed by RNA polymerase (RNAP) and regulated by transcription factors. The RNAP sigma factor is responsible for promoter selectivity. Many bacteria synthesize several different sigma factors, with different promoter selectivity, thus directing RNAP to a discrete set of genes, which results in the control of a set of genes needed for a certain response [10]. B. burgdorferi has only three sigma factors (RpoD, RpoS and RpoN), a relatively small number compared to other bacteria, which can encode up to eighteen [10]. Moreover, transcription of rpoS is regulated by RpoN, effectively decreasing the regulatory breadth of sigma factors in B. burgdorferi [11]. Transcription is also regulated via several characterized transcription factors in B. burgdorferi [5, 12–14]. However, little is known about post-transcriptional gene regulation in this spirochete.

Posttranscriptional gene regulation via a variety of regulatory RNAs has emerged in the past decade as a major mechanism of modulating gene expression [15–18]. The most extensively studied regulatory RNAs in bacteria are the trans-encoded small RNAs (sRNAs) that are predominately encoded in intergenic regions between two annotated genes. Most trans-encoded sRNAs regulate gene expression by imperfectly base pairing with a target mRNA and affecting the translation and/or the stability of the mRNA. The RNA chaperone Hfq is often required for sRNA:mRNA base-pairing [19]. Cis-encoded antisense RNAs (asRNAs), which are RNAs transcribed opposite to annotated genes, have been ubiquitously reported [20–22]. asRNAs have complete complementarity to their sense mRNA counterpart and were originally identified opposite to phage, toxin and transposon genes [23, 24]. asRNAs act in a manner similar to trans-encoded sRNAs via binding their cognate mRNA and influencing its translation and/or stability [20, 21]. However, asRNAs can also regulate their sense mRNAs transcription via transcriptional interference [20, 21, 24, 25]. Riboswitches are another cis-acting class of regulatory RNAs encoded in long 5′ UTRs of the mRNAs they regulate. Riboswitches usually function as sensors of metabolic cues or temperature changes [26–28]. The structure of the riboswitch is affected either by binding of a metabolite or changes in temperature stimulating or inhibiting transcription or translation of the gene. Non-coding RNAs can also act as regulators by binding proteins and sequestering them or affecting their activity [29, 30]. Furthermore, intragenic sRNAs are a new class of sRNAs that are encoded within annotated genes [31]; relatively little is known about their function. Finally, dual-function RNAs encode both a regulatory RNA and a protein. Staphylococcal RNAIII regulates the translation of several genes via imperfect base pairing and encodes a small (25 amino acid) hemolysin peptide [32, 33].

sRNAs are recognized as important regulators of many adaptive and physiological gene expression changes in pathogenic bacteria [34–36]. However, despite the pervasive nature of regulatory RNAs, only one regulatory RNA has been identified and characterized to date in B. burgdorferi [37]. B. burgdorferi encodes two characterized RNA-binding proteins, a unique Hfq protein (Hfq_Bb) (39), which is required for murine infection via needle inoculation, and a homolog of CsrA (CsrA_Bb), although there is controversy regarding its function in infection of the mammalian host [38–41]. CsrA normally acts in a concerted manner with two non-coding RNAs, which have not been identified in B. burgdorferi. We hypothesized that B. burgdorferi has a large sRNA network that is required for transducing the enzootic life cycle and pathogenesis. Here we specifically identified the sRNA transcriptome of B. burgdorferi at 23 °C and 37 °C, temperatures that mimic the tick vector and vertebrate host, respectively, and we found a large sRNA network. This study is the first transcriptome-wide analysis of sRNAs in the Lyme disease spirochete.

Transcriptome-wide identification of small RNAs

The main goal of this study was to identify small RNAs in B. burgdorferi that are important for gene regulation associated with the enzootic cycle of the spirochete [4, 5]. Temperature is one of the key environmental stimuli that modulate gene expression in B. burgdorferi [6–8]. For this reason, we shifted the temperature of B. burgdorferi growing in liquid culture from 23 °C to 37 °C to mimic transmission from the tick vector to the vertebrate host and deep-sequenced the small RNA transcriptome from these cultures.

Stranded cDNA libraries were prepared from the ribosomal RNA-depleted size-selected (50–500 nt) RNA. Three independent biological replicates were sequenced. To obtain the most complete coverage of the B. burgdorferi genome and capture also lowly expressed sRNAs, the first biological replicates (23 °C and 37 °C) were each sequenced on a single lane of an Illumina HiSeq 2000 and mapped to the B. burgdorferi strain B31 reference genome (replicate 0, rep0). This resulted in 170 and 190 million mapped 50 bp single-end (SE) reads, respectively, which corresponds to very deep theoretical genomic coverages of about 5600X and 6200X. Two additional replicates (rep1, rep2) were sequenced to lower coverage (27–44 million mapped reads, 900-1400X, see Table 1) and were used for validation and differential gene expression analysis.

For the accurate identification of B. burgdorferi sRNAs, we extracted strand-specific coverage signals from our deep sequencing datasets and developed a simple computational method to search these data for sRNAs-derived coverage peaks (see Methods). This method identified peaks stemming potentially from sRNAs in both deep data sets (23 °C and 37 °C) that were then merged in order to get a final set of genomic candidate intervals (sRNAs) based on the peak borders, see Methods. This strategy and the depth of coverage of our data sets allowed us to identify sRNAs expressed at low levels as well as sRNAs expressed in a temperature-dependent fashion. We acknowledge, however, that our method may have missed some very lowly expressed sRNAs due to the thresholds we configured for reducing the number of false positive calls (e.g., a minimum peak height of 500 reads), see Methods. As proof of principle, our method identified known and annotated tRNAs, RNase P, tmRNA, 4.5S RNA, 6S RNA (S. Samuels personal communication) and DsrA_Bb. For example, coverage maps illustrate the peaks called for two phenylalanine tRNAs (Additional file 1: Figure S1).

The resulting list of 5,600 genomic intervals identified by our peak caller were then manually curated by multiple members of our research group by visual inspection of the normalized coverage data in the IGV genome browser [42] resulting in 1,005 putative sRNAs (Additional file 2: Table S1) that were classified based on their genomic context (Fig. 1a). 116 of the called genomic intervals overlapped the 5′ end of an annotated open reading frame and were classified as a 5′ UTR sRNAs. 156 of sRNAs were transcribed between annotated open reading frames (ORFs) and were classified at intergenic sRNAs (IG-sRNAs). 357 sRNAs were found to be transcribed from the opposite strand to an annotated open reading frame and were classified as cis-encoded antisense RNAs (asRNA). 339 sRNAs were identified within annotated ORFs and were classified as intragenic sRNA (intraRNAs). The remaining 37 sRNAs were previously known and annotated accordingly (tRNAs, 5S RNAs, ffx (4.5S RNA), SsrA, DsrA and RNase P). For the categorization of the sRNAs, we took a simplified approach (the as-, intra- and 5′ UTR RNAs only needed one base of the interval to overlap with the open reading frame) and acknowledge that, depending on the genomic context, intraRNAs, IG RNAs and asRNAs may be 5′ UTRs for proximal genes; however, this cannot be determined from our sRNA sequencing. In summary, our search strategy successfully detected all 37 previously annotated sRNAs and additionally identified 968 novel sRNAs. In addition, an in silico analysis of the nucleotide context surrounding the detected 5′ ends of the sRNAs identified an adenine and thymine rich −10 region correlating to the Pribnow box providing additional evidence for the accuracy of our approach (Additional file 3: Figure S7).

Genomic location, categorization and temperature-dependence of small RNAs. a sRNAs were categorized according to their genomic context as 5′ UTR sRNAs (blue arrow), antisense (AS) RNAs (red arrow), intergenic (IG) RNAs (yellow arrows), or intragenic (intra) RNAs (green arrow). Annotated coding features (ORFs) are represented as black arrows. The numbers of sRNAs identified in each category are in parenthesis next to the category. b A bar chart illustrates the numbers of temperature-dependent and -independent sRNAs identified in each category. Gray bars are temperature-independent sRNAs, blue bars are sRNAs up-regulated at 23 °C compared to 37 °C and red bars are sRNAs up-regulated at 37 °C compared to 23 °C

Full size image

RNA-seq results including normalized coverage tracks and all identified sRNAs are available at http://www.cibiv.at/~niko/bbdb/. The data can be downloaded or directly viewed in the UCSC Genome Browser [43]. Note that the normalization method used to generate the coverage tracks for visualization differs from the method used by differential gene expression analysis programs (EdgeR and DESeq). Therefore, visual inspection of sRNA coverage signals cannot be used to determine statistically significant gene expression changes.

Temperature-dependent sRNAs

We identified the sRNA transcriptome at both 23 °C and 37 °C to elucidate sRNAs involved in temperature-dependent gene regulation associated with the enzootic cycle. Stranded sRNA transcriptomes from the two biological low-coverage replicates at both temperatures were sequenced and annotated sRNAs from our peak-finding algorithm were used to test for statistically significant differential expression at the two temperatures using EdgeR. This method identified 431 (43%) of the sRNAs as temperature-dependent and differentially regulated at the two temperatures: 128 sRNAs were up-regulated at 23 °C, while 303 were up-regulated at 37 °C. 22 sRNAs were validated by Northern blot analyses (Additional file 4: Table S2 and Additional file 5: Figure S2, Additional file 6: Figure S3, Additional file 7: Figure S4, Additional file 8: Figure S5 and Additional file 9: Figure S6) and temperature-dependent and temperature-independent sRNAs were identified in all categories of sRNAs (Fig. 1b).

Cis-encoded small RNAs

asRNAs are transcribed opposite to annotated ORFs and have some portion completely complementary to the corresponding sense mRNA. asRNAs vary in size from a few hundred nucleotides to 6.5 kb [20, 21]. We identified 357 small asRNAs: 147 are temperature-dependent; 25 are up-regulated at 23 °C, while 122 are up-regulated at 37 °C (Fig. 1b). Overall, 296 out of 1569 annotated genes have at least one asRNAs encoded opposite to them and 53 of these have more than one associated asRNA (Additional file 10: Table S3). asRNA-dependent gene regulation occurs via base pairing with its cognate sense mRNA and inhibiting or stimulating expression via different mechanisms. Predominately, genes with asRNAs associated with them are hypothetical ORFs with unknown functions (Fig. 2a). Notably, asRNAs were identified opposite to key genes active in the maintenance of the B. burgdorferi life cycle and cell division, such as: bb0420 (hk1), bb0827 (hrpA), bbb03 (resT), bbb17 (guaB), and bba66 [4, 13, 44–50]. Hk1 is the histidine kinase in the two-component system necessary for B. burgdorferi survival in the tick [45]. There are two as-hk1 RNAs, which are both up-regulated at 37 °C. The hk1 transcript is expressed at low levels in vivo in both the flat tick, fed tick and in B. burgdorferi grown in dialysis membrane chambers (mimicking the mammalian host). However, the levels were the highest in larvae fed to repletion and flat ticks [45]. We hypothesize that the as-hk1 RNAs play a role in the fine-tuning of gene regulation of hk1. There are also two as-resT RNAs, one is temperature independent and the other is up-regulated at 23 °C. ResT is the telomere resolvase required for telomere resolution during replication of the linear and circular genetic elements [51]. Finally, there are two asRNAs, both up-regulated at 37 °C, encoded opposite to bba66, which is required for mouse infection via tick transmission [44].

Identification and confirmation of small asRNAs. a Genes with asRNAs encoded opposite to them are functionally categorized using the following abbreviations: CD, cell division; CE, cell envelope; CM, cell motility; MT, metabolism; PD, protein degradation and folding; SR, stress response; TL, translation; TP, transporter; TR, transcription; RM; RNA metabolism; DM, DNA metabolism; VM, vitamin metabolism and U, unknown. asRNAs (SR0493 and SR0104) identified by sRNA-seq were confirmed via Northern blot analyses for as-bb0155 and as-bb0612. b and d The deep sequencing results are displayed in a coverage map of an overlay of the two rep0 libraries sequenced deeply for sRNA peak calling (Peak) and the overlay of the two biological replicates rep1, rep2 at both 23 °C and 37 °C. The height at each position indicates the normalized number of reads that mapped to that base. The – strand coverage is shown in blue. Note that the y-axis scale is different between the peak calling libraries (peak) and the biological replicates used for differential expression analyses (indicated by the grey numbers in brackets). The genomic context is illustrated below the coverage maps: black arrows indicate the annotated ORFs, the yellow box indicates the region called as a small asRNA by our peak finder and the wavy line is the putative transcript determined by Northern blot. The red line represents the location of the oligonucleotide probes used for the Northern blot. c and e Northern blot analyses of total RNA fractionated on a denaturing polyacrylamide gel, blotted to a nylon membrane, and hybridized with oligonucleotides. The predicted transcripts are denoted and marked with the appropriate band (^ for sRNA) in the Northern blot. Two independent experiments were performed and representative data are shown

Full size image

In 2002, microarray analysis identified 79 genes up-regulated and 15 genes down-regulated at 37 °C compared to 23 °C [52]; 17 of these 94 genes have asRNAs associated with them. Eight of these asRNAs are temperature-dependent in our data; two of the sense/as RNA pairs are reciprocally regulated, while six are co-regulated (Additional file 11: Table S4). These data hint at a molecular mechanism for gene regulation by the asRNAs. For example, bb0385 is down-regulated at 37 °C compared to 23 °C in the microarray data [52], while the as-bb0385 RNA (as-0385) is up-regulated at 37 °C compared to 23 °C in our RNA-seq data. We hypothesize that the asRNA, when up-regulated at 37 °C, binds to the mRNA and inhibits transcription or initiates degradation effectively lowering the levels of the mRNA. For the mRNAs and asRNAs that are up or down-regulated together, we propose that the asRNA stabilizes the mRNA and/or stimulates transcription of the gene, or acts in trans on a different mRNA.

RNA-seq data were validated by Northern blot analyses of five of the asRNAs (Additional file 4: Table S2 and Additional file 5: Figure S2). We identified asRNAs encoded opposite to all regions of their corresponding mRNAs; for example, the bb0612 (clpX) gene has an asRNA transcribed opposite to its 5′ end and the hypothetical open reading frame bb0155 has an antisense RNA encoded opposite to it’s 3′ end. Northern blot analyses validate the RNA-seq data and demonstrate that both as-0155 and as-0612 RNA steady-state levels are higher at 37 °C compared to 23 °C (Fig. 2b–d).

Stable small RNAs originating from 5′ and 3′ UTRs have been reported in several different pathogenic bacteria [18, 31, 53–59]. Cis-acting transcriptional riboswitches in 5′ UTRs terminate transcription in the absence or presence of metabolites and regulate gene expression of the associated mRNA. Early termination of transcription results in small RNA byproducts, which can act as trans-encoded sRNAs regulating gene expression of another mRNA. We identified 116 sRNAs associated with the 5′ UTR’s of annotated ORFs, suggesting they may act as riboswitches and/or trans-encoded sRNAs. Riboswitches often control expression of genes for transport or synthesis of key metabolic compounds and genes involved in physiological changes, virulence and stress responses. The majority of the genes with 5′ UTR sRNAs associated with them are hypothetical ORFs of unknown function and genes associated with metabolism and RNA and DNA metabolism (Fig. 3a). Of note, the transcriptional regulator BosR [14] has a 5′ UTR up-regulated at 37 °C. 44 of the 5′ UTR sRNAs were differentially regulated by temperature, suggesting these RNA elements may act as a type of transcriptional thermosensor, riboswitches that respond to temperature changes, rather than a ligand. To validate the RNA-sequencing data, we performed Northern blot analyses on two of the 5′ UTR sRNAs (Additional file 2: Table S1 and Additional file 6: Figure S3) Northern blot analysis of the 5′ UTR sRNA associated with bba57 confirms the presence of an sRNA at 37 °C and mRNA at both 23 °C and 37 °C (Fig. 3b and c). We hypothesize that this sRNA may act in-trans to regulate a different mRNA or could be a transcriptional RNA thermometer. RNA thermometers are highly structured RNAs in the 5′ UTR of mRNAs that usually mask the ribosome-binding site at low temperatures inhibiting translation. At higher temperatures the structure melts and the ribosome-binding site is accessible allowing for efficient translation [60, 61]. To our knowledge no transcriptional RNA thermometers have been described. We strictly categorized 5′ UTR sRNAs, requiring the intervals we call in this category to overlap with the 5′ end of an annotated ORF. However, some of the intergenic, intragenic or antisense sRNAs may be associated with the 5′ UTRs of downstream genes and still function as riboswitches.

Identification and verification of 5′ UTR small RNAs. a Genes with small RNAs encoded in their 5′ UTRs are functionally categorized using the following abbreviations: CD, cell division; CE, cell envelope; CM, cell motility; MT, metabolism; PD, protein degradation and folding; SR, stress response; TL, translation; TR, transcription; RM; RNA metabolism; DM, DNA metabolism; VM, vitamin metabolism and U, unknown. b A 5′ UTR small RNA (SR0902) was confirmed via Northern blot analyses for bba57. The deep sequencing results are displayed in a coverage map as described in the Fig. 2 legend, except the yellow box indicates the region called as a small 5′ UTR RNA by our peak finder. c Northern blot analyses of total RNA fractionated on a denaturing polyacrylamide gel, blotted to a nylon membrane, and hybridized with oligonucleotides. The predicted transcripts are denoted and marked with the appropriate band (* for sRNA and ^ for mRNA) in the Northern blot. Two independent experiments were performed and representative data are shown

Full size image

Intragenic small RNAs

Intragenic small RNAs (intraRNAs) are a relatively new class of sRNAs encoded from within protein coding regions. We identified 339 intraRNAs encoded within 287 annotated genes. 243 genes have one intraRNA encoded within them, while 44 genes have more than intraRNA associated with them. Intragenic RNAs have been identified via genome-wide transcriptional start site identification or co-immunoprecipitations with Hfq in several different bacteria. The genome size of Helicobacter pylori is similar to B. burgdorferi and 439 internal transcriptional start sites were identified globally. Our data identified 143 temperature-dependent intraRNAs: 86 were up-regulated at 37 °C, while 57 were up-regulated at 23 °C (Fig. 1b). Several genes important for maintenance of B. burgdorferi in its enzootic cycle also encode intraRNAs including [4]: bb0382 (bmpA), bb0365, bb0419 and bb0420 (rrp1 and hk1, respectively), bb603, bbb18 (guaA), bbk17, bbk32, bba16 (ospB), bba64, and bba66. For instance, we identified a temperature-independent intraRNA encoded at the 3′ end of the hypothetical open reading frame bbq07 and confirmed it via a Northern blot (Fig. 4). bbq07 potentially encodes two small ORFs; one is a truncated version of the bbq07 ORF with a canonical AUG start codon, while the other would utilize an alternative start codon UUG and encode a different ORF. Both putative proteins would be small peptides, 27 and 33 amino acids, respectively. Relatively little is known about the function of intraRNAs and they could encode non-coding regulatory RNAs or mRNAs for small peptides. Most intraRNAs have been identified through genome-wide transcriptional start site analyses or co-immunoprecipitation with Hfq and few have been functionally characterized [31, 55, 59, 62]. However, the only characterized sRNA in B. burgdorferi, DsrA_Bb, is an intragenic RNA encoded from within bb0577 that post-transcriptionally regulates the alternative sigma factor RpoS [37]. Recently, a small RNA processed from the 3′ UTR of an mRNA in Salmonella was shown to act in trans to regulate another mRNA [54]. The identification of intraRNAs greatly increases the coding potential of the relatively small genome of B. burgdorferi.

Identification and confirmation of small intragenic RNAs. The intraRNA (SR0921) was confirmed via Northern blot analyses for intra-bbq07. a The deep sequencing results are displayed in a coverage map as described in Fig. 2, except the yellow box indicates the region called as a small intraRNA by our peak finder. Putative ORFs are indicated by black lines with the canonical start codon (AUG) illustrated by a green box, the non-canonical start codon illustrated by the green box with a black outline, and stop codons indicated by a red box. b Northern blot analyses of total RNA fractionated on a denaturing polyacrylamide gel, blotted to a nylon membrane, and hybridized with oligonucleotides. The predicted transcripts are denoted and marked with the appropriate band (* for sRNA) in the Northern blot. Two independent experiments were performed and representative data are shown

Full size image

Intergenic sRNAs

The most well-studied class of sRNAs are encoded intergenically (IG-sRNA) and act in trans by base pairing imperfectly with an mRNA, encoded from another genomic region, and regulating its expression. Currently the number of IG-sRNAs is approximately 300 in both E. coli and Salmonella enterica. Moreover, in Yersinia pseudotuberculosis the sRNA transcriptome identified 150 novel intergenic sRNAs. Here, we identified 156 intergenic sRNAs; 42 are up-regulated at 37 °C and 38 are up-regulated at 23 °C (Fig. 1b). Considering the small genome size of B. burgdorferi, 156 IG-sRNAs seems on par with other well-studied bacteria. Northern blot analyses validate 13 IG-sRNAs (Additional file 4: Table S2, Additional file 8: Figure S5 and Additional file 9: Figure S6). The sRNA encoded between bbb13 and bbb14 is not differentially expressed at 23 °C or 37 °C in the RNA-seq data and has similar steady-state levels at both temperatures (Fig. 5a and b). In contrast, the IG-sRNA encoded between bba34 and bba36 is up-regulated at 37 °C in the RNA-seq data and has higher steady-state levels at 37 °C in the Northern blot analyses (Fig. 5c and d). Like other intergenically encoded sRNAs, we hypothesize that many of the IG-sRNAs regulate their target mRNAs via a variety of mechanisms requiring imperfect base pairing between the sRNA and the target mRNA. However, these sRNAs could also code for peptides. For example, the sRNA encoded between bbb13 and bbb14 has a putative small ORF (28 amino acids) (Fig. 5a).

Identification and verification of intergenic small RNAs. IG-sRNAs (SR0961 and SR0897) were confirmed by Northern blot analyses. a and c The deep-sequencing results for the IG-sRNAs between bbb13 and bbb14 and between bba34 and bba36 are displayed in a coverage map as described in Fig. 2. A putative ORF is indicated by a black line with a canonical start codon (AUG) illustrated by a green box and the stop codon indicated by a red box. b and d Northern blot analyses of total RNA fractionated on a denaturing polyacrylamide gel, blotted to a nylon membrane, and hybridized with oligonucleotides. The predicted transcripts are denoted and marked with the appropriate band (*,^ for sRNA) in the Northern blot. Two independent experiments were performed and representative data are shown

Full size image

Small open reading frames

sRNAs are considered primarily non-coding and function as riboregulators. However, there are several reports of dual-function sRNAs that act as a riboregulator and code for a small peptide. The sRNA RNAIII in Staphylococcus aureus codes for a regulatory RNA and a 26 amino acid peptide [32, 33]. In addition, small ORFs are emerging as important components of cells in both bacteria and eukaryotes [63, 64]. We examined the coding-potential of all of the 1005 sRNAs we identified, searching for any ORFs at least 25 amino acids in length using canonical start codons. 988 of the sRNAs putatively encode at least one 25 amino acid ORF. Although this is not an accurate indicator of protein-coding sRNAs, it does demonstrate the capacity of B. burgdorferi sRNAs to code for small peptides.

sRNAs are now universally considered to be major regulators of gene expression in pathogenic bacteria [18]. However, almost nothing was known about sRNAs in the Lyme disease spirochete, B. burgdorferi [37, 39]. Here we report a large temperature-dependent and -independent sRNA transcriptome in B. burgdorferi. The identification of massive antisense and intragenic sRNAs impacts the broadly used genetic approach to studying gene function in B. burgdorferi. There are now profound implications for mutagenesis experiments: deletion of a gene could include the loss of an intraRNA and/or asRNA with pleiotropic consequences. The abundance of sRNAs in the B. burgdorferi and their association with genes required for maintenance of its enzootic cycle suggest post-transcriptional gene regulation plays a significant role in pathogenesis.

Bacterial strains

Low passage B. burgdorferi strain B31-5A4 was grown at 23 °C and shifted to either 23 °C or 37 °C in BSK-H medium (Sigma) and grown to a low cell density 2 to 5 × 10⁷ cells/mL. Infectious B. burgdorferi strain B31-5A4 is a clone of the B31 strain used in Revel et al. Linear plasmid 5 was the only plasmid missing from the B31-5A4 strain used in this study.

RNA isolation and library preparation

Total RNA was isolated using a hot phenol protocol. Total RNA was treated with DNase I (Roche) following the manufacturer’s protocol. RNA integrity was measured using the Agilent 2100 Bioanalyzer. RNA with an RNA Integrity Number (RIN) above 9.0 was used for cDNA library construction. Directional (strand-specific) RNA-seq cDNA libraries were constructed with a ligation based protocol as previously described with an initial size-selection instead of fragmentation [22, 65]. Briefly, ribosome-depleted (Ribo-Zero RNA removal kit for gram negative bacteria; Epicentre) total RNA was size fractionated on an 8% TBE-UREA gel. RNA was eluted from gel slices correlating to 50 to 500 nucleotides. RNA was treated sequentially with tobacco acid pyrophosphatase (Epicenter) and calf intestinal phosphatase (New England Biolabs) per the manufacturer’s protocols to remove 5′ tri- and monophosphates. A 3′ RNA adaptor, based on the Illumina multiplexing adaptor sequence (Oligonucleotide sequences © 2007–2014 Illumina, Inc. All rights reserved) blocked at the 3′ end with an inverted dT (5′-GAUCGGAAGA GCACACGUCU [idT]-3′), was phosphorylated at the 5′ end using T4 PNK (New England Biolabs) per the manufacturer’s protocol. The 3′ multiplex RNA adaptor was ligated to the 3′ ends of the RNA using T4 RNA ligase I (New England Biolabs). RNA was incubated at 20 °C for 6 h in 1X T4 RNA ligase reaction buffer with 1 mM ATP, 20 μM 3′ RNA adaptor, 1 μl DMSO, 5 U of T4 RNA ligase I, and 40 U of RNasin (Promega) in a 10 μl reaction. RNA was gel purified and size-selected (75–550 nt) and purified over a denaturing 8% polyacrylamide/8 M urea/TBE gel. Gel slices were incubated in RNA elution buffer (10 mM Tris–HCl, pH 7.5, 2 mM EDTA, 0.1% SDS, 0.3 M NaOAc) with vigorous shaking at 4 °C overnight. The supernatant was subsequently ethanol precipitated using glycogen as a carrier molecule. The RNAs were phosphorylated at the 5′ ends using T4 PNK (New England Biolabs) per the manufacturer’s protocol to allow for subsequent ligation of the 5′ RNA adaptor. The Illumina small RNA 5′ adaptor (5′-GUUCAGAGUU CUACAGUCCG ACGAUC-3′) was ligated to the libraries using T4 RNA ligase I (New England Biolabs). The ligated RNAs were size-selected (100–600 nt) and gel-purified over a denaturing 8% polyacrylamide/8 M urea/TBE gel (as described above). The di-tagged RNA libraries were reverse-transcribed with SuperScript®II reverse transcriptase (Invitrogen) using random nonamers per the manufacturer’s protocol. cDNA was amplified in PCR performed using Phusion® High-Fidelity Polymerase (New England Biolabs). cDNA was amplified with Illumina-compatible PCR primers by 18 cycles of PCR. The products were analyzed on an Agilent 2100 Bioanalyzer.

Deep-seq analyses

Three independent biological replicates were sequenced on an Illumina HiSeq2000 (single-end, read length: 50 bp) to different genomic coverages (Table 1). The high-coverage replicate (rep0) was used for the identification of B. burgdorferi sRNAs, the two low-coverage replicates (rep1, rep2) were used to identify sRNAs that are differentially expressed at the two probed temperatures. Sequencing adapters were removed from the data using cutadapt v1.2.1 [66]. The resulting reads were mapped to the B. burgdorferi B31 genome (GenBank Ids: AE000783, AE001583, AE000793, AE001582, AE000785, AE000794, AE000786, AE000784, AE000789, AE000788, AE000787, AE000790, AE001584, AE000791, AE000792, AE001575, AE001576, AE001577, AE001578, AE001579, AE001580, and AE001581) using NextGenMap v 0.4.5 with default parameters and minimum identity set to 90%. Finally, multi-mapped reads were removed and strand-specific coverage data was extracted using bedtools v2.25 [67]. We developed a simple peak-calling algorithm for identifying potential sRNA peaks in these coverage signals. Briefly, the strand-specific coverage signal is read and smoothed using a moving Gaussian kernel. Then, peaks are called with a wavelet-transform based approach that uses a simple mirrored sawtooth kernel for the detection of potential peak boundaries. Peaks were called only if they complied to predefined minimum/maximum dimensions, based on our size selection (peak width: 45–500 bp, minimum peak height: 500 reads). Additionally, peaks were filtered based on their shape as we expected sRNA-derived peaks to present a sharp rise at the 5′ end and an overall “boxy” shape (i.e., a minimum height/width ratio), as seen with the positive control tRNAs (Additional file 1: Figure S1). sRNA libraries were generated from non-fragmented, size-selected RNAs and sequenced from a single end; therefore, we expected and observed a strong 5′ end bias in the sRNA coverage, which we took advantage of to determine the 5′ end of the sRNA. However, coverage of the entire sRNA is unlikely and due to capturing the degradation products of RNAs, thus making it difficult to accurately and consistently determine the 3′ end of the sRNAs. Our peak-calling algorithm, together with manual curation accurately identified intervals that correlate to sRNAs and their 5′ ends, but not their 3′ ends (Additional file 3: Figure S7). In this manner, peaks were called in both deep data sets (rep0, 23 °C and 37 °C) and corresponding genomic intervals were merged between data sets and between intervals if they overlapped more than 80%.

The resulting list of 5,600 genomic candidate intervals (sRNAs) was then subjected to differential expression (DE) analysis. We extracted count tables for these intervals from the two low-coverage replicates (rep1, rep2) and calculated DE with edgeR and DEseq [68, 69]. 5 intervals could not be tested due to too-low coverage in (one of) the replicates. Results were then filtered by adjusted P-value ≤ 0.05 (P-values were adjusted using Benjamini and Hochberg’s algorithm to control the false discovery rate). For further analyses we continued with only the edgeR data. We do, however, provide the adjusted DEseq P-values in the final result table for completeness.

Finally, all genomic intervals were manually assessed and curated by multiple members of our research group by visual inspection of the normalized coverage data in the IGV [70] genome browser. Besides obvious false positive calls and likely degradation products (Additional file 12: Figure S8), we also removed all peaks associated with annotated open reading frames that were 600 nucleotides or shorter from the list of sRNAs.

Analysis of repetitive regions

The B. burgdorferi genome consists of a linear chromosome and up to 21 or more linear (lp) as well as circular (cp) plasmids. A considerable fraction of the plasmid DNA is highly repetitive, particularly some regions on the cp32s, lp56, lp21 and lp5 are nearly sequence-identical [71], which makes unambiguous read mapping in these regions difficult or impossible. This constituted two major problems for our analysis: first, we observed artificial coverage peaks around the borders to such repetitive regions that were wrongly classified as sRNA peaks by our method, resulting in false-positive calls. Secondly, we were aware that this would also lead to false negative (i.e., missed) peaks in repetitive regions (as we removed all multi-mapped reads from our datasets).

For this reason, we conducted an additional analysis that guided our manual curation of sRNA candidate regions. First, we remapped all multi-mapped reads (i.e., reads with mapping quality zero (MQ0), mapping accurately to more than one place in the genome) from our alignments using NextGenMap, this time configuring the mapper to additionally output up to 100 alignments that share the maximum alignment score per read (i.e., a read that was sequenced from a genomic region that has three perfect copies in the genome would be represented with three entries in the resulting BAM file). Then, we calculated the read coverage signals from these data as we did before, but this time each alignment contributed to the coverage at a particular genomic position only with a weight 1/X0, where X0 is the number of optimal alignments for this read in the dataset. In other words, a read that aligned to three different (repetitive) regions with the same maximum score would contribute with a “weight” of 1/3 to the each of them, thereby equally dividing its contribution among all possible (optimal) mapping locations in the genome. The resulting coverage tracks were converted to the BigWig [72] format and loaded into IGV along with our standard tracks for use in our manual curation. A peak was considered a false positive (due to MQ0 reads) if it had reads surrounding it in the MQ0 coverage map. For example, an intragenic peak was called in the bbs02 gene on cp32-3, but manual inspection of the MQ0 reads shows good coverage of the proximal area around the peak, suggesting it is a false positive peak due to repetitive sequence surrounding it (Additional file 13: Figure S9). In addition, we likely missed sRNAs in these repetitive regions as well. For example, manual inspection of the MQ0 reads on cp32-1 revealed a putative antisense RNA opposite the bbp17 gene. The uniquely mapped coverage maps do not have any reads mapped to this genomic region because it is repetitive with the other cp32 plasmids, but the MQ0 coverage map suggests at least one of these cp32 derived sequences are transcribed (Additional file 14: Figure S10).

Northern blots

For Northern blot analysis 15 μg of DNase I (Roche) treated RNA was separated under denaturing conditions either by a 6–8% TBE-Urea (8 M) polyacrylamide gels in 1X TBE (for small transcripts) or a 1% formaldehyde/MOPS agarose gels in 1X MOPS (for larger transcripts). RNA was initially denatured in 2X RNA load dye (Fermentas) and heated to 65 °C for 15 min before loading on a gel. RNA was transferred to HybondXL membranes either by electroblotting at 12 V for 1 h in 0.5X TBE (polyacrylamide TBE-Urea gels) or capillary action (formaldehyde-agarose gels). The membranes were UV cross-linked and probed with DNA oligonucleotide probes (Additional file 15: Table S5). DNA oligonucleotide probes were end-labeled with gamma-32P [ATP] and T4 PNK (New England Biolabs) per the manufacturer’s protocol.

Nucleotide context analysis

To validate the accuracy of our 5′ end peak-calling method on a genome-wide scale we analyzed the nucleotide composition around all identified 5′ ends. We calculated averaged nucleotide fractions in genomic windows centered at the most 5′ positions of identified sRNAs and plotted the results in Additional file 3: Figure S7. Overall, these data demonstrate an adenine and thymine (A and T) rich sequence at the −10 region correlating to the Pribnow box and an enrichment in thymine exactly 35 nucleotides upstream of the putative transcriptional start sites (TSS) across all novel sRNA categories described in this paper (as, IG, Intra, 5′ UTR), suggesting indeed these 5′ ends are accurate. Moreover, we didn’t detect these sequence elements associated with the annotated group of sRNAs, which are primarily tRNAs. tRNAs are processed from primary transcripts to their active and stable form and should not have a Pribnow box at the −10 region. The same plots were generated for all other sRNA subcategories used in this manuscript (as, IG, Intra, 5′ UTR), as well as for all sRNAs on the two different strands individually, and the respective signals were similar to the overall signal (the −35 T-peak being slightly less prominent for intra-RNAs). We have also generated respective plots for randomly chosen positions in the genome as a control, which showed, as expected, no major deviations from the genomic average (data can be found at http://www.cibiv.at/~niko/bbdb/). Together, these data provide additional in silico evidence for the ability of our method to accurately detect 5′ ends of small RNAs on a genome-wide scale.

3′ UTR:

3′ untranslated region

5′ UTR:

5′ untranslated region

asRNA:

Antisense small RNA

IG-RNA:

Intergenic small RNA

intraRNA:

Intragenic small RNA

sRNA:

Small RNA

We thank Scott Samuels and members of the Lybecker and Schroeder laboratories for thoughtful and critical reading of the manuscript and useful discussions. We also thank Johanna Stranner for excellent technical assistance and Maarja Lepamets for implementing an early version of the peakfinder software. We also thank the NGS Facility of the VBCF for providing excellent sequencing support.

Funding

This research was supported by the Austrian Science Fund (FWF) (Grants FWF F4301, and F4308 to R.S.), the University of Vienna and the University of Colorado Colorado Springs. PR acknowledges support by the RNA-DK Biology (W1207-B09). The funders had no role in the study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Availability of data and materials

Sequences have been deposited at the National Center for Biotechnology Information Sequence Read Archive (Study accession number SRP078488: experiment accession numbers SRX1948239 and SRX1948242). In addition, RNA-seq results including normalized coverage data and genomic intervals of all identified small RNAs are available at http://www.cibiv.at/~niko/bbdb/. The data can either be downloaded or directly viewed in the UCSC Genome Browser [43].

Authors’ contributions

ML, IB, NP, PR and RS designed experiments. ML performed experiments. NP and PR performed bioinformatics analyses. ML, NP, IB, PR interpreted experimental results. ML and NP wrote the manuscript. All authors critiqued and edited the final manuscript. All authors read and approved the final manuscript.

Competing interests

The authors declare that they have no competing interests.

Consent for publication

Not applicable.

Ethics approval and consent to participate

Not applicable.

Affiliations

1. Wellcome Trust Centre for Human Genetics, University of Oxford, Oxford, UK

   • Niko Popitsch

2. Department of Biochemistry and Cell Biology, Max F. Perutz Laboratories, University of Vienna, Vienna, Austria

   • Ivana Bilusic
   •  & Renée Schroeder

3. Center for Integrative Bioinformatics Vienna, Max F Perutz Laboratories, Medical University of Vienna, University of Vienna, Vienna, Austria

   • Philipp Rescheneder

4. Department of Biology, University of Colorado Colorado Springs, Colorado Springs, CO, USA

   • Meghan Lybecker

Corresponding author

Correspondence to Meghan Lybecker.

Reprints and Permissions