Deep sequencing whole transcriptome exploration of the σE regulon in Neisseria meningitidis.

Bacteria live in an ever-changing environment and must alter protein expression promptly to adapt to these changes and survive. Specific response genes that are regulated by a subset of alternative σ(70)-like transcription factors have evolved in order to respond to this changing environment. Recently, we have described the existence of a σ(E) regulon including the anti-σ-factor MseR in the obligate human bacterial pathogen Neisseria meningitidis. To unravel the complete σ(E) regulon in N. meningitidis, we sequenced total RNA transcriptional content of wild type meningococci and compared it with that of mseR mutant cells (ΔmseR) in which σ(E) is highly expressed. Eleven coding genes and one non-coding gene were found to be differentially expressed between H44/76 wildtype and H44/76ΔmseR cells. Five of the 6 genes of the σ(E) operon, msrA/msrB, and the gene encoding a pepSY-associated TM helix family protein showed enhanced transcription, whilst aniA encoding a nitrite reductase and nspA encoding the vaccine candidate Neisserial surface protein A showed decreased transcription. Analysis of differential expression in IGRs showed enhanced transcription of a non-coding RNA molecule, identifying a σ(E) dependent small non-coding RNA. Together this constitutes the first complete exploration of an alternative σ-factor regulon in N. meningitidis. The results direct to a relatively small regulon indicative for a strictly defined response consistent with a relatively stable niche, the human throat, where N. meningitidis resides.

Introduction

Bacteria live in an ever-changing environment and must alter protein expression promptly to adapt to these changes and survive. Specific response genes that are regulated by a subset of the alternative σ70-like sigma factors have evolved in order to respond to this changing environment [1]. In reaction to specific external stimuli, these sigma factors recruit RNA polymerases to the appropriate response genes. The name extracytoplasmic function (ECF) sigma factors or σE factors refer to the fact that most of the genes under control of these σ factors encode proteins residing in the outer membrane or periplasmic space [2].

Recently, using a proteomic approach, we have described the existence of a σE regulon in the obligate human bacterial pathogen Neisseria meningitidis [3]. The σE operon so far identified encompasses an operon consisting of 6 genes (NMB2140-NMB2145) among which the gene encoding σE itself (NMB2144). In addition msrA/msrB (NMB0044) was also found to be subjected to regulation by σE, it encodes for methionine sulfoxide reductase, an enzyme repairing proteins exposed to reactive oxygen species (ROS) [4], [5]. Furthermore, the anti-σE factor MseR (for Meningococcal SigmaE Regulator; NMB2145) was identified. Deletion of mseR results in overexpression of the other 5 genes in the σE operon, thereby increasing the expression of σ itself as well as msrA/msrB [3]. The proteomic exploration of the σE regulon of meningococci and a microarray gene expression analysis for the closely related human pathogen Neisseria gonorrhoeae performed by others have revealed a surprisingly small regulon for the σE factor in Neisseria [6]. However, both methods have several inherent and technical limitations. Identification of proteins is limited to those that can be isolated efficiently, are within a certain size range resolvable with SDS-PAGE, and have to be sufficiently expressed to be both visible and able to be confidently identified using mass spectrometry. DNA microarrays are only able to detect products that can sufficiently hybridize with probes that have been created based on available genomic data, usually limited to in silico predicted open reading frames (ORFs) of protein coding genes in their sense orientation. Moreover, reliability and reproducibility remains a concern [7]–[9]. Both approaches are furthermore unsuitable to detect an important class of RNA regulators called small RNAs (sRNA) [10], [11]. Previous investigations in Salmonella enterica and Escherichia coli have shown several sRNAs to be regulated by their respective σE orthologues, making it essential to be able to detect these transcripts in any complete exploration of σ-factor regulation [12], [13]. RNA-seq using next generation sequencing has the potential to overcome all mentioned limitations by offering an approach to identify all RNA species expressed within a cell without pre-selection [14], [15]. Therefore, in order to unravel the complete σE regulon in N. meningitidis, we sequenced total RNA transcriptional content of wild type meningococci and compared it with that of mseR knock-out cells in which σE is highly expressed.

Materials and Methods

Bacterial strains and growth conditions

