Transcriptome profile of a bovine respiratory disease pathogen: Mannheimia haemolytica PHL213

RNA isolation

M. haemolytica PHL213 was cultured in brain heart infusion (BHI) to mid-log phase (OD₆₂₀ = 0.8). Cells from a single culture were treated with RNAprotect reagent (Qiagen, Valencia, CA) and stored at -80°C for subsequent RNA isolation. Total RNA from this single culture was extracted using the RNeasy mini kit (Qiagen, Valencia, CA), following manufacturer's protocols. It is to be noted that this kit allows for the extraction of transcripts that are at least 200 nucleotides and larger. RNA preparations were treated with RNase-free DNAse (Invitrogen, Carlsbad, CA) and the integrity of the RNA was determined using Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA). RNA sample with RNA Integrity Number (RIN) of 8 was used for the RNA-Seq experiment. From total RNA, mRNA was enriched by removing the rRNAs using MICROBExpress™ kit (Ambion, Austin, TX). This enrichment step specifically removes large rRNAs; small RNAs (i.e., tRNA and 5S rRNA) are not removed. In the first step of the MICROBExpress™ kit procedure, total RNA was mixed with an optimized set of capture oligonucleotides that bind to the bacterial 16S and 23S rRNAs. Next, the rRNA hybrids were removed from the solution using derivatized magnetic microbeads. The mRNA remained in the supernatant and was recovered by ethanol precipitation and quantified by Bioanalyzer 2100. Our RNA preparation did not include entities < 200 nucleotides in length.

RNA-Seq

A cDNA library was constructed using the Illumina mRNA-Seq sample prep kit (Illumina, San Diego, CA) with 100 ng RNA enriched for mRNA isolated from a single in vitro culture, following manufacturer's instructions. mRNA was chemically fragmented and randomly primed for reverse transcription and second-strand synthesis. The resulting cDNA was end-repaired and 'A' overhangs were added. Illumina paired-end sequence adaptors were ligated to the cDNA fragments. Fragments of approximately 200 bp were isolated from a 2% agarose gel and amplified (18 cycles) according to the Illumina protocol. Bioanalyzer 2100 (Agilent) was used to quantify and confirm the fragment size of each library. 1 nM of mRNA-seq library sample prepared for sequencing on the Illumina GAII (San Diego, CA) was denatured and diluted to 6 pM for clustering (v2) according to the manufacturer's protocol. Single read sequencing of the clustered flow cell was performed using Illumina's SBS chemistry (v3) and SCS data analysis pipeline v2.4. Flow-cell image analysis and cluster intensity calculations were carried out by Illumina Real Time Analysis (RTA v1.4.15.0) software. Subsequent base-calling was performed using the Illumina GA Pipeline v1.5.1 software. The resulting Illumina reads were quality-filtered by removing reads containing Ns.

Alignment

The sequencing experiment produced 9,055,826 reads. FASTQ reads generated by Illumina were converted to Sanger FASTQ format using Perl scripts from the Mapping and Assembly with Qualities (MAQ) software package [16]. Reads (Sanger FASTQ format) were mapped to the 2.6 Mb M. haemolytica PHL213 [GenBank: AASA00000000] reference genome using Bowtie [17]. The parameters in Bowtie that control the speed and sensitivity were adjusted as follows: reads with no more than 2 mismatches per read (n = 2) were aligned, and any reads mapped to more than one location across the genome (ambiguous reads) were discarded (m = 1). Post alignment, a human-readable sequence alignment/map (.SAM) format file was converted to a "pileup" format file using SAMTools [18]. This pileup file contains the count of reads per base aligned to each location across the length of the genome. The SAM file was also converted into a binary alignment/map (.BAM) format. These BAM formatted files are necessary for visualization of read alignments in Artemis viewer. The Artemis browser enabled the visual/manual inspection of alignment results in the context of the existing genome annotation. The pileup file, in conjunction with the annotation information of M. haemolytica PHL213, was processed using in-house Perl [19] scripts. Data generated from the RNA-Seq experiment was submitted to the NCBI Sequence Read Archive [SRA049621.1] as reads in FASTQ format [SRR402063.1] and the .BAM alignment file [SRR402079.4] generated by aligning the reads to the reference genome.

Analysis of expressed intergenic regions

