Whole-transcriptome, high-throughput RNA sequence analysis of the bovine macrophage response to Mycobacterium bovis infection in vitro

Summary statistics for the RNA-seq data

All 14 RNA-seq libraries were sequenced across seven lanes of one Illumina® flow cell with a mean of 25.5 million reads (range: 20.7 million to 29.3 million reads) generated per lane. Deconvolution and filtering of sequence reads to remove adapter-dimer sequences yielded a mean of 11.3 million reads (range: 2.7 million to 18.0 million reads) per individual RNA-seq library. Subsequent alignment of the filtered RNA-seq reads to the B. taurus reference genome (Btau 4.0.63 genome release) yielded a mean of 7.2 million reads (63.6%) for each RNA-seq library that mapped to unique locations in the bovine genome; a mean of 3.3 million reads (29.3%) for each library that mapped to multiple locations in the genome; and a mean of 0.8 million reads for each library (7.1%) that did not map to any genome location (Figure 1A). The number of reads per individual RNA-seq library is provided in Additional file 1: Table S1.

Further analysis of the individual library reads mapping to unique locations in the B. taurus reference genome (7.2 million reads) revealed that a mean of 5.3 million reads (73.6%) aligned to exonic regions. 1.6 million (22.2%) and 0.2 million (0.3%) reads were associated with exon-3^′ UTR and exon-5^′ UTR sequences, respectively; 1.4 million reads (19.4%) mapped to intergenic locations and 0.4 million reads (5.5%) mapped to intronic regions (Figure 2). The filtered sequence reads for all deconvoluted libraries were also aligned to the complete genome sequence of the M. bovis AF2122/97 strain (GenBank accession number NC_002945.3) to assess the presence of mycobacterial RNA contamination in individual libraries―a mean of 234 reads (0.004%) mapped to locations in the M. bovis genome. The reads that mapped to the M. bovis genome did not map to the B. taurus reference genome. Following preliminary filtering and quality checks, only sequence reads that mapped to unique locations in the B. taurus reference genome were used for downstream bioinformatics and IPA analyses.

Quantification of the number of reads that exclusively mapped to genes with bovine Ensembl IDs was performed using the HTseq-count software package. This analysis demonstrated that a mean of 5.0 million reads (43.8%) per library mapped to Ensembl gene IDs based on sense strand sequence information, with a mean of 0.4 million sequence reads for each library (3.5%) mapping to Ensembl gene IDs based on antisense strand sequence information; only these two sets of reads were used to separately derive gene expression values for sense and antisense strand transcription, respectively. Of the remaining reads—which were not used to derive gene expression values—a mean of 1,600 reads (0.01%) for each library were associated with multiple Ensembl gene IDs, while a mean of 1.8 million reads for each library (16.3%) were not associated with Ensembl gene IDs (Figure 1B).

Prior to multi-dimensional scaling (MDS) and differential gene expression analysis, density plots (Additional file 2: Figure S1) displaying the number of sequence reads per gene were constructed and analysed. Two RNA-seq libraries comprising a control and M. bovis-infected sample from the same animal (animal number 700; Additional file 1: Table S1) showed skewed distributions compared to all other libraries; consequently these samples were removed from all further downstream analyses. All downstream analyses, including differential gene expression, IPA analysis based on sense and antisense strand data, and technical comparisons between RNA-seq and microarray platforms were performed using data from the remaining six animals (i.e. 12 RNA-seq libraries).

Gene expression and IPA analysis of sense strand transcription

Analysis of the gene coverage using reads mapping to unique locations of the B. taurus reference genome based exclusively on sense strand sequence information, showed that of the 25,669 annotated B. taurus genes, 15,422 genes (60.1%) had at least one sequence read count (i.e. one mapped read) in at least one library. The 15,422 detectable genes were further filtered by removing lowly expressed genes. Consequently, only genes displaying more than one count per million reads for at least three libraries were used for subsequent analyses. This yielded 11,131 genes (43.4% of annotated B. taurus genes) that were suitable for downstream analyses.

Prior to differential gene expression analysis, the data from the 11,131 filtered genes were used for MDS analysis (Figure 3). MDS analysis demonstrated that MDM samples were differentiated according to treatment status (i.e. infected versus non-infected control samples) along dimension 1, while dimension 2 separated the MDM samples according to animal ID.

Statistical analysis of all 11,131 genes that passed the filtering process identified a total of 2,584 differentially expressed genes (adjusted P-value ≤ 0.05), of which 1,392 were upregulated and 1,192 were downregulated in the M. bovis-infected MDM relative to the non-infected control MDM. In addition, expression fold-change values were markedly higher for upregulated genes compared to downregulated genes. A list of the differentially expressed genes based on sense strand data is presented in Additional file 3: Table S2.