N. meningitidis strain H44/76, B:P1.7,16:F3-3: ST-32 (cc32), is closely related to serogroup B strain MC58, belonging to the same clonal complex [16]–[20]. N. meningitidis H44/76 in which mseR was replaced with the erythromycin resistance cassette ermC (H44/76ΔmseR) was described previously [3]. Meningococci were grown on GC plates (Difco), supplemented with 1% (vol/vol) Vitox (Oxoid) at 37°C in a humidified atmosphere of 5% CO2. Broth cultures were incubated in GC medium supplemented with 1% (vol/vol) Vitox on a gyratory shaker (180 rpm) at 37°C. Where appropriate, plates were supplemented with erythromycin (5 µg/mL). Growth was monitored by measuring optical density of cultures at 530 nm (OD530, Pharmacia Biotech Ultraspec 2000) at regular intervals. For RNA isolation, used in both RNA-seq and qRT-PCR experiments, cells were grown to exponential growth phase (OD530 = 0.5–0.6).

Total RNA preparation

Immediately after sample removal, RNA was stabilized by adding 1/10th volume of stop solution (95% ethanol/5% phenol pH 4.3) and samples were further processed after rapid cooling on ice [21]. Isolation of total RNA was performed using an acidic hot phenol method adapted from [22]. RNA molecules smaller than 200 nucleotides (nts), including 5S and tRNAs, were retained using classic acidic phenol/chloroform extractions for sample clean-up [23]. The majority of 16S and 23S rRNA molecules were selectively removed using the MICROBExpress kit (Life Technologies, 2008). RNA quality was assessed by resolving samples on 1% agarose gels and measuring RNA Integrity Number using a 2100 Bioanalyzer (Agilent Technologies) [24]. For RNA-seq experiments RNA with a RIN value ≥9.0 was used. Quantification was assessed by UV spectroscopy (Nanodrop 1000, Thermoscientific).

SOLiD sequencing

Transcriptome sequencing was performed according to the manufacturer's protocol using 1 µg RNA isolated from H44/76 wild type (H44/76 wt) or H44/76ΔmseR cells (SOLiD Whole Transcriptome Analysis Kit, Life Technologies, 2009 Rev. F). Adaptor mix A was used for RNA-adaptor ligation, ensuring strand specificity. In short, this is achieved by ligating the fragmented RNA molecules using double stranded adapters with single stranded random hexamers overhangs either on the 5′ or the 3′ end. Subsequent sequencing is performed using primers that anneal to only one of these two adapters. Bead preparation was performed according to the SOLiD v3 Plus system's instructions (Life Technologies, 2009). The capacity of a single full slide was used, using barcodes to individually tag two H44/76 wt and two H44/76ΔmseR technical replicates, consisting of separately constructed cDNA libraries.

qRT-PCR

After DNase treatment (Roche, 2008) 2.5 µg RNA was used as input for reverse transcription, performed as per ThermoScript kit instructions (Invitrogen, 2010) with a single adaptation: all procedures were performed in duplo, substituting RT with nuclease-free water in 1 reaction to serve as no-RT controls. The resulting RT+ and RT- cDNA samples were serially diluted and PCR amplification of the constitutively expressed outer membrane protein rmpM was used to titrate relative total cDNA concentrations of the H44/76 wt and H44/76ΔmseR samples. The cDNA samples were adjusted accordingly for use in subsequent real-time PCRs (Lightcycler, Roche) (see Table S1 for primers used). Crossing point (Cp) values were calculated manually for all genes using the Fit Points Method in the LightCycler Software package (Roche, version 4.05) [25], [26]. All experiments were carried out in triplo. To normalize target gene Cp values, mean Cp values of the internal control rmpM of both H44/76 wt and H44/76ΔmseR were deducted from the Cp value of the respective target genes. Unpaired two-sample t-tests were performed for assessing the statistical significance of the difference between the means of these normalized Cp values. This was followed by a one-tailed test to derive p-values with 0.05 chosen as the threshold for statistical significance.

Data analysis