Identifying expressed regions within the genome that have not been previously annotated will improve the existing structural annotation of the M. haemolytica PHL213. Prior to the analysis of expressed regions in the genome, we determined the signal to noise ratio cutoff for background expression using the pileup file. Coverage depth (reads per base) greater than the lower tenth percentile of all reads was considered to be expressed and in this dataset, this corresponded to 7 reads/base [8, 20]. Based on this read/base cutoff, expressed intergenic regions (EIRs) were identified by applying an additional length cutoff of 70 bp. Shorter regions (less than 70 bp) were discarded to reduce the number of false positives. Custom Perl scripts were written to parse the pileup file and the existing genome structural annotation to identify (i) expressed annotated regions, (ii) expressed regions previously not annotated and, (iii) regions that are annotated but are not expressed. All EIRs were further analyzed using BLASTX [21] searches to determine their protein coding potential. If an EIR was found to be a perfect match (~100% coverage) for a protein, it was classified as a putative novel protein coding region. All EIRs with partial BLASTX hits were evaluated for the presence of an alternate start site or mutation in the start or stop codon associated with the annotated region. If the BLASTX search revealed a frameshift mutation, the EIR and the gene associated with the frameshift mutation were classified as a frameshift. EIRs with poor BLASTX hits and without any association to genes containing annotation errors were excluded from further analysis. EIRs without BLASTX hits were considered to be potential small non-coding RNA.

The Prokaryotic Promoter Prediction (PPP) program [22] (from PePPER suite [23]) and Transterm HP [24] were used to predict promoters and rho-independent terminators, respectively, in the forward and reverse strands of the M. haemolytica PHL213 genome. The locations of promoters and terminators were organized into .GFF files. A Perl script was written to identify putative sRNA i.e. EIRs with promoters or terminators associated to their loci. EIRs with no computationally-predicted promoters or rho-independent terminators were searched against the Rfam database [25] to determine whether these sequences were annotated in Rfam. EIRs that could not be classified as sRNA by Rfam were excluded from further analysis.

Analysis of expressed annotated regions

Using the annotation information (gene loci) of M. haemolytica PHL213 and the pileup file, all annotated regions that were expressed above the background signal to noise ratio cutoff with at least 60% coverage were considered to be expressed, which accounts for uniform evaluation of varying gene lengths. Similar measures have been used in other transcriptome profiling studies [8, 26, 27]. Annotated regions below 60% coverage were considered as 'not expressed' under the current experimental conditions. After having identified expressed genes, operon structures within the genome were also defined. The first step towards identifying an operon was to identify co-expressed pairs of coding regions. Two regions were considered to be co-expressed when they were identified as expressed on the same strand (5' to 3' or 3' to 5') and the region between them was also expressed. After such co-expressed pairs were identified, they were extended to construct operons by including additional co-expressed pairs in the vicinity satisfying the same conditions for co-expression as described earlier. Operon structures identified by RNA-Seq were compared to the computationally-predicted operon structures described by the Database for prOkaryotic OpeRons (DOOR) [28] for cross validation.

Read alignment to the M. haemolytica PHL213 genome

The M. haemolytica PHL213 is a 2.6 Mb draft genome containing 2,837 annotated regions of which 2,695 are protein coding with a 40% G+C content [29]. For structural annotation of M. haemolytica at the RNA level, the transcriptome of M. haemolytica PHL213 was sequenced using RNA-Seq. Sequencing-based analysis of the transcriptome overcomes the limitations of the hybridization-based microarray approach. Head-on comparison of RNA-Seq with microarrays has shown that RNA-Seq has negligible technical variability [30], making it possible to obtain a reliable estimate of gene expression without replicate analysis. Therefore, we applied RNA-Seq for re-annotation of M. haemolytica and conducted the analysis from a single in vitro experiment. Reads with an average length of 76 bp generated on the Illumina platform were mapped onto the reference genome using the Bowtie read alignment program. Bowtie is an ultrafast, memory efficient alignment program that uses the Burrows-Wheeler transform [31] with a novel quality backtracking algorithm that permits mismatches. Bowtie performs better than Short Oligonucleotide Analysis Package (SOAP) [32] and MAQ, and its sensitivity at aligning reads is as good as both SOAP and MAQ. Of the 9,055,826 reads generated by Illumina, 3,917,458 reads (43.26%) that mapped uniquely to the genome were used for downstream analysis. 2,989,603 reads (33.01%) failed to align due to mismatches. The remaining 2,148,765 reads (23.73%), which mapped to more than one location in the genome (ambiguous reads), were excluded from analysis. For annotation purposes, reads that map to unique locations alone are used [8, 33–37]. The cutoff value for true-positive expression of a coding region of 7 reads/base was calculated from the expression (number of reads per base) in the tenth percentile of all reads [8, 20], as we did earlier for RNA-Seq based re-annotation of another BRD pathogen Histophilus somni [8].