Functional categorisation of the 2,584 differentially expressed genes using IPA revealed an enrichment of genes with roles in immunology, cell development and proliferation, and apoptosis (Additional file 4: Table S3). IPA was also used to identify canonical pathways enriched for differentially expressed genes, many of which were had immune and apoptotic functions. (Additional file 5: Table S4). Based on the well-documented role of apoptosis in the host response to mycobacterial infection [25] the Death receptor signalling canonical pathway was overlaid with RNA-seq gene expression results and is presented in Figure 4.

Gene expression analysis of antisense strand read data

Analysis of the RNA-seq reads that mapped to antisense strand gene sequences revealed that 13,796 genes of the 25,669 annotated B. taurus genes had at least one sequence read count (i.e. one mapped read count) for at least one library. Filtering out of lowly expressed genes (as detailed for the sense strand data above) led to the exclusion of 6,954 genes, leaving a total of 6,842 genes for differential expression analysis based of antisense strand data. 757 of the genes were differentially expressed (adjusted P-value ≤ 0.05), of which 558 showed increased expression and 199 showed decreased expression in the M. bovis-infected MDM relative to the non-infected control MDM. As with the sense strand data, the fold-change in expression for the genes showing increased expression was markedly higher than that for the genes decreasing in expression. A list of significant differentially expressed genes derived from the antisense strand expression data is presented in Additional file 6: Table S5. Additional file 7: Figure S2 shows the distribution of sense and antisense reads that mapped to the spectrin, beta, erythrocytic gene (SPTB). As these reads do not map to the exact same locations in the SPTB gene, we are confident that they represent actual antisense transcripts based on the current annotation of the B. taurus reference genome and are not due to technical artefacts introduced during RNA-seq library preparation.

A systematic comparison was performed between significantly differentially expressed genes detected using the sense strand RNA-seq data (n = 2,584) and those identified from the antisense strand data (n = 757). This generated a list of 694 differentially expressed genes that were common to both the sense strand and the antisense strand data sets. The common 694 genes were further subdivided according to the direction of expression (i.e. up- or downregulated) in the M. bovis-infected MDM relative to the control MDM. 520 genes displayed an increase in expression in both the sense and antisense strand data sets, while 173 genes showed decreased expression from both strands. In addition, no gene displayed increased expression from the sense strand and decrease expression from the antisense strand. Finally, one gene showed decreased expression from the sense strand and increased expression from the antisense strand (Figure 5).

Functional categorisation of the 757 differentially expressed genes based on antisense strand data using IPA revealed an enrichment of genes with roles in immunology, cell development and proliferation, and response to infection (Additional file 8: Table S6). Four of the top five ranking GO categories identified from the antisense and sense strand data were identical (Additional file 4: Table S3 and Additional file 8: Table S6). Similarly, IPA analysis of the antisense strand data generated a list of canonical pathways with immune related function (Additional file 9: Table S7). Four of the five top-ranking canonical pathways identified from the antisense and sense strand data were identical (Additional file 5: Table S4 and Additional file 9: Table S7).

Comparison of differential gene expression profiles obtained from RNA-seq and microarray platforms

The total RNA samples extracted from M. bovis-infected and non-infected control MDM samples described here, were previously used by us for gene expression profiling with the Affymetrix® GeneChip® Bovine Genome Array [20]. This array contains 24,027 probe sets representing more than 23,000 transcripts and includes approximately 19,000 UniGene clusters [26].

In order to directly compare gene expression profiles from these previously generated microarray data with the RNA-seq data generated in the current study, a number of steps were performed. Firstly, we removed the same two samples which were excluded from the RNA-seq data analysis (i.e. the non-infected control and M. bovis-infected MDM-extracted RNA from animal number 700), then re-analysed all microarray data using all 12 RNA samples from the remaining six animals. Secondly, we mapped Affymetrix® GeneChip® Bovine Genome Array probes to genes with bovine Ensembl gene IDs based on data from the B. taurus reference genome. Of the 24,072 probe sets represented on the array, 11,790 probe sets passed filtering; with 8,560 unique probe sets mapping to a B. taurus Ensembl gene ID, representing 6,807 unique genes. Prior to differential gene expression analysis, the data from the 11,790 filtered probes was used for MDS analysis (Additional file 10: Figure S3). As with the data from the RNA-seq, the MDS analysis demonstrated that MDM samples were differentiated according to treatment status (i.e. infected versus non-infected control samples) along dimension 1, while dimension 2 separated the MDM samples according to animal ID. It is important to note that joint inspection of Figure 3 and Additional file 10: Figure S3 shows that the samples are spatially arranged in a very similar pattern for the two different gene expression platforms.

Comparison of the number and fold-change of differentially expressed genes identified from RNA-seq and microarray