The SOLiD Accuracy Enhancer Tool (SAET) was used to correct miscalled bases and subsequently the P2 sequencing adapter was trimmed using a tool provided by Nicholas Socci (www.pyrodigm.com). The corrected and trimmed reads were mapped to the recently published H44/76 Whole Genome Shotgun (WGS) sequence (GenBank accession AEQZ00000000) using BWA [19], [27]. This mapping was performed in two steps. First the reads were mapped with default BWA settings (version 0.5.7 r1310). Subsequently the unmapped but high quality reads are mapped in a second step with less stringent settings allowing for 2 mismatches in the first 25 nts and 4 mismatches for the total read (up to 50 bp). These mismatches were allowed to compensate for limited genome sequence errors, i.e., the reference H44/76 WGS genome sequence was obtained by 454 pyrosequencing prone to inaccurate sequencing of long homopolymer stretches, and to allow possible mapping of post-transcriptionally altered RNA molecules [19]. The results of these mappings were merged for further analysis using SamTools [28]. For every ORF and intergenic region (IGR) we counted the number of aligned reads using the Rsamtools package [29], [30]. Reads that overlapped with two ORFs and/or IGRs both contributed to the read count of these two ORFs/IGRs. Reads that could not be unambiguously mapped to a single location were ignored in subsequent analysis. This concerned an estimated 2% of the genome involving repetitive sequences and duplicated genes like transposase elements (73 annotated), phage related proteins (13 annotated), iron regulated proteins FrpA/B/C/D (includes several pseudogenes), the iron uptake system TpbA/B (duplicated), maf type adhesins (duplicated), Type IV pili, and rRNA genes (one copy annotated) as H44/76 does not contain large duplicated regions [19]. Library normalization and differential expression analysis for all IGRs larger than 29 nts and all ORFs on both strands was performed using DEseq [31]. For ORFs this resulted in the separate analysis of both sense and anti-sense transcription. DEseq settings were optimized for the most conservative analysis of technical replicates as described by Anders et al. [31]. The library sizes of the different samples were adjusted based on the total number of reads per library, without the use of size factor estimation. The p-values resulting from the pairwise comparison between H44/76 wt and H44/76ΔmseR were adjusted by the Benjamini-Hochberg procedure to control the false discovery rate (FDR) [32]. Visual inspection of mapped reads was facilitated with the Integrative Genomics Viewer (16).

In silico prediction of transcriptional start sites (TSSs) was performed using the Neural Network Promoter Prediction program (prokaryote organism setting), ρ-independent terminators were identified manually using Mfold or adapted from those predicted by TransTerm HP in N. meningitidis MC58 [33]–[35]. All reported nt positions are according to the published H44/76 whole genome shotgun sequence [19]. Potential mRNA targets (in the reference strain MC58) of sRNAs were predicted in silico using TargetRNA, its default settings adjusted only to allow for G:U pairs in seed [36].

Results

Output SOLiD

The total number of reads varied from 50×106 to 92×106 (Table 1). The proportion of mapped reads ranged from 50% to 68%, with the second less stringent mapping accounting for 4–6% of all mapped reads. Thirty percent to 47% of the reads mapped to tRNA or rRNA genes. In one technical replicate, H44/76ΔmseR replicate A, the proportion of reads mapped to tRNA/rRNA genes was lower compared to the other 3, while the proportion of unmapped reads was higher. This might be attributable to slightly different gel size fractioning; however this has no observable effect on the total reads mapped to the rest of the genome. Raw read mapping analysis (before library normalization) showed at least 1 read could be mapped to 2593/2598 (100%) ORFs and 1635/2599 (63%) IGRs, with 43% of the positive and 48% of the negative strand covered by 10 or more reads.

Table 1

Mapping statistics of SOLiD output.

	Number of reads ×106 in cells of
	H44/76 wt		H44/76ΔmseR
Reads	replicate A (%)	replicate B (%)	replicate A (%)	replicate B (%)
Mapped to genome
Excluding tRNA/rRNA	20 (24%)	12 (21%)	18 (20%)	12 (24%)
tRNA/rRNA genes only	37 (44%)	27 (47%)	28 (30%)	22 (44%)
Unmapped	27 (32%)	19 (33%)	46 (50%)	16 (32%)
Total	84 (100%)	58 (100%)	92 (100%)	50 (100%)

Two technical replicates were performed for both H44/76 wt and H44/76ΔmseR, resulting in a total of 280×106 reads.

Differentially expressed genes

Eleven coding genes and one non-coding gene were found to be differentially expressed between H44/76 wt and H44/76ΔmseR cells with a corrected p-value≤0.05 (Table 2).

Table 2

Differentially expressed loci, ranked by locus number, between N. meningitidis H44/76 wt and H44/76ΔmseR.

		reads/locusα