Expressed intergenic regions

We used the existing annotation of open reading frames in M.haemolytica PHL213 i.e., locus of each gene in the genome and reads identified as expressed by RNA-Seq, to identify expressed intergenic regions (EIRs). We identified 630 EIRs, previously un-annotated as expressed, of a minimum length of 70 bp. Each EIR was further characterized by adding computationally-predicted promoter and rho-independent terminators. Prokaryotic Promoter Prediction (PPP) identified 11,847 promoter regions and Transterm HP identified 1,204 rho-independent terminator regions, in forward and reverse strands of the genome. Identified EIRs, in conjunction with existing gene annotation information and loci of regulatory signals, were subjected to the analysis workflow described in Figure 1.

Artemis is a genome browser and annotation tool that allows visualization of sequence features, next generation sequencing data, and the results of the analyses within the context of the genome sequence [38]. The Artemis genome browser illustrates all the six reading frames of the genome sequence along with the translated amino acid sequences, start and stop codons, as well open readings frames across the length of genome (Figure 2). We visualized the alignment file generated by Bowtie, the gene annotation file, promoter and terminator loci, and EIRs in Artemis. Artemis generated a base coverage graph, giving a pictorial representation of the expression in various regions of the genome.

Novel protein coding regions

The protein coding potential of EIRs was determined by conducting BLASTX searches with the translated nucleotide sequence of EIRs, against the protein database containing all bacterial species. BLASTX results showed that 14 EIRs had full length matches to target sequences, indicating their potential for coding proteins. The Artemis browser was used to identify the boundaries of these 14 potential novel protein coding regions (Figure 3). These novel protein coding regions had an average G+C content of approximately 46%. The length of these regions was between 37 to 200 amino acids. While the RNA-Seq experiment itself was not strand specific, strand specificity of novel protein coding regions was inferred from the proteins identified as ~100% matches to these EIRs in BLASTX. EIR MHP4 aligned to PG1 protein of Lactobacillus crispatus ST1 while MHP12 aligned to serine acetyltransferase of Haemophilus influenzae NT127. The rest of the EIRs (Table 1) aligned to proteins classified as hypothetical.

Corrections made to the existing genome annotation

Artemis genome browser creates open reading frames (ORFs) of a desired minimum length. It identifies ORFs as regions between two consecutive stop codons with the specified minimum length. Thus ORFs corresponding to EIRs can be generated and visualized in this browser. RNA-Seq based expression in relation to the existing genome annotation, when visualized in Artemis, enabled the identification of the actual locus and length for some of the annotated proteins. We identified 4 genes with a mutated start codon. This anomaly could be the result of computational gene prediction programs identifying the next available "AUG" as the start codon (Figure 4). Our observation is substantiated by the consecutive expression of an identified EIR preceding the 5' region of these genes. 4 genes had a mutation that led to the replacement of start codon by a leucine (L).

Since Artemis allows marking the start codons within an ORF, it is possible to identify alternate start sites, if any, for any ORF associated with an EIR. ORFs created from EIRs in Artemis revealed possible alternate starts sites for 7 genes. BLASTX searches of the EIR and its translated protein revealed that the actual start site varied with respect to previous annotation (Figure 5). Where there was a discrepancy between the existing annotation and the current transcriptome based identification of start site, the consensus of start site of similar proteins identified in BLASTX was used to determine the actual start site. The suggested revisions to existing annotation are in Table 2 (detailed results in Additional file 1). BLASTX searches of EIRs also revealed mutations that lead to disruption of a protein coding region, resulting in a frameshift (Figure 6). Two such EIRs which had BLASTX alignments revealed frameshifts which would otherwise be protein coding regions (Additional file 1).

Small RNA