Subsequent analysis of the microarray data showed that 2,015 unique genes were differentially expressed, with 917 genes and 1,098 genes displaying up- and downregulation, respectively, in the infected MDM relative to the control MDM (FDR ≤ 0.05; Additional file 11: Table S8). Comparison of the differentially expressed genes with Ensembl gene IDs identified from the microarray and RNA-seq sense strand data sets revealed 964 differentially expressed genes displaying the same direction of expression that were common to both the microarray and RNA-seq data sets; these comprised 607 upregulated and 357 downregulated common genes. One gene was found respectively downregulated and upregulated in the microarray and RNA-seq sense strand data sets (Additional file 12: Table S9). Of the remaining differentially expressed genes detected, 1,050 genes (comprising 310 upregulated genes and 740 downregulated genes) were unique to the microarray data set, while 1,619 genes (comprising 784 upregulated genes and 835 downregulated genes) were unique to the RNA-seq data set (Figure 6). The differentially expressed genes that were common to both platforms and displayed the same direction of expression, represent 47.8% (i.e. 964/2,015 genes) and 37.3% (i.e. 964/2,584 genes) of all differentially expressed genes detected by microarray and RNA-seq analysis, respectively.

Spearman rank analysis of the log₂ RPKM values (RNA-seq data) with the log₂ hybridisation intensity signals (microarray data) for all 964 differentially expressed genes that were common to both platforms and had the same direction of expression revealed a correlation coefficient of r = 0.713 for infected MDM and r = 0.660 for control MDM (P ≤ 0.01), indicating that there is a positive relationship between raw signal intensities for each platform/treatment group. Finally, comparison of the log₂ fold-change in gene expression between the two platforms for these 964 genes yielded a correlation coefficient of r = 0.929 (P ≤ 0.001).

To investigate the effect of platform dynamic range on the number differentially expressed genes unique to each platform we plotted and compared the density of the log₂ of the mean hybridisation intensity and the density of the log₂ of the mean CPM for the differentially expressed genes unique to each platform and for the differentially expressed genes common to both platforms. This analysis showed that differentially expressed genes uniquely detected by RNA-seq are largely characterised by low levels of expression for both the control and infected MDM groups. This pattern was not observed for the microarray data (Additional file 13: Figure S4). This analysis shows that the sensitivity of RNA-seq for the detection of lowly expressed transcripts is greater than the microarray.

IPA analysis of differentially expressed genes common to both microarray and RNA-seq

To validate the IPA analysis performed using differentially expressed genes based on RNA-seq sense strand data only, we uploaded the 965 differentially expressed genes common to both RNA-seq and the microarray into IPA (964 genes plus the one gene that displayed reciprocal expression for the two platforms). We hypothesised that the 965 common genes were more likely to reflect true biological changes induced following infection with M. bovis, compared to the differentially expressed genes identified by either platform independently. This approach has been previously used to validate GO category analysis in a study where RNA-seq and microarray data were available for the same samples [28]. Of the 63 GO categories identified for the common differentially expressed genes, 56 (88.9%) and 60 (95.2%) were observed in the RNA-seq and microarray IPA data sets, respectively (Additional file 4: Table S3, Additional file 14: Table S10 and Additional file 15: Table S11). Similarly, of the 198 canonical pathways identified for the common differentially expressed genes, 160 (80.8%) and 166 (83.8%) were observed in the RNA-seq and microarray IPA data sets, respectively (Additional file 5: Table S4, Additional file 16: Table S12 and Additional file 17: Table S13). This analysis suggests that meaningful biological interpretation of data from either RNA-seq or microarrays is possible.

Validation of differentially expressed genes using real time quantitative reverse transcription PCR

For method validation purposes, we quantified a panel of 16 genes via real time quantitative reverse transcription PCR (qRT-PCR) using the same MDM-extracted RNA samples that were used for the RNA-seq and microarray platforms as described above. The methods used for real time qRT-PCR analysis of these genes have been described by us elsewhere [20]. Using these data, we estimate the dynamic range of the real time qRT-PCR platform has a log₂ dynamic range of 15—this is based on C_t value ranges between 20 and 35 for a number of cytokine genes.

Real time qRT-PCR analysis showed that 12 of these genes (CCL4, CCL5, CCL20, CD40, CFB, CXCL2, IL15, IL1B, IL6, IRF1, NFKB2 and TNF) were upregulated, while two genes (AREGB and FOS) were downregulated in the infected MDM (P ≤ 0.05). The remaining two genes, PIK3IP1 and SPRY2 were not differentially expressed based on analysis of the real time qRT-PCR data.

Summary of RNA-seq results

For each individual RNA-seq library analysed here, we obtained a mean of 7.2 million 69 bp reads that mapped to single unique locations in the B. taurus reference genome. This yielded a mean of 496.8 Mb of sequence per individual library. Examination of the gene coverage from these uniquely mapping reads demonstrated that 15,422 genes (based only on sense strand data) from a total of 25,669 annotated bovine genes (60.1%) gave at least one sequence read count in total. However, studies have shown that gene expression analysis using RNA-seq is dependent on sequencing depth and that genes with low sequencing coverage are more susceptible to the generation of false-positive differentially expressed genes. In this regard, we only considered genes displaying more than one count per million reads in at least three libraries for differential gene expression analysis; consequently a total of 11,131 genes (43.4% of the total bovine gene content) based on sense strand data were used in the analysis presented here. Although this stringent threshold lowers the number of genes for which differential expression can be performed, we believe that this criterion is sufficient for quantification and analysis of highly expressed genes with a corresponding reduction in the number of Type I errors due to lowly expressed genes [35–37].