Locusβ	Geneγ	wt	ΔmseR	Fold change	p-valueδ	Descriptionε
NMH_0518		1998	5686	2.8	0.010	pepSY-associated TM helix family protein
NMH_1358		3740	1403	0.4	0.021	outer membrane OpcA family protein
NMH_2149		5216	187292	35.9	<0.0001	doxX family protein, σE operon
NMH_02153		5042	60516	12.0	<0.0001	Hypothetical protein in σE operon
NMH_2154		4807	63164	13.1	<0.0001	Hypothetical protein in σE operon
NMH_2155		3071	40218	13.1	<0.0001	Hypothetical protein in σE operon
NMH_2156	σE	1587	9997	6.3	<0.0001	RNA polymerase σE, σE operon
NMH_2157	mseR	331	-ζ	n/a	n/a	Meningococcal σE Regulator, σE operon
NMH_2253	nirK	638	117	0,2	<0.0001	aniA (MC58)
NMH_2270	opA50	1098	418	0,4	0.027	nspA (MC58)
NMH_2477	msrA	1726	138968	80.5	<0.0001	methionine sulfoxide reductase msrA/msrB
σE sRNA		7	114	16.3	<0.0001	σEsRNA

Reads/locus after library normalization by DEseq.

Locus number as annotated in H44/76 WGS (AEQZ00000000).

Gene name as annotated in H44/76 WGS (AEQZ00000000).

All loci with Benjamini-Hochberg adjusted p-value≤0.05.

Where available N. meningitidis strain MC58 orthologues are noted.

mseR artificially deleted.

In contrast to H44/76 wt cells, in H44/76ΔmseR cells no reads were detected that mapped to the mseR locus of the H44/76 reference sequence (Table 2) [19]. In addition, RNA from H44/76 wt cells did not contain reads that mapped to ermC, while H44/76ΔmseR cells in which mseR was replaced by ermC yielded reads that mapped to this locus (data not shown). Therefore, differential transcription of artificially introduced and removed genes was confidently identified. Previously, using RT-PCR and proteomics, we identified the σE operon and msrA/msrB to be upregulated in H44/76ΔmseR cells [3]. In the current RNA-seq analysis the expression of msrA/msrB (Figure 1) and the genes comprising the σE operon (Figure 2) was 81-fold and between 6-fold and 36-fold increased, respectively (Table 2). Of interest, four differentially expressed genes were identified which were hitherto not associated with σE. Of these, one was upregulated: NMH_0518 (analogous to NMB1721; encoding a hypothetical protein of the PepSY-associated TM helix family of proteins). Three genes, NMH_2270 (analogous to NMB0663; nspA), nirK (analogous to NMB1623; aniA), and NMH_1358 (analogous to NMB1053; opcA outer membrane protein), appeared to be downregulated (Table 2). However, the promoter of the latter gene contains a homopolymeric cytidine tract located at the -35 region which can be extended or shortened during DNA replication by a mechanism called slipped strand DNA replication [37], [38]. As a result, opcA displays variable expression [37]. Indeed, classic Sanger re-sequencing of this location revealed the promoter of H44/76 wt opcA contained a tract of 14 cytosines, while that of H44/76ΔmseR contained 15, the latter causing severe reduced transcriptional efficiency [37]. Therefore, we considered differential expression of NMH_1358 most likely to be the result of unrelated phase-variation independent of σE expression. A complete list of differentially expressed loci with uncorrected p-values≤0.01 can be found in Table S2. Among the 30 loci with differential sense expression, four clusters of genes and their associated IGRs, either physically arranged in an operon or with a functional relationship, were recognized; cluster 1 comprising three tRNA-Leu genes (NMH_1084, NMH_1085, and NMH_2191); cluster 2 consisting of genes encoding a restriction modification system (NMH_1821-NMH_1823), cluster 3 encoding a leukotoxin type I secretion system (NMH_1851, NMH_0534-NMH_0535) and cluster 4 encoding genes involved in a thiamin biosynthesis operon (NMH_2073, NMH_2074, NMH_2076, thiG). In these clusters the difference in expression between H44/76 wt and H44/76ΔmseR cells varies between 2-fold and 3.4-fold. Statistical analysis of the total expression levels of these clusters (i.e., using the total read count of all genes in a cluster) did not result in additional significant differential expression considering the corrected p-value. Visual inspection of pairs of 100 randomly picked non-differentially expressed genes using the Integrative Genomics Viewer software revealed consistently highly similar patterns of read distribution between H44/76 wt and H44/76ΔmseR normalized libraries, indicating that variation between strains is limited.