Small RNA are known to have regulatory roles in Escherichia coli, Staphylococcus aureus, Pseudomonas aeruginosa, Vibrio cholera and many other bacterial pathogens [39]. Genome-scale identification of sRNA using RNA-Seq is reported for E. coli [40] and Vibrio cholerae [41], among other pathogenic bacteria. The identification of the loci of sRNA in the genome is an important pre-requisite for understanding their role in modulating bacterial physiology and virulence [42]. sRNA are synthesized by RNA polymerase (RNAP) in a manner analogous to the synthesis of any RNA in bacteria (mRNA, rRNA, tRNA); sRNA promoters could be regulated by transcription factors or use of alternative sigma factors [43]. Therefore, the presence of promoters and terminators for potential sRNA [8, 44] identified by experimental approaches like RNA-Seq, increases the confidence in their identification. Although the RNA extraction protocol used in this study does not facilitate extraction of smaller transcripts and strand specificity is lost during cDNA synthesis, we identified potential sRNA. EIRs with no protein coding potential, as observed via BLASTX searches, were considered to be candidate sRNA. It is possible that EIRs with no BLASTX matches are non-conserved ORFs; since there are no in silico methods to validate this assumption, we chose to consider all EIRs with no BLASTX as candidates for small RNA analysis. Candidate sRNA loci were searched for the presence of a promoter or terminator. For 44 EIRs that had no BLASTX matches, a promoter or a rho-independent terminator was identified either on the forward or the reverse strand (Figure 7) of their locus. Promoters/terminators were present in the transcriptional regulatory regions, i.e., a promoter was present in the -1 to -35 region or a terminator was present in the +1 to +20 position at the end of the EIR. Therefore, we classified the 44 candidate sRNA as potential novel small RNA in the M. haemolytica PHL213 genome (Table 3). The average length of the identified novel sRNA was approximately 100 bp and ranged between 70 to 253 bp. The average G+C content of sRNA was 34.35%, which is relatively lower than the G+C content of the genome. All identified sRNA had a promoter associated with their locus and sRNA MHS17 also had an associated terminator. When sequences of the identified sRNA were searched in the Rfam [25] database to identify their function, no matches were found.

Gene expression and operons

The M. haemolytica PHL213 genome consists of 2,837 annotated genes, 2,695 of which code for proteins. Genes were considered to be expressed if 60% of the gene length had at least 7 reads aligned/nucleotide. Based on this criteria, 2,506 of all annotated regions in the genome (87.63%) were identified as expressed with 95.25% coverage i.e. approximately 95% of the sequence of the annotated region had at least 7 reads aligned/nucleotide. Expressed annotated genes and their coverage are documented in Additional file 2.

Functional analysis of the expressed annotated regions was based on the existing annotation of M. haemolytica genome available at NCBI. It is interesting to note that genes that are described as virulence factors are also expressed under normal culture conditions. For example, genes related to leukotoxin (MHA_0253, MHA_0254, MHA_0255 and MHA_0266), an important virulence factor, were all expressed. Also, the 12 capsule forming genes whose role in virulence includes adherence to host and resistance to serum-mediated killing and phagocytosis [45, 46] were all found to be expressed. In addition to these, we also found that 40 genes associated with lipopolysaccharide or lipoproteins and contribute to virulence by initiating an inflammatory cytokine response [45, 47] to be expressed. Genes responsible for forming the type IV pilus associated with M. haemolytica that is responsible for DNA uptake, adhesion, and motility [48] were expressed. Filamentous hemagglutinin genes of M. haemolytica (MHA_0866, MHA_0867), responsible for adhesion to host mucosa [49], were expressed. Adhesins play an important role in virulence, and all annotated genes related to this function, such as MHA_2262, MHA_0708, MHA_2492, MHA_2701, MHA_1367, MHA_0563 and MHA_2800, among others, were all identified as expressed in our experiment. Genes responsible for resistance towards antibiotics such as β-lactams, tetracycline, streptomycin, and sulfonamides [45] in M. haemolytica were also expressed. Annotated regions that were not expressed had coverage of only 30%. Of the 331 annotated regions that were not expressed 236 were annotated as "hypothetical proteins" and 26 were "hypothetical bacteriophage proteins."

Using expression patterns of coding regions, we identified paired gene expression and operon structures. RNA-Seq based operon structures were compared to the computationally predicted structures using DOOR [28]. We identified 1,086 co-expressed pairs of genes that could be organized into 518 potential operons. DOOR predicted 1,295 co-expressed pairs forming 599 operons (Additional file 3). The overlap between RNA-Seq based and DOOR-based co-expressed pairs was 854. Our study identified relatively fewer co-expressed pairs as compared to DOOR. This could be due to the fact that 331 of the 2,837 annotated regions were not expressed in our dataset. Furthermore, this method cannot detect genes whose expression is suppressed by polar mutations. The single nucleotide resolution map enabled the identification of co-expressed pairs and definition of operon structures and regulatory patterns. Availability of operon structures will facilitate understanding the coordinated regulation of genes in M. haemolytica to moderate metabolic pathways under different environmental conditions.