Sense strand gene expression and IPA analyses

Gene expression analysis of the sense strand information in the current study detected a total of 2,584 differentially expressed genes (adjusted P-value ≤ 0.05), of which 1,392 and 1,192 were upregulated and downregulated, respectively, in the M. bovis-infected MDM relative to the non-infected control MDM. This finding contrasts with previous work performed by us and others showing that mycobacterial infection in vitro and in vivo results in a higher number of downregulated genes relative to upregulated genes based on microarray analysis [17, 20, 21, 29, 38]. This discrepancy is most obviously explained by differences in sensitivity and technical biases between the two transcriptomic platforms used. This assertion is supported in the current study by the increased dynamic range of RNA-seq compared to the microarray platform. In agreement with previous studies, however, the expression fold-change values were markedly higher for upregulated genes compared to downregulated genes [17, 19–21].

Analyses performed at the gene level revealed that all the top upregulated and downregulated genes have roles and functions involved in or related to immune response and infectious disease. One of the most highly upregulated genes identified was ADRB3; the protein encoded by this gene has been shown to have functions in carbohydrate metabolism, energy reserve metabolism, positive regulation of the MAPK cascade and regulation of apoptosis [39, 40]. IL17A was also upregulated and has known role in the pro-inflammatory response, expansion and recruitment of innate immune cells, production of defensins and antimicrobial peptides, and linking innate and adaptive immune responses [41]. Also upregulated was SAA3, which was the most statistically significant gene in the data set, and has a major multifunctional role in the acute-phase response [42].

Among the most highly downregulated genes was GABBR2, which has multiple roles and functions, including involvement in pulmonary disorder and regulation of the coughing process [43–45]. In addition, CDH26 was downregulated and cadherins have been shown to have a role in Ca²⁺-dependent cell-to-cell adhesion as well as cell polarisation and cell migration [46]; however, the specific role of CDH26 has not been characterised. The SLC16A12 gene was also detected, a member of a gene family containing genes recently shown to be associated with tuberculosis susceptibility in cattle [47]. FAIM2 was highly downregulated, a gene with a role in the inhibition of apoptosis, a key process associated with mycobacterial infection of macrophages [25, 48]. Finally, SH3TC2 was downregulated, which interacts with members of the RAB11 protein family, leading to regulation of endocytic recycling, an important process during mycobacterial infection [49–51].

IPA was used for identification of significantly over-represented GO categories and canonical pathways. Notably, these analyses demonstrated that the processes of cell death, apoptosis and cell survival were over-represented as GO categories and within canonical pathways levels. These pathways included: Death receptor signalling; Apoptosis signalling; TNFR2 signalling; TNFR1 signalling; PTEN signalling; and TWEAK signalling. Inspection of these pathways revealed that the genes encoding TNF-α and the tumour necrosis factor (ligand) superfamily, member 10 (TNFSF-10) protein are highly upregulated. These proteins can, via their respective receptors—tumour necrosis factor receptor superfamily member 1A (TNFRSF-1A) and TNFRSF-10A/TNFRSF-10B—trigger an activation cascade through Fas (TNFRSF6)-associated via death domain (FADD) and TNFRSF1A-associated via death domain (TRADD) [both genes encoding these proteins were upregulated in the present study]. This cascade concludes with the activation of several caspases that are associated with the induction of apoptosis (upregulation of CASP7 and CASP8 was observed here), further supporting the role of apoptosis in the host response to mycobacterial infection [25].

IPA analysis also demonstrated the importance of gene products involved with immune cell communication and chemotaxis of immune cells, including the canonical pathways Communication between innate and adaptive immune cells and The role of cytokines in mediating communication between immune cells. For example, it is well established that macrophage recognition of bacterial pathogen-association molecular patterns (PAMPs) via specific Toll-like receptors (TLRs) [upregulation of TLR3; TLR5; and TLR10 and downregulation of TLR9 was observed here] leads to increased expression of a wide range of NF-κB-inducible cytokines and chemokines. In this regard, CCL4, CCL5, CXCL10 and IL8 were upregulated in the present study. Indeed, the overall pattern of immune gene expression changes observed here using RNA-seq data was consistent with previous work from our group using microarray technology [20].

Antisense strand transcript expression