Figure 1

Read coverage visualization of the upregulation of msrA/msrB in H44/76ΔmseR (blue) versus H44/76 wt (red).

Transcription on the + strand is visualized on the positive x-axis, - strand transcription on the negative x-axis. Nucleotide position refers to contig9 of H44/76 WGS (AEQZ01000046.1). The relatively high peak of the last 50 nucleotides of msrA/msrB is present in both H44/76 wt and H44/76ΔmseR in both replicates. Expression of msrA/msrB through its ρ-independent terminator can be seen (see text for details), explaining the upregulation of anti-sense transcription of NMH_2475 and NMH_2476 shown in table 3.

Figure 2

Read coverage visualization of the upregulation of the σE operon in H44/76ΔmseR (blue) versus H44/76 wt (red).

Transcription on the + strand is visualized on the positive x-axis, - strand transcription on the negative x-axis. Nucleotide position refers to contig6 of H44/76 WGS (AEQZ01000043.2). The largest absolute transcriptional differences can be found in the first 2 genes of the σE operon, NMH_2149 and NMH_02153.

Table 3

Differentially expressed loci in anti-sense direction, ranked by locus number, between N. meningitidis H44/76 wild type and N. meningitidis H44/76ΔmseR.

		reads/locusα
Locusβ	Geneγ	wt	ΔmseR	Fold change	p-valueδ	Descriptionε
NMH_0763	nqrE	30	121	4.0	0.020	NADH: ubiquinone oxidoreductase
NMH_2154		157	19	0.1	<0.001	Hypothetical protein in σE operon
NMH_2156	σE	26	141	5.4	0.001	RNA polymerase σE, σE operon
NMH_2475		17	461	27.1	<0.001	mechanosensitive ion channel family protein
NMH_2476		6	725	120.8	<0.0001	Competence/damage-inducible CinA C-terminal
NMH_2477	msrA/msrB	771	118	0.2	<0.0001	methionine sulfoxide reductase msrA/msrB

Reads/locus after library normalization by DEseq.

Locus number as annotated in H44/76 WGS (AEQZ00000000).

Gene name when annotated in H44/76 WGS (AEQZ00000000).

All loci with Benjamini-Hochberg adjusted p-value≤0.05.

Where available N. meningitidis strain MC58 orthologues are noted.

Differential transcription in intergenic regions

IGRs may contain non-coding RNAs involved in the riboregulation of target mRNAs [10], [11]. Of all IGRs of H44/76, only the IGR between NMH_1566 and NMH_1568 (analogous to NMB1826 and dnaE in MC58) showed differential expression with 16 times higher expression in H44/76ΔmseR versus H44/76 wt cells (Table 1). A gene encoding a sRNA of 74 nts can be identified in this IGR (Figure 3). Its TSS deduced from read mapping in RNA-seq is at contig4-nt position 59241 (AEQZ01000037.1), which is 11 nts downstream of an in silico predicted -10 region of a putative promoter and 3 nts downstream of the predicted TSS (Figure 4). The transcript ends within the T-tract of a ρ-independent terminator, with the last reads ending at contig4 nt position 59314 (Figure 4). The gene encoding this sRNA was found to be conserved in all currently sequenced complete and WGS genomes of N. meningitidis, N. gonorrhoeae and N. lactamica strains (query coverage 100%, identity ≥94%) and in several but not all commensal Neisseria species. Therefore, putative interaction between this sRNA and the 5′ UTR of target mRNAs was assessed using TargetRNA in the related strain MC58 [36]. This way, NMB0205 (fur), NMB0792 (nadC), NMB0810 (encoding a putative TetR family transcriptional regulator) and NMB1224, NMB1914, NMB2014, and NMB2110, all encoding hypothetical proteins, were identified as putative targets of this novel sRNA.

Figure 3

Read coverage visualization of the upregulation of the σEsRNA in H44/76ΔmseR (blue) versus H44/76 wt (red).

Transcription on the + strand is visualized on the positive x-axis, - strand transcription on the negative x-axis. Nucleotide position refers to contig4 of H44/76 WGS (AEQZ01000037.1). The black flag indicates a predicted TSS, the closed circle indicates a predicted ρ-independent terminator. The ORF of hypothetical protein NMH_1566 shows no expression. dnaE (NMH_1568), truncated in this picture with only 150 of 3435 nts shown, is not differentially expressed.