The present study identified 0.4 million sequence reads that mapped to the antisense strand and it is important to note that the phenomenon of natural antisense transcription has been observed previously in mammals [52, 53]. We identified 6,842 putative natural antisense transcripts (NATs) that were suitable for differential expression analysis of the M. bovis-infected MDM relative to the non-infected control MDM. This represents 26.7% of the 25,669 currently annotated B. taurus genes, which corresponds to results obtained from previous studies of antisense strand transcription on different mammals including cattle [54–58]. Of the 6,842 NATs detected, 757 were significantly differentially expressed; these comprised 558 that were upregulated and 199 that were downregulated. NATs have been shown to have varying regulatory roles in eukaryotic cells, including transcriptional and post-transcriptional control, splicing event regulation, allele-specific transcript expression, RNA editing and RNA translocation [59–64]. Further analyses identified 694 genes that were differentially expressed from both the sense and antisense strand.

IPA analysis of these 757 antisense strand transcripts revealed enrichment of genes involved in several immune processes including the inflammatory response, pattern recognition receptor signalling, cell death and apoptosis, immune cell movement and interaction, and antigen presentation. These findings suggest that NATs may play a significant role in regulating the host macrophage response following infection with pathogenic agents, including mycobacteria [52, 65].

Comparison with microarray and real time qRT-PCR

The MDM samples infected with M. bovis used for the present study had previously been analysed using the Affymetrix® GeneChip® Bovine Genome Array. Consequently, it was possible to directly compare gene expression results from microarray and sense strand RNA-seq data. The percentage of overlapping genes between the two platforms was estimated at 47.8% for the microarray and 37.3% for the RNA-seq data. Similar estimates of correspondence between RNA-seq and microarrays platforms have been reported in studies of the rat and Candida parapsilosis transcriptomes [28, 66]. However, it should be noted that compared to the study of Su et al.[28] we observed a greater correlation between the fold-change in expression for: (1) all common genes that passed filtering detected by both platforms (6,183 genes; Spearman correlation coefficient of 0.858; P ≤ 0.001), and (2) the differentially expressed genes common to both platforms (965 genes; Spearman correlation coefficient of 0.935; P ≤ 0.001).

There are a number of possible reasons for the differences observed between the two gene expression platforms used in the current study. Firstly, there are differences in the dynamic range of the two platforms: our results show that RNA-seq has a larger dynamic range and sensitivity than the microarray; similar results have recently been obtained by Chen and colleagues [27]. Secondly, the stringent expression thresholds that we have applied to the RNA-seq data in order to reduce the number of false-positive differentially expressed genes, have also reduced the number of genes common to both platforms. Thirdly, the probes on the Affymetrix® GeneChip® Bovine Genome Array are based on sequences from the 3^′ end of genes and are therefore 3^′ biased [67]; RNA-seq reads are more randomly distributed across gene transcripts. Fourthly, microarrays are susceptible to cross-hybridisation, particularly among members for the same gene family that can result in elevated false-positive rates [30–32]. Fifthly, there are differences in the statistical models used to detect differentially expressed genes for the two platforms (for example, analysis of the microarray data involve moderated t-tests; analysis of the RNA-seq data using edgeR involved a negative binomial distribution) [28, 68–71]. Sixthly, differences in the bovine genome resources used to design the Affymetrix® GeneChip® Bovine Genome Array (May 2005, Gene Expression Omnibus platform accession number GPL2112) and the B. taurus reference genome (Btau 4.0.63 genome release) used to analyse the RNA-seq data may also impact on the analysis as recently proposed in a study of the yeast transcriptome [72]. Seventhly, it has been demonstrated that increased sequencing depth also contributes to greater correspondence between RNA-seq and microarray platforms, such as that recently observed in a study of the Saccharomyces cerevisiae transcriptome [36, 72].

It is important to note, however, that despite the moderate concordance in the number of differentially expressed genes common to both transcriptomic platforms, we did observe a high concordance between the GO categories and canonical pathways identified by IPA analysis. This indicates that the RNA-seq and microarray platforms both provide gene expression data that can be used for meaningful biological interpretation.

MDM preparation, MDM infection and RNA purification

The materials and methods used to isolate, purify and infect bovine MDM with M. bovis have been previously described in detail by us [20]. Briefly, MDM were purified from peripheral blood mononuclear cells prepared from whole blood extracted from seven age-matched (four-years old) Holstein-Friesian females (Additional file 1: Table S1). All seven animals were selected from a herd without a recent history of bovine tuberculosis infection, which was confirmed using the single intradermal tuberculin test. MDM were cultured over an eight-day period with routine changing of culture media every two days. MDM were visually inspected using microscopy, counted and seeded in 24-well tissue culture plates at a density of 2 × 10⁵ cells per well prior to infection.

For the MDM infections, 4 × 10⁵ M. bovis cells (as determined from bacterial cell counts performed using a Petroff Hausser chamber and confirmed by colony-forming unit [cfu] counts) were added to each tissue culture plate well (MOI 2:1). The M. bovis M2137 strain bearing the SB0142 spoligotype was used for the infection experiments and non-infected control MDM received culture media only. Both non-infected control and infected MDM were prepared in adjacent duplicate tissue culture plate wells. The culture media in all wells for both non-infected control and infected MDM was replaced 2 h post-infection with fresh culture media and plates were then re-incubated at 37°C, 5% CO₂ until the MDM were harvested. Culturing of M. bovis and the MDM infections were performed in a Biosafety Containment Level 3 (CL3) laboratory.

Infected and non-infected control MDM were harvested using RLT buffer from the RNeasy Mini kit (Qiagen Ltd., Crawley, UK) supplemented with 1% β-mercaptoethanol (Sigma-Aldrich Ireland Ltd., Dublin, Ireland) 24 h post-infection. For each control and treatment, MDM lysates from duplicate culture plate wells were pooled and stored at -80°C until required for RNA extraction. All RNA extractions were performed in the CL3 laboratory using an RNeasy kit incorporating an on-column DNase treatment step according to the manufacturer’s instructions (Qiagen). RNA quantity and quality was ascertained using a NanoDrop™ 1000 spectrophotometer (Thermo Fisher Scientific Ltd., Waltham, MA, USA) and an Agilent 2100 Bioanalyzer with the RNA 6000 Nano LabChip kit (Agilent Technologies Ltd., Cork, Ireland). All samples displayed a 260/280 ratio greater than 2.0 and RNA integrity numbers ≥ 8.5 [73].

Strand-specific RNA-seq library preparation

In total, 14 strand-specific RNA libraries for high-throughput sequencing were prepared (seven libraries for each treatment: M. bovis-infected and control samples) using 200 ng of total RNA. Total RNA was first heated at 65°C for 5 min to disrupt any secondary structure. Purification of poly(A) RNA was performed using a Dynabeads® mRNA DIRECT™ Micro Kit according to the manufacturer’s instructions (Invitrogen™/Life Technologies Ltd., Paisley, UK). Purified poly(A) RNA was then fragmented using 1× RNA Fragmentation Reagent (Ambion®/Life Technologies Corporation, Warrington, UK) for 5 min at 70°C and precipitated using 68 mM sodium acetate pH 5.2 (Ambion), 227 ng/μl glycogen (Ambion) and 30 μl of 100% ethanol (Sigma-Aldrich Ltd., Dublin, Ireland). Pellets were washed with 80% ethanol, air-dried for 10 min at room temperature and re-suspended in 10.5 μl DNase- and RNase-free water.

Synthesis of first strand cDNA was performed by incubating fragmented RNA with 261 mM Random Hexamer Primers (Invitrogen), 1× first strand buffer (Invitrogen); 10 mM DTT (Invitrogen); 0.5 mM dNTPs; 20 U RNaseOUT™ Recombinant Ribonuclease Inhibitor; and 200 U SuperScript® II Reverse Transcriptase (Invitrogen) at 25°C for 10 min, at 42°C for 50 min, and 70°C for 15 min. First strand synthesis reaction mixtures were purified using MicroSpin™ G-50 columns according to the manufacturer’s instructions (GE Healthcare UK Ltd., Little Chalfont, Buckinghamshire, UK).

Second strand cDNA synthesis, involving the incorporation of uracil, was performed by adding the first strand cDNA synthesis reaction to a second strand reaction mix consisting of 0.065× first strand buffer (Invitrogen); 1× second strand buffer (Invitrogen); a dNTP mix consisting of a final concentration of 0.3 mM dATP, dCTP, dGTP (Sigma-Aldrich) and 0.3 mM dUTP (Bioline Reagents Ltd., London, UK); 1 mM DTT (Invitrogen); 2 U RNase H (Invitrogen) and 50 U E. coli DNA Polymerase I (Invitrogen). Reactions were incubated at 16°C for 2.5 h. The double stranded cDNA was subsequently purified using a QIAquick PCR Purification kit (Qiagen) according to the manufacturer’s instructions and eluted in 30 μl of the provided elution buffer.

Blunt-end repair of cDNA was performed in a 100 μl reaction containing 1× T4 DNA ligase buffer with 10 mM dATP (New England Biolabs® Inc., MA, USA), 0.4 mM of each dNTP (Invitrogen), 15 U T4 DNA polymerase (New England Biolabs), 5 U DNA Polymerase I Large [Klenow] Fragment (New England Biolabs) and 50 U T4 polynucleotide kinase (New England Biolabs). Reactions were incubated at 20°C for 30 min and the cDNA was then purified using a QIAquick PCR Purification Kit (Qiagen) according to the manufacturer’s instructions and eluted in 32 μl of the provided elution buffer.

To facilitate Illumina® GA adaptor ligation, a single ‘A’ base was added to the 3^′ ends of the blunt-end repaired cDNA samples. 32 μl of purified phosphorylated blunt end-repaired cDNA was included in a final 50 μl reaction mixture containing: 1× Klenow fragment buffer (New England Biolabs); 0.2 mM dATP (Invitrogen), and 15 U Klenow fragment with 3^′-to-5^′ exonuclease activity (New England Biolabs). Reactions were incubated at 37°C for 30 min, after which cDNA was purified using a QIAquick MinElute Kit (Qiagen) according to the manufacturer’s instructions and eluted in 21 μl of the provided elution buffer.