Figure 4

Sequence and predicted secondary structure of the σEsRNA.

(A) Mfold prediction of the secondary structure of the σEsRNA. (B) Sequence of the σEsRNA gene including promoter. The boldfaced -35 and -10 regions and transcriptional start site are predicted in silico. The boxed +1 nucleotide is the experimentally validated transcriptional start site of the in vivo present transcripts.

Differential expression between H44/76 wt and H44/76ΔmseR of NMH_0518, NMH_2270, nirK, msrA and NMH_2149 (representing the σE operon) was also assessed by qRT-PCRs on reverse transcribed RNA. The internal control rmpM showed very similar Cp values for both cDNA samples (11.7±1.4 vs 11.9±1.5; p = 0.4573), consistent with RNA-seq results showing less than 6% difference between total normalized reads mapped to rmpM (data not shown). Four genes showed similar direction of differential expression (three significant), confirming the RNA-seq data, although fold differences show varying order of magnitude in comparison (Table S3) [39]. Detection of transcription of NMH_2270 and the σE sRNA was insufficient to yield reliable results.

Antisense transcription

Cis-antisense transcripts of genes have been implicated to interact with its mirror image mRNA molecule, regulating transcription, translation or degradation [15], [40]. We investigated differential expression by considering the number of reads per locus in the anti-sense direction of all ORFs and detected 6 genes to have differentially expressed anti-sense transcripts (Table 3). In H44/76 wt cells 771 reads were mapped to the anti-sense strand of msrA/msrB, while in H44/76ΔmseR cells this number was 6.4 fold reduced. In H44/76ΔmseR cells, antisense transcription of both NMH_2475 and NMH_2476 was enhanced 27 and 120 fold, respectively, possibly due to read-through from the strongly transcribed msrA/msrB locus (Figure 1). In addition, an 8.3 fold reduction in antisense transcription of NMH_2156 (part of the σE operon) was observed in H44/76ΔmseR cells. Lastly, antisense transcription of nqrE (NMB0565, Na+-translocating NADH-quinone oxidoreducase subunit E) and σ itself was 4.0 and 5.4 fold increased, respectively.

Discussion

Positive feedback loop of the σE operon and control of msrA/msrB expression by σE

Previous reports indicated a relatively small σE operon regulon in pathogenic Neisseria when compared to other species with a much wider range of habitats, e.g. Bacillus subtilis, where it directs the transcription of many genes that contribute to the formation of mature spores [3], [6], [41]. Here we used RNA-seq to explore this regulon in N. meningitidis in greater depth. We found 4 novel genes of which their expression was regulated by σE. We did not detect any similarities between the promoters of these genes and the promoter consensus of the σE operon and msrA/msrB described previously, implicating an indirect effect of σE [3]. An additional gene, the OpcA superfamily protein NMH_1358, also showed differential expression between H44/76 wt and H44/76ΔmseR cells, but this was explained by phase variation. Evidently, the σE operon and msrA/msrB quantitatively show the most dramatic upregulation of sense transcription. The magnitude of differences in expression of the newly identified genes is smaller, which may indicate that the major function of σE is to provide positive feedback to its own operon and to induce a strong expression of msrA/msrB. However, even relatively small transcriptional differences can have a major impact on genome gene expression, e.g. change in transcription factors may propagate among many genes.

Discovery of novel differentially expressed genes

In general, σE is found to be involved in the regulation of genes encoding proteins residing in the outer membrane or periplasmic space [2]. Of the 4 newly identified σE regulated genes, 3 encode outer membrane proteins or secreted proteins. Of these, only one, NMH_0518 (NMB1721), encoding a protein of the PepSY-associated TM helix family of proteins, was upregulated and its protein structure suggests that it is held at the cell wall or is secreted [42]. PepSY domains are present in a diverse family of secreted and cell-wall associated proteins and it has been suggested that they act as protease inhibitors, and as such regulating protease activity in the local environment to the cell, with possible significance for pathogenicity [42]. The 2 other σE regulated protein coding genes are nirK (aniA), and NMH_2270 (nspA); both appeared to be downregulated. NirK is an anaerobically induced outer membrane protein producing nitric oxide through denitrification [43]. Its gene is activated by the fumarate and nitrate reduction regulator FNR when iron concentration is low [44]. When iron concentration is not limiting than nirK is repressed by the nitrite-sensitive repressor NsrR, replacing FNR at the nirK promoter. NsrR contains a labile Fe-S cluster and when intact it binds to the nirK promoter and represses nirK [44], [45]. Methionine sulfoxide reductase, ubiquitously present across all organisms and encoded by msrA/msrB, might stabilize iron-sulfur clusters as has recently been shown in Saccharomyces cerevisiae [4], [46]. We suggest that, considering the pervasive similarities of msrA/msrB in both eukaryotes and prokaryotes, high expression of the methionine sulfoxide reductase may favorably affect the stability of the Fe-S complexes of NsrR increasing the relative NsrR-S-Fe concentration shifting the balance towards more repression of nirK.

NMH_2270, encoding the Neisserial surface protein NspA is a vaccine candidate that can bind to human factor H and enhances complement resistance [47]. nspA is activated by the Ferric Uptake Regulator protein (Fur) when iron concentration is not limiting and downregulated when iron concentrations are low [48]. In conditions with sufficient iron, Fur also downregulates itself via a negative feedback loop. The Fur protein contains Fe-S clusters and in the complexed form it binds to its own promoter and represses transcription of itself. The Fur protein contains Fe-S clusters and again with high expression of methionine sulfoxide reductase the Fe-S clusters in Fur may be stabilized [48]–[50]. Reduced transcription of fur was detected in H44/76ΔmseR cells albeit not significant, therefore a(n) alternative mechanism(s) may explain the downregulation of nspA.

As described by Anders et al., DEseq was designed to improve on previous tools like edgeR, which tends to overestimate differential expression in the lower range of expression, while underestimating differential expression when genes are relatively highly expressed [31], [51], [52]. Nevertheless, examples of the first may still be present, with σE and nqrE, listed in Table 2 as being upregulated in their antisense direction. Both have very low expression in both the H44/76 wt and H44/76ΔmseR. Upon closer visual inspection all reads are scattered within these relatively large loci without any apparent organization. They are therefore either near the detection limit of this experiment or not likely to represent any biological relevance. Examples representing the other side of the quantitative transcriptional spectrum are the genes encoding tRNA-Leu; NMH_1084, NMH_1085 and NMH_2191, which show 2.0 to 3.4 fold upregulation (Table S2). Especially, NMH_2191 has a high number of uniquely mapped reads with 224294 readsn the H44/76 wt and 771702 reads in H44/76ΔmseR cells. This extremely high amount of reads is seen only with tRNAs, 5S RNA and several mRNAs coding for ribosomal proteins. Increased expression of tRNA-Leu genes can be caused by an increased demand for leucine incorporation in leucine-rich proteins. Analysis of the σE operon using Composition Profiler reveals relatively high usage of leucine (p<0.001) compared to the average usage of this amino acid in the total protein coding content of the H44/76 genome [53]. Therefore, artificially created constitutional overexpression of this operon increases the intracellular demand for leucine dramatically, which provides a biological explanation for the observed increase of reads mapped to three tRNA-Leu loci (NMH_2191, NMH_1084 and NMH_1085).

There are three additional clusters of genes which showed differential expression without reaching significance when adjusted according to Benjamini-Hochberg to control for FDR, comprising 1) a restriction modification system, 2) the leukotoxin type I secretion system and 3) the thiamin biosynthesis operon (Table S2). The increased expression of the restriction modification locus in H44/76ΔmseR cells may be caused by phase variation. The methyltransferase (mod) protein NMH1820 (analogous to NMB1375, pseudogene part 1) contains tetrameric AGCC repeats causing a frameshift that leads to a premature stop codon [54]. After a stop codon, ribosomes dissociate from mRNAs, making that part of mRNA less protected against degradation and lowering the number of reads of the mRNA posterior of the stop codon in H44/76 wt [55]. However, Sanger re-sequencing of this repeat region in H44/76 and H44/76ΔmseR using genomic DNA extracted from the samples used for RNA-seq did not show differences in the number of tetrameric repeats. Different mechanisms for differential expression, e.g. transcriptional slippage might be involved, but the lack of genomic differences indicate in favor of the current finding of non-significant differential expression [56].