Illumina® RNA-seq adaptor ligation reactions (50 μl volumes) involved incubation of 21 μl of phosphorylated blunt-ended cDNA containing a 3^′-dATP overhang with 1× Quick DNA ligase buffer (New England Biolabs); 30 nM custom indexed single-read adaptors [Additional file 1: Table S1] and 15 U T4 DNA ligase (Invitrogen). Reaction mixes were incubated at room temperature for 15 min and purified using a QIAquick MinElute Kit according to the manufacturer’s instructions (Qiagen) and eluted in 10 μl of the provided elution buffer. Adaptor-ligated cDNA was gel-purified using 2.5% agarose gels stained with 1 μg/ml ethidium bromide (Invitrogen). Gels were electrophoresed at 100 Volts using 1× TAE buffer (Invitrogen) for 75 min at room temperature. Size fractionated bands corresponding to 200 bp (+50 bp) were excised from each sample and purified using a QIAquick Gel Extraction kit (Qiagen) according to the manufacturer’s instructions and eluted in 30 μl of elution buffer.

To generate strand-specific RNA-seq libraries, the second strand of the gel-purified adapter-ligated cDNA containing uracil was digested enzymatically in 30 μl reaction volumes containing 1× Uracil-DNA Glycosylase buffer and 1 U Uracil-DNA Glycosylase (Bioline). Reactions were incubated at 37°C for 15 min followed by 94°C for 10 min.

PCR enrichment amplifications (50 μl) containing 15 μl of second strand-digested, adaptor-ligated cDNA; 1× Phusion® High-Fidelity DNA polymerase buffer (New England Biolabs); 334 nM each Illumina® PCR primer (Illumina® Inc., San Diego, CA, USA); 0.4 mM each of dATP, dCTP, DGTP and dTTP (Invitrogen) and 1 U Phusion® High-Fidelity DNA polymerase (New England Biolabs). PCR amplification reactions consisted of an initial denaturation step of 98°C for 30 seconds, 18 cycles of 98°C for 10 seconds, 65°C for 30 seconds and 72°C for 30 seconds, followed by a final extension step of 72°C for 5 min. PCR products were visualised following electrophoresis on a 2% agarose gel stained with ethidium bromide (0.6 μg/ml; Invitrogen) and purified to remove PCR-generated adaptor-dimers using an Agencourt AMPure XP kit (Beckman Coulter Genomics, Danvers, MA, USA) according to the manufacturer’s instructions with final elution in 30 μl of 1× TE buffer.

All RNA-seq libraries were quantified using a Qubit® Fluorometer and Qubit® double stranded DNA High Sensitivity Assay Kit (Invitrogen). RNA-seq library quality was assessed using an Agilent Bioanalyzer and Agilent High sensitivity DNA chip (Agilent) and confirmed that library insert sizes were ~200-250 bp for all individual libraries. Individual RNA-seq libraries were standardised and pooled in equimolar quantities (10 μM for each individual library). The quantity and quality of the final pooled library was assessed as described above prior to sequencing.

Cluster generation and sequencing of the pooled RNA-seq library was carried out on an Illumina® Cluster Station and Illumina® Genome Analyzer IIx sequencer according to the manufacturer’s instructions (Illumina). The pooled library was sequenced as single-end read 84-mers using Illumina® version 4.0 sequencing kits and the standard Illumina® Genome Analyzer II_x pipeline. The Illumina® Sequencing Control Software version 2.9 and Real Time Analysis version 1.9 software packages were used for real-time tracking of the sequencing run, real-time image processing, the generation of base intensity values and base calling. These RNA-seq data have been deposited in the NCBI Gene Expression Omnibus (GEO) database with experiment series accession number GSE45439.

Statistical analysis of RNA-seq data

Sequence reads obtained from seven lanes of the Illumina flow cell were deconvoluted into 14 individual libraries using the unique indexed barcoded adapters. A Perl script was used to screen adapter artefacts, removing any reads containing a full-length match to the 33 nucleotide Illumina® adapter sequence, allowing up to four mismatches. Read quality was then assessed using FastQC software [version 0.9] (http://www.bioinformatics.bbsrc.ac.uk/projects/fastqc), revealing low Phred scores at the 3^′ end of sequence reads. Consequently, nine nucleotides were trimmed from the 3^′ ends of all sequence reads using the FASTX Toolkit [version 0.0.13] (http://hannonlab.cshl.edu/fastx_toolkit) generating usable reads of 69 nucleotides.

Deconvoluted quality-checked sequence reads were aligned to the B. taurus reference genome (Btau 4.0.63 genome release) with the TopHat splice junction mapper [version 1.3.0] [74, 75], which aligns reads using the Bowtie aligner (version 0.12.7). The Bowtie alignment procedure was configured for strand-specific libraries and to filter only unique hits to the reference genome sequence. In addition, all sequence reads were aligned to the M. bovis AF2122/97 chromosome, complete genome (accession number NC_002945.3) as detailed above. To obtain raw counts per transcript, HTSeq package [version 0.5.3p1] (http://www-huber.embl.de/users/anders/HTSeq/doc/overview.html) was used on alignment files in BAM format. HTSeq-count was used with option overlap resolution mode set to intersection non-empty. Counts of uniquely-mapped reads were obtained for all bovine Ensembl genes and transcripts, with separate counts obtained for sense and antisense DNA strands.

Once raw counts were obtained from HTseq-count, analysis of differential expression was performed in the R statistical programming environment [76] using the edgeR (version 2.2.5) Bioconductor package [70, 77] and lattice [version 0.19-30] (http://lattice.r-forge.r-project.org). First, quality checks were performed by plotting the density of counts per feature (a feature being a gene or transcript generated respectively from the sense or antisense strand) for each sample and also by generating a multidimensional scaling plot of the RNA-seq data. The edgeR software package was used to determine differential expression using the paired-sample statistical test. Filtering of lowly expressed features was performed by retaining only features with at least one count per million in three or more libraries. A normalisation factor was calculated using the default trimmed mean of M values (TMM) method [70, 77], and the dispersion parameter for each feature was estimated as the Cox-Reid common dispersion method in the edgeR package. Differential expression was evaluated by fitting a negative binomial generalized linear model for each feature and then adjusting the P-value for multiple testing using the Benjamini-Hochberg correction [78] with a false discovery rate (FDR) of 0.05.

For RNA-seq dynamic range analysis we used two different normalisation strategies: (1) RPKM values generated by the edgeR package; and (2) quantile normalisation of RNA-seq data performed using the Linear Models for Microarray Data (LIMMA) package [71]. RPKM and quantile-normalised counts were not used for differential expression analysis.

Reanalysis and reannotation of Affymetrix® GeneChip® Bovine Genome Array data

The RNA samples analysed in the current study were previously analysed by us using the Affymetrix® GeneChip® Bovine Genome Array [20]. These microarray data have been deposited in the NCBI GEO database with experiment series accession number GSE33309. In order to compare gene expression profiles from these samples using RNA-seq and microarray technologies, we first reanalysed the microarray data using a series of Bioconductor packages.

Firstly, quality control analysis was performed on all microarrays using the Simpleaffy software package [79]. Then the raw microarray expression values were normalised using the gcRMA package [80]. The raw data followed an additional normalisation step using the Factor Analysis for Robust Microarray Summarization (FARMS) algorithm to remove probe sets with high noise:signal ratios [81], these normalised data were then further subjected to filtering for informative probe sets using informative/non-informative calls (I/NI-calls) package [82]. The obtained I/NI filtering list was then applied to the gcRMA normalised data to remove all non-informative probes. Differentially expressed genes were identified from the filtered gcRMA normalised data using the LIMMA package. The Benjamini-Hochberg multiple testing correction method [78] was applied to all differentially expressed genes to minimise the FDR and adjusted P-values for differentially expressed genes were calculated. Probe sets displaying differential expression between control and infected samples were annotated to Ensembl gene IDs using biomaRt package [83].

Supplementary statistical analyses were performed using the SPSS software package (version 20; http://www-01.ibm.com/software/analytics/spss).

IPA analyses

Ingenuity® Systems Pathway Analysis (IPA, Ingenuity Systems, Redwood City, CA, USA; release date November 2012) was used to identify canonical pathways and functional processes of biological importance within the list of differentially expressed genes identified with RNA-seq and microarray platforms. The Ingenuity® Knowledge Base contains the largest database of manually-curated and experimentally-validated physical, transcriptional and enzymatic molecular interactions. Furthermore, each interaction in the Ingenuity® Knowledge Base is supported by previously published information.

Functional analysis of genes was performed using IPA to characterise the GO categories of differentially expressed genes between the control and M. bovis-infected MDM. For this, IPA performed an over-representation analysis that categorises differentially expressed genes into functional groups using the Ingenuity® Knowledge Base. Each category in IPA is ranked based on the number of differentially expressed genes in each functional group. The right-tailed Fisher’s exact test was used to calculate a P-value for each GO category assigned to differentially expressed genes.

Ingenuity® Systems Pathway Analysis contains a large library of known canonical pathways that were overlaid with the differentially expressed genes to identify major biological pathways associated with M. bovis infection in MDM. The significance of the association between differentially expressed genes and each canonical pathway was assessed using two methods. Firstly, a ratio was estimated from the number of molecules from the differentially expressed gene data set that map to each pathway, compared to the total number of molecules that map to the canonical pathway based on the reference gene list; and secondly, a Fisher’s exact test that generates a P-value for the assignment of the differentially expressed genes to a particular canonical pathway compared to the reference gene list. Canonical pathways were then overlaid with the expression values of the differentially expressed genes.