It remains unclear why the leukotoxin type I secretion system is 2-fold downregulated and genes encoding proteins involved in thiamin biosynthesis were 2.2 to 2.8 fold upregulated in H44/76ΔmseR cells. Of note, the thiamin biosynthesis pathway also contains enzymes with iron-sulfur clusters which might suggest a relation with the strong upregulation of msrA/msrB [57]. Interestingly, oxidative stress leads to auxotrophic requirements for thiamin in Salmonella enterica and Erwinia chrysanthemi. [58], [59]. In E. chrysanthemi, thiamin auxothrophy can also results from depletion of the Isc system, which is known to assist Fe-S cluster biogenesis, resulting in drastic alterations in virulence [59]. In meningococci, increased expression of MsrA/MsrB can prevent Fe-S cluster disintegration, thus supporting increased thiamin biosynthesis.

Differential expression of non-coding RNA transcripts

Only one IGR showed differential expression between H44/76 wt and H44/76ΔmseR cells. The expression of a sRNA transcript located in this IGR between NMH_1566 and NMH_1568 is upregulated in the H44/76ΔmseR mutant (Figure 3). This σEsRNA contains no open reading frame and a BLAST search shows no similarity on the nucleotide level with species other than Neisseria. Putative mRNA targets of this σEsRNA are, in addition to several genes coding for hypothetical proteins, Fur, a NadC transporter family protein and a TetR family transcriptional regulator. By regulating the expression of fur it may also effect the transcription of nspA. Moreover, regulation of fur and the TetR family transcriptional regulator would extend the impact of σE on the global protein expression pattern. Interaction of many sRNAs and their mRNA targets has been shown to be enhanced by the bacterial sm-like protein Hfq by facilitating base pairing and stabilizing the double stranded RNA complex [60]–[62]. Two proteomic analyses of the influence of deletion of hfq on the protein expression patterns of N. meningitidis strains MC58 and H44/76 have shown an attenuated phenotype with an extensive profile of differentially expressed proteins [63], [64]. None of the putative mRNA targets of this σEsRNA have been shown to be differentially expressed in these studies. This might imply that the targets identified in silico are not relevant in vivo or the interaction between σEsRNA and its targets is Hfq independent. Since the sRNA expression is increased in H44/76ΔmseR cells one might have anticipated differential expression of its target genes. Nevertheless, none of the putative target genes of the sRNA show differential expression between H44/76 wt and H44/76ΔmseR cells. When sRNAs interact with target mRNAs, the latter are often destabilized due to the activity of the RNase E and RNAse III dependent degradosome [55], [65]. As a result, a target mRNA is degraded to smaller fragments which can be observed as a less intense or fragmented signal in northern blots or with a reduced or absent signal in RT-PCRs [62], [66]. RNase III is also used to process total RNA extracted from the bacterial cultures during SOLiD library preparation. In addition, with the SOLiD sequencing protocol employed here, reads between 25 and 50 nts were generated from all RNA molecules with a native or processed size of 50 nts or larger. Together this makes RNA-seq by SOLiD sequencing with the currently used protocol less suitable for detection of sRNA-mRNA interactions. Further study is underway to assess the function of this interesting non-coding sRNA.

ρ-independent terminator read-through

NMH_2475 and NMH_2476 show increased antisense transcription in H44/76ΔmseR cells. This can be explained by read-through past the ρ-independent terminator of the highly expressed msrA/msrB locus located diametrically oriented just downstream of these genes. This is visualized by small peaks of the blue H44/76ΔmseR line on the negative x-axis (representing negative strand expression) at the NMH_2475 and NMH_2476 loci located on the positive strand, while the red H44/76 wt line shows no antisense expression (Figure 1). DEseq analysis shows that sense expression of NMH_2475 and NMH_2476 is not significantly altered; therefore, the potential stereometrical obstruction of opposing RNA polymerases does not seem to influence the expression of these genes. Similar inefficient termination of transcription has been observed in 454 pyrosequencing RNA-seq experiments in Helicobacter pylori, confirmed by northern blot probe hybridization with larger RNA molecules than expected [67].

This report provides an in-depth analysis of the concise but complete σE factor regulon in the Beta-proteobacterium Neisseria meningitidis. This relatively small regulon implies a strictly defined response suitable for the relatively stable niche, the human throat, where N. meningitidis resides.

List of primers used in qRT-PCR experiments.

(XLS)

Click here for additional data file.

List of differentially expressed loci with uncorrected -values≤0.01.

(XLS)

Click here for additional data file.

Results of qRT-PCR experiments.

(XLS)

Click here for additional data file.