Analysis of the Bovine Monocyte-Derived Macrophage Response to Mycobacterium avium Subspecies Paratuberculosis Infection Using RNA-seq

Johne’s disease, caused by infection with Mycobacterium avium subsp. paratuberculosis, (MAP), is a chronic intestinal disease of ruminants with serious economic consequences for cattle production in the United States and elsewhere. During infection, MAP bacilli are phagocytosed and subvert host macrophage processes, resulting in subclinical infections that can lead to immunopathology and dissemination of disease. Analysis of the host macrophage transcriptome during infection can therefore shed light on the molecular mechanisms and host-pathogen interplay associated with Johne’s disease. Here, we describe results of an in vitro study of the bovine monocyte-derived macrophage (MDM) transcriptome response during MAP infection using RNA-seq. MDM were obtained from seven age- and sex-matched Holstein-Friesian cattle and were infected with MAP across a 6-h infection time course with non-infected controls. We observed 245 and 574 differentially expressed (DE) genes in MAP-infected versus non-infected control samples (adjusted P value ≤0.05) at 2 and 6 h post-infection, respectively. Functional analyses of these DE genes, including biological pathway enrichment, highlighted potential functional roles for genes that have not been previously described in the host response to infection with MAP bacilli. In addition, differential expression of pro- and anti-inflammatory cytokine genes, such as those associated with the IL-10 signaling pathway, and other immune-related genes that encode proteins involved in the bovine macrophage response to MAP infection emphasize the balance between protective host immunity and bacilli survival and proliferation. Systematic comparisons of RNA-seq gene expression results with Affymetrix^® microarray data generated from the same experimental samples also demonstrated that RNA-seq represents a superior technology for studying host transcriptional responses to intracellular infection.

Introduction

Johne’s disease, caused by infection with Mycobacterium avium subsp. paratuberculosis (MAP) is a chronic granulomatous enteritis of ruminants, both domestic and wild, including cattle, sheep, deer, and other mammalian species (1). Furthermore, there is some evidence, albeit contentious, suggesting that infection with MAP may be associated with Crohn’s disease in humans (2–4). While prevalence figures of Johne’s disease in cattle are difficult to determine – due, in part, to limited sensitivity and specificity of MAP diagnostic tests – current estimates in European countries vary from 31 to 71% (5–8). In the United States, Johne’s disease is estimated to cost the economy between $200 million and $1.5 billion annually, with that figure rising concurrently with herd-level MAP prevalence (9, 10).

The primary route of MAP transmission is believed to be fecal-oral or through ingestion of infected colostrum (11, 12). Once internalized, infectious bacilli cross the intestinal mucosa by penetrating specialized microfold cells (M cells) or enterocytes, which are located in the epithelium lining of the dome areas of Peyer’s patches (13–15). The bacilli then traverse the M cells by transcytosis and migrate to the basolateral side of the cell where they are recognized and phagocytosed by intestinal macrophages. Macrophage recognition of MAP bacilli is mediated by host pathogen recognition receptors (PRRs), including cell-surface Toll-like receptors (TLRs) and intracellular NOD-like receptors (NLRs) (16, 17); indeed, it has been demonstrated that TLR2, TLR4, and NOD2 can independently recognize MAP cellular components (18). Infected macrophages secrete pro-inflammatory cytokines, such as IL-1B and TNF, which activate an early protective T_H1 response characterized by the release of IFN-γ from T-cells. IFN-γ activates the antimicrobial mechanisms of the macrophage that destroys the internalized pathogen and also induces the development of granulomas that actively contain infection in the majority of animals such that clinical signs do not usually manifest (19–21).

The outcome of MAP infection is dependent on the interaction between infected macrophages and T-cells; progression to clinical infection is thought to develop in animals that fail to eradicate the pathogen with a concomitant shift in the immune system from a protective cellular response to a non-protective humoral response. Consequently, both humoral and cellular immune responses can exist simultaneously in infected individuals and it is possible for MAP bacilli to latently infect animals by persisting in host macrophages for prolonged periods and later become reactivated if, for example, the animal subsequently becomes immunosuppressed (22). MAP has the capacity to survive and subvert the macrophage response to ensure its survival and replication (20, 21, 23, 24). In general, the interactions between the macrophage and MAP upon infection are comparable to those observed for other pathogenic mycobacteria such as M. tuberculosis and M. bovis (22). In this context, MAP prevents phagosome maturation, thus facilitating bacterial survival in phagosomes, which in turn provide a niche for further bacterial growth (25). The mechanisms used by MAP to do this are complex but primarily involve the modulation of various cell signaling pathways through interaction with cell membrane receptors, inhibiting phagosome acidification and phagolysosome fusion, and reducing antigen presentation to the immune system (26). MAP, in common with other mycobacterial pathogens also subverts cell death processes, particularly apoptosis to inhibit antigen presentation and the subsequent development of an effective immune response (25). It has also been suggested that inhibition of apoptosis may contribute to the large numbers of infected macrophages that persist in affected tissues (10, 25). Persistence of MAP in macrophages underlies the progression to clinical disease, which is characterized by immunopathology, proliferation of the pathogen, dissemination infection through the host, and ultimately fecal shedding of the pathogen from the host, thus maintaining the cycle of infection (11, 12, 27).

Through modulation and subversion of the bovine host macrophage, MAP promotes its short- and long-term survival. Therefore, analysis of the macrophage transcriptome in response to MAP infection can shed light on the cellular processes underlying pathogen–macrophage interactions and how the perturbation of these pathways is associated with the pathogenesis of Johne’s disease. In recent years, RNA sequencing (RNA-seq) has provided unprecedented opportunities for gene expression analysis of host response to infection, including unbiased whole-transcriptome profiling, sense and antisense transcription analysis, the characterization of new classes of RNA, and the identification of novel mRNA splice variants (28, 29).

Previously, we used the Affymetrix^® GeneChip^® Bovine Genome Array to study host gene expression in RNA extracted from MAP-infected and non-infected control bovine monocyte-derived macrophages (MDM) across a 24 h time course (30). Our analysis revealed a marked reduction in the number of differentially expressed (DE) genes at the 24 h time point compared to the two earlier infection time points; indeed, these results indicated that majority of transcriptional changes induced by infection occur within the first 6 h of infection, with differential gene expression having largely abated 24 h post-infection (hpi). Consequently, for the present study, we describe analysis of the same RNA samples from the 2 and 6 hpi time points using RNA-seq to enhance detection of host macrophage mRNA transcripts and molecular pathways perturbed and modulated by MAP infection.

Materials and Methods

Ethics Statement

All animal procedures were carried out according to the provisions of the Cruelty to Animals Act (Irish Department of Health and Children license number B100/3939) and ethical approval for the study was obtained from the UCD Animal Ethics Committee (protocol number AREC-P-07-25-MacHugh).

Animals

Seven age-matched (4-year old) Holstein-Friesian females were used for this study and have previously been described by our group. These animals had been maintained under uniform housing conditions and nutritional regimens at the UCD Lyons Research Farm (Newcastle, County Kildare, Ireland). The animals did not have a recent history of Johne’s disease and were also negative for infection with M. bovis (30).

MDM Preparation and Infection and RNA Purification

The methods used to isolate, purify, and infect bovine MDM with MAP have been previously described by our group (29–32). MDM from seven age-matched, female Holstein-Friesian cattle were infected in vitro with a clinical isolate of MAP [multiplicity of infection (MOI) of 2 bacilli:1 MDM] and parallel non-infected control MDM samples were also generated.

Total RNA was extracted from each individual sample and purified individually at 0, 2, and 6 hpi and used to prepare pooled strand-specific RNA-seq libraries as previously described by us (29, 33). RNA was extracted using an RNeasy kit incorporating an on-column DNase treatment step (Qiagen Ltd., Crawley, UK) according to the manufacturer’s instructions. The quantity and quality of the RNA was assessed using a NanoDrop™ 1000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA) and an Agilent 2100 Bioanalyzer with an RNA 6000 Nano LabChip kit (Agilent Technologies Ltd., Cork, Ireland). A260/280 ratios >2.0 and RNA integrity numbers (RINs) >8.5 were obtained for all total RNA samples purified across the infection time course.

Strand-Specific RNA-seq Library Preparation

The protocol used for RNA-seq library preparation was adapted from a protocol previously published by our group (29). Thirty-five strand-specific Illumina^® RNA-seq libraries were generated (seven libraries for the MAP-infected and control groups at the 2 and 6 hpi time points and seven 0 h time point control samples) using 150–200 ng of total RNA. Samples were heated at 65°C for 5 min to disrupt RNA secondary structure and purification of poly(A)+ RNA was performed using the 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 Ltd., 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). The RNA pellets obtained were then washed with 80% ethanol, air-dried for 10 min at room temperature and then re-suspended in 10.5 μl of DNase- and RNase-free molecular biology-grade H₂O.

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™) using the following temperature profile: 25°C for 10 min, 42°C for 50 min, and 70°C for 15 min. First strand synthesis reaction mixtures were then purified using MicroSpin™ G-50 columns according to the manufacturer’s instructions (GE Healthcare UK Ltd., 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 solution 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 samples was performed in 100 μl reaction volumes containing 1× T4 DNA ligase buffer with 10 mM dATP (New England Biolabs^® Inc., Ipswich, 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.

Illumina^® RNA-seq adaptor ligation reactions (50 μl volumes) were performed using 21 μl of each of the phosphorylated blunt-ended cDNA (with 3′-dATP overhangs) samples and 1× Quick DNA ligase buffer (New England Biolabs^®); 30 nM custom indexed single-read adaptors (Table S1 in Supplementary Material) 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 V 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. For generation of strand-specific RNA-seq libraries, the second strand of the gel-purified adapter-ligated cDNA containing uracil was enzymatically digested in 30 μl reaction volumes containing 1× Uracil-DNA Glycosylase buffer and 1 U Uracil-DNA Glycosylase (Bioline). These reactions were incubated at 37°C for 15 min followed by 94°C for 10 min.

PCR enrichment amplifications (25 μl) containing 9 μ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 were performed with the following temperature cycling profile: 98°C initial denaturation for 30 s; 18 cycles of 98°C for 10 s, 65°C for 30 s, and 72°C for 30 s; and 72°C final extension step for 5 min. PCR products were visualized 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 Colter 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 (Invitrogen™). RNA-seq library quality was assessed using an Agilent Bioanalyzer and Agilent High sensitivity DNA chip (Agilent) and confirmed that insert sizes were 200–250 bp for all individual libraries. Individual RNA-seq libraries were standardized 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 libraries was performed on an Illumina^® HiSeq 2000 sequencer according to the manufacturer’s instructions. These RNA-seq data have been deposited in the NCBI Gene Expression Omnibus (GEO) database with experiment series accession number GSE62048.

Bioinformatics and Statistical Analysis of RNA-seq Data

All of the bioinformatics pipeline information and associated scripts used for computational analyses are available in a GitHub repository at https://github.com/mauracasey/RNA-sequencing. These analyses were performed on a 32-node Compute Server running Linux Ubuntu (version 12.04.2) with 256 GB of RAM and 24 TB of hard disk drive storage.

Initial quality checks were performed on each of the raw reads data files using the FastQC software (version 0.10.1)¹ to determine the most appropriate read quality filtering methodology. Consequently, a custom perl script was used to deconvolute sequence reads obtained from the flow cell into 35 individual libraries using the indexed barcoded adapters (the script was optimized to work with single-end reads and a six nucleotide barcode at the 5′-end of each read).

For initial sequence adapter removal and quality filtering, appropriate parameters were used with the custom perl script to filter out reads containing adapter sequence (allowing up to three mismatches) and reads below a sequence quality threshold (discard reads with more than 25% bases with a phred score <20); all reads were also trimmed of 20 nucleotides at the 3′-ends.

The FastQC package was used to further assess the filtered individual fastq files, revealing that no further filtering steps were required. The STAR RNA-seq aligner software package (version 2.3.0) (34) was used to align filtered sequence reads to the most recent version of the Bos taurus reference genome [UMD3.1.73; (35)]. Aligned sequence reads in individual SAM files were then used for a final FastQC quality check step to detect quality score biases in the aligned reads and all samples successfully passed.

The featureCounts tool, which is part of Subread software package (36, 37), was used to perform count summarization of sense genes. Reads were assigned to a gene if they were not multi-hit reads and if the mapped location was associated with a unique gene on the sense strand. Differential gene expression analysis was performed using the Bioconductor edgeR package (38) with the gene raw counts obtained from featureCounts. The BioMart tool was used first for gene annotation with Ensembl gene IDs (39). Ribosomal RNA genes were filtered out and lowly expressed genes were also removed with a minimally set threshold of one count per million (CPM) in at least seven individual libraries (the choice of seven libraries is based on the sample size of each treatment group) (38). For each library, a normalization factor was calculated based on RNA composition among libraries (computed using trimmed mean of M-values). For the present study, at this stage the seven 0 h control samples were removed from the data set and not used for any subsequent bioinformatics, differential gene expression, or downstream data analyses.

Using the edgeR package (38), DE genes between MAP-infected versus non-infected control MDM samples for each time point post-infection (2 and 6 hpi) were obtained using paired-sample statistics by fitting a negative binomial generalized linear model to each gene. Multiple-testing correction was performed using the Benjamini–Hochberg method (40) with a false discovery rate (FDR)-adjusted threshold of ≤0.05.

Functional Analyses of DE Genes Obtained Using RNA-seq

The RNA-seq DE gene lists obtained for each time point post-infection were used for downstream systems analysis to identify important cellular pathways with the Ingenuity^® Systems Pathway Analysis Knowledgebase (IPA²;Summer Release, June 2014). This approach was used to identify canonical pathways that were overrepresented based on the list of DE genes at each of the two time points post-infection using Fisher’s exact test (FDR-adjusted P value threshold ≤0.05).

The GOseq Bioconductor package (41) was used to determine gene ontology (GO) biological process functions that were enriched based on the RNA-seq DE gene lists obtained for each time point post-infection (Bonferroni-adjusted P value threshold ≤0.05).

Comparative Analysis of Microarray Data

The raw microarray data generated from the 35 total RNA samples used for the RNA-seq DE gene and downstream analyses (MDM from seven animals at 0 h, 2, and 6 hpi with the corresponding control samples) were retrieved from the NCBI GEO repository (42) with the accession number GSE35185 (30). The Affymetrix^® GeneChip^® Bovine Genome Array used to generate these data contains 24,072 probe sets representing more than 23,000 gene transcripts. The retrieved microarray data were then analyzed with a number of different Bioconductor packages (43) using the UMD3.1.73 build of the bovine genome (35). The Factor Analysis for Robust Microarray Summarization (FARMS) algorithm was used to normalize the microarray data (44) and these normalized data were then filtered for informative probes sets using the FARMS informative/non-informative (I/NI) calls unsupervised feature selection method (45).

To compare data generated using the two different gene expression technologies, microarray probe sets were annotated with bovine Ensembl gene IDs from the B. taurus reference genome build used to annotate the RNA-seq data [UMD3.1.75; (35)] using the Bioconductor biomaRt package (39). DE genes were detected between experimental groups using the Linear Models for Microarray Data (LIMMA) Bioconductor package (46). A Benjamini–Hochberg multiple-testing correction of P ≤ 0.05 was used for all DE genes (40) and the Euclidean distance was used as the distance metric for MDS plotting.

Results

Preliminary Analysis and Summary Statistics for RNA-seq Data

The 35 RNA-seq libraries used for the present study were sequenced across six lanes of an Illumina^® HiSeq 2000 sequencing apparatus and generated mean values per library of 26.72 million raw reads, of which 20.02 million reads (74.94%) remained after adapter sequence and poor quality reads filtering (Figure S1A in Supplementary Material). Alignment of the filtered RNA-seq reads to the B. taurus reference genome (UMD3.1.73) yielded mean values per library of 16.19 million reads (80.85%) mapping to unique locations in the bovine genome, 1.66 million reads (8.31%) mapping to multiple locations in the genome, and 2.17 million reads (10.84%) not mapping to any genome location (Figure S1B in Supplementary Material). Further analysis, focusing on the uniquely mapping reads demonstrated that a mean of 11.91 million reads (73.60%) per library were assigned to annotated sense regions of the genome. Only these sequence reads were then used to calculate raw counts per sense gene and for downstream differential gene expression and systems biology analyses. In addition, a mean value per library of 4.27 million reads (26.40%) could not be assigned to annotated genome locations or were assigned to overlapping annotated genomic regions (Figure S1C in Supplementary Material). The detailed number of reads per individual RNA-seq library at each stage of the analysis is provided in Table S1 in Supplementary Material.

Analysis of the gene coverage based solely on sense sequence information, demonstrated that of the 24,616 B. taurus genes annotated in Ensembl (release 73), 17,571 of these genes (71.4%) had at least one sequence read count (i.e., one mapped read) in at least one of the 35 individual RNA-seq libraries. These 17,571 genes were further filtered by removing lowly expressed genes, yielding 11,813 sense-strand genes (48% of annotated B. taurus genes) that were considered for downstream analyses.

Analysis of Differential Gene Expression from RNA-seq Data

Following preliminary RNA-seq analysis, the sequence reads that mapped to unique locations in the B. taurus reference genome were used to generate lists of DE genes between the MAP-infected and control MDM groups at 2 and 6 hpi (the 0 h control MDM samples were not used for this phase of the analysis). Using an FDR threshold of ≤0.05, at 2 hpi 209 genes were significantly upregulated and 36 genes were significantly downregulated (Table S2 in Supplementary Material). It is important to note that the number of DE genes observed between MAP-infected and control MDM samples at 2 hpi was markedly higher for upregulated genes (209) compared to the downregulated genes (36). Inspection of the list of DE genes in Table S2 in Supplementary Material at the 2 hpi time point reveals that many of the top-ranked DE genes by FDR-adjusted P value have immune-related functions; for example, the v-maf avian musculoaponeurotic fibrosarcoma oncogene homolog F gene (MAFF); the nuclear factor of kappa light polypeptide gene enhancer in B-cells inhibitor, delta gene (NFKBID), the chemokine (C–C motif) ligand 3 gene (CCL3), and the chemokine (C–C motif) ligand 4 gene (CCL4). Table 1 shows the top 10 upregulated and top 10 downregulated DE genes between MAP-infected and control MDM samples at 2 hpi ranked by fold-change and with FDR-adjusted P values ≤0.05.

TABLE 1

Table 1. The top 10 upregulated and downregulated DE genes (FDR ≤ 0.05) for MAP-infected versus control MDM samples at 2 hpi as ranked by fold-change.

Notably, the difference between the numbers of upregulated and downregulated DE genes (FDR ≤ 0.05) was not as marked at the 6 hpi time point (342 upregulated versus 232 downregulated genes). These results are in broad agreement with the previous microarray analysis (590 upregulated genes and 384 downregulated genes with an FDR-adjusted P value ≤0.10) (30). Top ranking genes by FDR-adjusted P value for 6 hpi (Table S2 in Supplementary Material) included the mucosal vascular address in cell adhesion molecule 1 gene (MADCAM1), the family with sequence similarity 129, member A gene (FAM129A), the CD40 molecule, TNF receptor superfamily member 5 gene (CD40), and the phospholipid transfer protein gene (PLTP). Table 2 shows the top 10 upregulated and top 10 downregulated DE genes between MAP-infected and control MDM samples at 6 hpi ranked by fold-change and with FDR-adjusted P values ≤0.05.

TABLE 2

Table 2. The top 10 upregulated and downregulated DE genes (FDR ≤ 0.05) for MAP-infected versus control MDM samples at 6 hpi as ranked by fold-change.

The DE genes were compared according to direction of expression between 2 and 6 hpi. As shown in Figure 1, 59 genes (54 upregulated and 5 downregulated) were DE at both time points while also displaying the same direction of expression. By comparison, 186 genes (155 upregulated and 31 downregulated) and 515 genes (288 upregulated and 227 downregulated) were observed to be uniquely DE at 2 and 6 hpi, respectively. The relatively low overlap of DE genes between the two time points most likely represents evolution of the MDM transcriptional response to MAP infection over the time course.

FIGURE 1

Figure 1. A Venn diagram showing the numbers of DE genes identified at 2 and 6 hpi. Overlap comparison of DE genes detected in MAP-infected MDM versus control non-infected MDM between 2 and 6 hpi using the RNA-seq dataset. Sets of upregulated genes are represented in red and sets of downregulated genes are shown in green.

Functional Categorization of DE Genes Detected with RNA-seq

Functional categorization of DE genes was performed using the Bioconductor GOseq package (41) at 2 and 6 hpi time points to identify enriched Biological Process GO functions. At 2 hpi, we identified 149 significantly overrepresented Biological Processes (Bonferroni-adjusted P value ≤0.05) (Table S3 in Supplementary Material). Among the top ranked (based on Bonferroni-adjusted P values) Biological Processes were inflammatory response, defense response, response to stimulus, response to stress, immune system process, signaling, and signal transduction (Figure 2A). In addition, at 6 hpi, there were 40 significantly over-represented Biological Processes (Bonferroni-adjusted P value ≤0.05) (Table S4 in Supplementary Material), including immune system process, regulation of signaling, regulation of cell communication, immune response, cell communication, regulation of response to stimulus, signaling, and defense response (Figure 2B). The significantly overrepresented Biological Processes are relatively similar between the two post-infection time points and are, for the most part, associated with immunobiology.

FIGURE 2

Figure 2. The 15 top-ranked overrepresented biological process GO functions identified using the GOseq package. (A) Pie chart of the enriched biological processes generated from DE genes at 2 hpi using the RNA-seq dataset. (B) Pie chart of the enriched biological processes generated from DE genes at 6 hpi using the RNA-seq dataset. The values below each function represent the ratio of DE genes versus the total gene set for each functional category and the Bonferroni-adjusted P value.

IPA was used to identify the canonical pathways that were enriched for DE genes at both post-infection time points. In the current study, we identified 155 and 177 canonical pathways that were significantly enriched (FDR-adjusted P value ≤0.05) at 2 and 6 hpi, respectively. It is notable that all of the top 10 ranking canonical pathways identified at 2 hpi have immunobiological functions (Table S5 in Supplementary Material). These canonical pathways include IL-10 signaling, the first ranked pathway, which is shown overlaid with gene expression results in Figure 3 and CD40 signaling, the fourth ranking pathway, which is presented in Figure 4. The top ranking canonical pathways at 6 hpi (Table S6 in Supplementary Material) included Interferon signaling (second ranked pathway), IL-15 signaling (third ranked pathway), and P13K signaling in B lymphocytes (fourth ranked pathway).

FIGURE 3

Figure 3. The top-ranked enriched canonical pathway identified using IPA at 2 hpi – the IL-10 signaling pathway. Red shading indicates increased expression in MAP-infected MDM relative to the non-infected control MDM. Green shading indicates decreased expression in MAP-infected MDM relative to the non-infected control MDM. White and gray shading indicates non-expression and non-differential expression, respectively.

FIGURE 4

Figure 4. The fourth-ranked enriched canonical pathway identified using IPA at 2 hpi – the CD40 signaling pathway. Red shading indicates increased expression in MAP-infected MDM relative to the non-infected control MDM. Green shading indicates decreased expression in MAP-infected MDM relative to the non-infected control MDM. White and gray shading indicates non-expression and non-differential expression, respectively.

Comparative Analyses of DE Genes Detected Using RNA-seq and Microarray Technologies

Total RNA samples purified from the MAP-infected and control non-infected MDM were analyzed previously by us using the Affymetrix^® Bovine Genome Array (30). To directly compare the gene expression results between the RNA-seq and microarray platforms, we re-analyzed the microarray data for the 0 h, 2, and 6 hpi time points (35 samples).

Of the 24,072 probe sets represented on the array, 11,259 probe sets were informative that represented 5,542 unique genes with Ensembl bovine gene ID. Prior to differential gene expression analysis, the data from the 11,259 informative probes was used to generate multi-dimensional scaling (MDS) plots at 2 and 6 hpi. Using the same procedure, MDS plots were also produced from the equivalent RNA-seq data at 2 and 6 hpi using all detectable genes (11,813 genes) (Figure S2 in Supplementary Material). Inspection of these MDS plots shows relatively small separation of individual samples according to their infection status by MDS dimension axis at either time point post-infection for the two technologies. This feature of both gene expression data sets may be due to the signal from DE genes being obscured by the background gene expression noise of the majority of detectable genes.

Further analysis of the microarray data showed that 315 genes (201 upregulated and 114 downregulated) were significantly DE at 2 hpi (FDR-adjusted P value ≤0.05). Comparison of overlapping DE genes between the microarray and RNA-seq data sets revealed 134 DE genes that displayed the same direction of expression with both technologies and 292 genes that were only DE on a single platform (181 DE genes unique to the microarray and 111 DE genes unique to RNA-seq) (Figure S3A in Supplementary Material).

At 6 hpi, 466 genes (307 upregulated and 159 downregulated) were DE in the MAP-infected relative to the control MDM based on the microarray data (FDR-adjusted P value ≤0.05). Comparison of common DE genes across the microarray and RNA-seq platforms revealed 189 DE genes displaying the same direction of expression for the two technologies. The remaining 662 genes were detected as DE using a single platform (277 DE genes unique to the microarray and 385 DE genes unique to RNA-seq) (Figure S3B in Supplementary Material).

Detailed information for all DE genes detected using the Affymetrix^® microarray in MAP-infected versus control non-infected MDM samples at 2 and 6 hpi is provided in Table S7 in Supplementary Material.

Estimation of Dynamic Range from RNA-seq and Microarray Data

To estimate the dynamic range of the RNA-seq and microarray platforms, the log₂ reads per kilobase per million mapped reads (RPKM) from the RNA-seq data and the log₂ intensities from the microarray data were analyzed as described by Nalpas et al. (29). The lowest gene expression value was subtracted from the highest gene expression value for each platform. For the RNA-seq platform, a log₂ dynamic range of 25.31 was estimated based on the FAT3 gene (ENSBTAG00000004081, log₂ RPKM = −9.15) and the FTH1 gene (ENSBTAG00000011184, log₂ RPKM = 16.16). For the microarray platform, this calculation yielded an estimated log₂ dynamic range of 13.56 based on the ZCCHC8 gene (ENSBTAG00000006114, log₂ intensity = 2.03) and the B2M gene (ENSBTAG00000012330, log₂ intensity = 15.59). These observations demonstrate for the present study that the dynamic range of the RNA-seq technology was 3,444-fold greater than that of the microarray platform.

Correlation of Observed log₂ Expression Values and log₂ Fold-Changes between the RNA-seq and Microarray Platforms

We next examined the correlation between the log₂ expression values generated using the RNA-seq (log₂ RPKM values) and microarray (log₂ intensity values) platforms for all genes that passed the filtering criteria and for which a definite gene length could be determined (RPKM values cannot be computed for genes with splicing events). Spearman rank correlation coefficient (ρ) values for the 4,844 filtered genes (common to both platforms) were then calculated separately for the MAP-infected and control groups at each post-infection time point. At 2 hpi, highly significant ρ values of 0.68 (P < 1.0 × 10⁻¹⁵) and 0.67 (P < 1.0 × 10⁻¹⁵) were observed for the MAP-infected and control sample groups, respectively. Similarly, at 6 hpi, highly significant ρ values were also observed: 0.68 (P < 1.0 × 10⁻¹⁵) for the MAP-infected sample group and 0.66 (P < 1.0 × 10⁻¹⁵) for the control sample group.

Following this, log₂ expression fold-changes were examined directly for the 5,419 genes that overlapped the RNA-seq and the microarray platform at both post-infection time points (this included the 4,844 gene transcripts detailed above, plus the 575 additional overlapping genes that exhibited alternative transcripts). Again, highly significant ρ values were observed for the correlation between log₂ expression fold-change values for RNA-seq and the microarray platform at both 2 hpi (0.46; P < 1.0 × 10⁻¹⁵) and 6 hpi (0.52; P < 1.0 × 10⁻¹⁵).

The results of these analyses, using both log₂ expression values and log₂ expression fold-changes from the two post-infection time points, support the reproducibility and robustness of gene expression studies on the same samples using RNA-seq and the Affymetrix^® microarray platform.

Discussion

In recent years, high-throughput functional genomics and systems analysis of the mammalian host response to a range of mycobacterial pathogens has greatly enriched scientific understanding of the immunobiology of these infections (47–51). In particular, transcriptomics and downstream systems biology analyses of the in vitro macrophage response to mycobacteria have been particularly informative regarding host–pathogen interactions that underlie pathogenesis and which are reflected in perturbation of host genes and cellular pathways (25, 29, 30, 32, 52–57).

Most of this research work has been performed using various types of microarray platform, which until recently, has been the technology of choice for transcriptomics studies of the host response to infection. However, during the last 6 years, RNA-seq has emerged as the most powerful tool for high-resolution interrogation of the eukaryotic transcriptome in response to external stimuli such as invasive pathogens (58–61). RNA-seq has significant advantages over microarrays for surveying the complex transcriptional landscape of multicellular organisms. For example, microarray construction and implementation requires pre-existing genome sequence information for probe design, while pre-selection of the genes and transcripts to be interrogated by the microarray may result in a biased representation of the transcriptome. In contrast, RNA-seq offers unbiased, genome-wide transcriptome profiling of host gene expression without the requirement of pre-existing genome sequence information prior to the initiation of experiments. In addition, compared to microarrays, which have a dynamic range constrained by technical factors (for example, probe saturation for highly expressed genes, or lack of detectable probe hybridization signal for lowly expressed genes), the dynamic range of RNA-seq is normally limited only by the depth of sequencing used for a particular experimental comparison, thereby leading to higher sensitivity for detection of lowly expressed transcripts. Also, where appropriate, RNA-seq data can be used to quantify alternatively spliced gene variants; identify novel transcribed genes; and study antisense transcription (28, 29, 62, 63). Consequently, for the present study, an RNA-seq approach was used to study the bovine host macrophage response to MAP infection in vitro across an experimental time course consisting of 2 and 6 hpi time points.

Differential Gene Expression and Functional Biology of RNA-seq Results

RNA-seq analysis demonstrated that of the 245 significantly DE genes detected at 2 hpi, 85.3% of these were upregulated in MAP-infected MDM compared to non-infected control MDM. Also at 6 hpi, 59.6% of the DE genes were upregulated in infected MDM relative to the controls (Table S2 in Supplementary Material). This pattern of a general increase in gene expression in bovine macrophages within the first 6 h of MAP infection in vitro has also been observed by our group and other workers (30, 64, 65). It is also noteworthy that the mean absolute log₂ fold-change in expression for upregulated genes in MAP-infected MDM was markedly higher than for downregulated genes at both 2 hpi (1.95 versus 1.07, respectively) and 6 hpi (1.50 versus 0.92, respectively). This is consistent with results obtained by Magee and colleagues using bovine MDM infected with M. bovis (32).

The most upregulated DE gene (ranked by fold-change) observed at 2 hpi from RNA-seq was the CSF3 gene (log₂ fold-change = +8.05, Table 1), which encodes a cytokine that controls the production, differentiation, and function of granulocytes and which has also been shown to be highly upregulated in MAP-infected MDM isolated from red deer (Cervus elaphus) (66). It is interesting to note that Marfell and colleagues also observed that upregulation of this gene was higher in susceptible animals compared to resistant animals. The most downregulated annotated gene at 2 hpi using RNA-seq was the RAB3A gene (log₂ fold-change = −2.10, Table 1), which plays an important role in intracellular vesicle and membrane trafficking (67). While this gene has not previously been shown to be associated with macrophage–mycobacteria interactions, its downregulation could reflect an aspect of the disruption of phagosome–lysosome fusion mediated by MAP to promote its survival (68).

The most upregulated DE gene (ranked by fold-change) detected at 6 hpi using RNA-seq was the LOXL4 gene (log₂ fold-change = +5.42, Table 2), which has not previously been associated with a functional role in macrophage–mycobacteria interactions, but has a primary role in connective tissue biogenesis (69). However, recent findings have suggested that the LOX family of proteins may also have an ancillary transcriptional regulatory function (70). The most downregulated gene at 6 hpi detected using RNA-seq was the OPRD1 gene (log₂ fold-change = −3.21, Table 2), which encodes an opioid receptor also not previously reported to be involved in the macrophage response to intracellular pathogens. However, it has been demonstrated that TNF-α and IL-1β can downregulate the expression of opioid receptors at the mRNA level (71).

The identification of DE genes that hitherto had no documented role in macrophage–mycobacterial interactions highlights the potential of RNA-seq for revealing novel layers of information regarding host cellular processes induced following MAP infection and the roles that these genes may play in the host immune responses to MAP infection.

Several pro-inflammatory cytokine and chemokine genes, including CCL20, CXCL2, CXCL3, IL1B, and TNF, were DE at the 2 hpi time point; previous studies have highlighted the important role played by the products of these genes in regulating the innate immune response to mycobacterial infection (15, 20, 23, 24). The pro-inflammatory response to infection is further supported by the perturbation of several immunological signaling pathways including CD40 signaling (Figure 4), which is required for activation of antigen-presenting cells such as the macrophage (72, 73); IL-15 signaling, which regulates pro-inflammatory cytokine production in the macrophage (74); and interferon signaling (75, 76).

Furthermore, IL-1 pro-inflammatory cytokine expression in the MAP-infected host is critical not only to protective immunity but also to MAP survival. IL-1 cytokines are key effector cytokines produced by macrophages in response to infection with MAP. Indeed, IL-1 cytokine expression was detected as early as 10 min after infection with MAP under experimental infection conditions and interestingly, co-culture systems have shown that the macrophages recruited as a result of epithelial cell-induced IL-1 cytokines can be exploited by MAP to enhance their survival within the host (77). It is noteworthy that in the present study, both IL1A and IL1B are significantly DE in MAP-infected MDM at 2 hpi (IL1A log₂ fold-change = +4.9; IL1B log₂ fold-change = +5.6).

In order to produce the mature forms of IL-1 cytokines, the inflammasome is required. This pro-inflammatory protein complex occurs in myeloid cells upon infection to coordinate the activation of effective anti-bacterial innate immunity (78). The exact composition of the inflammasome varies depending on the activator (e.g., bacterial toxin, bacterial components, flagellin, and dsDNA); however, it has not been well defined in bovine studies (79). Both NLRP3 (log₂ fold-change at 2 hpi = +2.7) and IRAK2 (log₂ fold-change at 2 hpi = +2.2) are important components of the NLRP3-inflammasome complex. Indeed, Nlrp3^−/− knockout mice do not produce IL-1 cytokines (80).

Other genes encoding proteins associated with induction and activation of the inflammasome include SAA3 (2 and 6 hpi) (81) – which encodes an important acute phase protein of macrophages – and CASP4 (6 hpi) (82) – which encodes a protease with a well-characterized role in programed cell death. In contrast, CASP8, which also exhibited increased expression at 6 hpi, encodes caspase 8, which has an inflammasome-blocking function (83). Therefore, CASP8 upregulation may reflect host-directed control of inflammasome activation or, possibly, immunoevasive modulation by mycobacterial factors. Previous work has demonstrated that mycobacteria, such as M. tuberculosis, can block inflammasome activation as a novel immune evasion strategy (79, 84). In addition, lung infection with M. tuberculosis generates increased NO expression levels, which negatively regulates the NLRP3-inflammasome, thereby decreasing IL-1β production (85). Therefore, the results for MAP infection of bovine MDM may be a useful avenue for future studies regarding the interplay between bovine macrophages and MAP.

The genes encoding IL-1RN, and the anti-inflammatory cytokine IL-10 – two important regulators of IL-1 cytokine family activity – are both DE [IL1RN log₂ fold-change = +1.1 (2 hpi), +1.9 (6 hpi); IL10 log₂ fold-change = +2.01 (2 hpi)]. Notably, IL-10 signaling is also the top ranked canonical pathway identified by IPA at 2 hpi (Figure 3). IL10 encodes an immunosuppressive cytokine that regulates the antimicrobial activity of the macrophage, thus limiting the level of cytokine-induced tissue damage. Upregulation of IL-10 expression induced by mycobacteria has been proposed to inhibit host innate immune responses during infection resulting in enhanced pathogen survival (86–88).

Our findings support the hypothesis that the immunomodulatory mechanisms employed by MAP are reflected in the host macrophage transcriptome. Ultimately, protection against mycobacterial infection is a balance between protection and pathology (89). While there is significant activation of a pro- inflammatory immune response in MDM at 2 hpi, it is clear that this response is quickly regulated as the pro-inflammatory mediators are no longer DE at 6 hpi. In this regard, the outcome of infection is decided by the balance between pro- and anti-inflammatory mediators (24, 90–92).

Technical Comparison of RNA-seq and Microarray Technologies for Gene Expression Analysis

Previously, the MDM-extracted RNA samples analyzed for the present study were examined with the pan-genomic Affymetrix^® GeneChip^® Bovine Genome Array microarray platform. Here, we have used new RNA-seq data and a re-analyzed microarray data set to perform a direct technical comparison of gene expression estimation between the two platforms. The number of DE genes identified 2 hpi was higher in the microarray compared to RNA-seq (315 versus 245), while conversely, the number of DE genes 6 hpi detected via RNA-seq analysis exceeded those detected by the microarray (574 versus 466). In total, across the two infection time points the number of DE genes was higher based on the RNA-seq data compared to the microarray data (819 versus 781; an increase of 5%). Although this increase in the number of RNA-seq-identified DE genes relative to microarray-identified DE genes is lower than that previously reported by us in a technical comparison of RNA extracted from M. bovis-infected and non-infected MDM, this finding is consistent with other studies demonstrating greater numbers of DE genes identified by RNA-seq compared to microarray analysis of the same samples (29, 93, 94).

The increased number of DE genes detected by RNA-seq may be attributed to the increased dynamic range of RNA-seq relative to the microarray, which permits the detection of lowly expressed DE genes between MAP-infected and non-infected control MDM (29, 94–96). Furthermore, the concordance between the two platforms, as estimated by the percentage of DE genes common to both platforms across both infection time points, was 41.36% (323/781 genes) for the microarray and 39.44% (323/819 genes) for RNA-seq. These estimates are in agreement with the concordance previously determined for a comparison of bovine MDM infected with M. bovis and control non-infected MDM (29). The differences observed in gene expression estimation between the two platforms may be explained by several technical and analytical factors including: (1) systematic differences in dynamic range; (2) differences in the statistical models used to analyze digital/count gene expression data such as that generated using RNA-seq and the analog/continuous data obtained from microarrays; and (3) differences in the mRNA transcripts analyzed by both platforms (for example, the probes on Affymetrix GeneChip arrays are predominantly based on sequences at the 3′ end of genes and are therefore 3′ biased, while RNA-seq read data are expected to be more equally distributed across gene transcripts) (38, 46, 97–101).

In summary, the present study describes a transcriptomics survey of the host macrophage response to MAP infection using bovine MDM as an experimental model. We have used RNA-seq data generated from MDM infected with a clinical strain of MAP across a 6 h infection time course and compared the results of the RNA-seq analysis to a comparable re-analysis of microarray data obtained using the same experimental samples. The results of this work provide new insights into macrophage-MAP interplay, highlighting potential functional roles for genes that previously have not been implicated in the host response to infection with MAP bacilli. Furthermore, the pro- and anti-inflammatory cytokines involved in the bovine MDM response to MAP infection, such as those associated with the IL-10 signaling pathway, emphasize the balance between protective host immunity and bacilli survival and proliferation. Finally, by directly comparing the performance of two transcriptomics platforms, we demonstrate that RNA-seq represents a superior technology to microarrays for in vitro analyses of gene expression using mammalian cells infected with intracellular bacterial pathogens.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

This work was supported by a Teagasc Walsh Fellowship to Maura E. Casey (No.: RMIS 6082), Investigator Grants from Science Foundation Ireland (No.: SFI/01/F.1/B028 and SFI/08/IN.1/B2038), a Research Stimulus Grant from the Department of Agriculture, Fisheries and the Marine (No: RSF 06 405), and a European Union Framework 7 Project Grant (No: KBBE-211602-MACROSYS). Kate E. Killick was supported by the Irish Research Council (IRC) funded Bioinformatics and Systems Biology thematic Ph.D. program (http://bioinfo-casl.ucd.ie/PhD). We thank Mr. Kevin Kenny, Mr. Eamon Costello, and staff at the Central Veterinary Research Laboratory (CVRL), Backweston, Co. Kildare for culturing and providing the field isolate of M. avium subsp. paratuberculosis. We also would like to thank the staff at the UCD Lyons Research Farm for assistance with blood collection.

Supplementary Material

The Supplementary Material for this article can be found online at http://www.frontiersin.org/Journal/10.3389/fimmu.2015.00023/abstract

Footnotes

1. ^www.bioinformatics.babraham.ac.uk/projects/fastqc
2. ^http://www.ingenuity.com

References

1. Greig A, Stevenson K, Perez V, Pirie AA, Grant JM, Sharp JM. Paratuberculosis in wild rabbits (Oryctolagus cuniculus). Vet Rec (1997) 140:141–3. doi:10.1136/vr.140.6.141

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

2. Over K, Crandall PG, O’Bryan CA, Ricke SC. Current perspectives on Mycobacterium avium subsp. paratuberculosis, Johne’s disease, and Crohn’s disease: a review. Crit Rev Microbiol (2011) 37:141–56. doi:10.3109/1040841x.2010.532480

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

3. Davis WC, Madsen-Bouterse SA. Crohn’s disease and Mycobacterium avium subsp. paratuberculosis: the need for a study is long overdue. Vet Immunol Immunopathol (2012) 145:1–6. doi:10.1016/j.vetimm.2011.12.005

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

4. Robertson R, Hill B, Cerf O, Jordan K, Venter P. A commentary on current perspectives on Mycobacterium avium subsp. paratuberculosis, Johne’s disease and Crohn’s disease: a review by Over et al. (2011). Crit Rev Microbiol (2012) 38:183–4. doi:10.3109/1040841x.2011.649712

CrossRef Full Text | Google Scholar

5. Boelaert F, Walravens K, Biront P, Vermeersch JP, Berkvens D, Godfroid J. Prevalence of paratuberculosis (Johne’s disease) in the Belgian cattle population. Vet Microbiol (2000) 77:269–81. doi:10.1016/S0378-1135(00)00312-6

CrossRef Full Text | Google Scholar

6. Muskens J, Barkema H, Russchen E. Prevalence and regional distribution of paratuberculosis in dairy herds in the Netherlands. Vet Microbiol (2000) 77:253–61. doi:10.1016/S0378-1135(00)00310-2

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

7. Nielsen S, Thamsborg S, Houe H. Bulk-tank milk ELISA antibodies for estimating the prevalence of paratuberculosis in Danish dairy herds. Prev Vet Med (2000) 44:1–7. doi:10.1016/S0167-5877(00)00098-2

CrossRef Full Text | Google Scholar

8. Good M, Clegg T, Sheridan H, Yearsely D, O’Brien T, Egan J, et al. Prevalence and distribution of paratuberculosis (Johne’s disease) in cattle herds in Ireland. Ir Vet J (2009) 62:597–606. doi:10.1186/2046-0481-62-9-597

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

9. Cho J, Tauer LW, Schukken YH, Gomez MI, Smith RL, Lu Z, et al. Economic analysis of Mycobacterium avium subspecies paratuberculosis vaccines in dairy herds. J Dairy Sci (2012) 95:1855–72. doi:10.3168/jds.2011-4787

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

10. Coussens PM, Sipkovsky S, Murphy B, Roussey J, Colvin CJ. Regulatory T cells in cattle and their potential role in bovine paratuberculosis. Comp Immunol Microbiol Infect Dis (2012) 35:233–9. doi:10.1016/j.cimid.2012.01.004

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

11. Sweeney RW. Transmission of paratuberculosis. Vet Clin North Am Food Anim Pract (1996) 12:305–12.

Google Scholar

12. Lombard JE. Epidemiology and economics of paratuberculosis. Vet Clin North Am Food Anim Pract (2011) 27:525–535,v. doi:10.1016/j.cvfa.2011.07.012

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

13. Stabel JR. Transitions in immune responses to Mycobacterium paratuberculosis. Vet Microbiol (2000) 77:465–73. doi:10.1016/S0378-1135(00)00331-X

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

14. Tessema M, Koets A, Rutten V, Gruys E. How does Mycobacterium avium subsp. paratuberculosis resist intracellular degradation? Vet Q (2001) 23:153–62. doi:10.1080/01652176.2001.9695105

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

15. Bannantine JP, Bermudez LE. No holes barred: invasion of the intestinal mucosa by Mycobacterium avium subspecies paratuberculosis. Infect Immun (2013) 81:3960–5. doi:10.1128/IAI.00575-13

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

16. Delbridge LM, O’Riordan MX. Innate recognition of intracellular bacteria. Curr Opin Immunol (2007) 19:10–6. doi:10.1016/j.coi.2006.11.005

CrossRef Full Text | Google Scholar

17. Killick KE, Ni Cheallaigh C, O’Farrelly C, Hokamp K, MacHugh DE, Harris J. Receptor-mediated recognition of mycobacterial pathogens. Cell Microbiol (2013) 15:1484–95. doi:10.1111/cmi.12161

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

18. Ferwerda G, Kullberg BJ, de Jong DJ, Girardin SE, Langenberg DM, van Crevel R, et al. Mycobacterium paratuberculosis is recognized by Toll-like receptors and NOD2. J Leukoc Biol (2007) 82:1011–8. doi:10.1189/jlb.0307147

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

19. Janagama HK, Jeong K, Kapur V, Coussens P, Sreevatsan S. Cytokine responses of bovine macrophages to diverse clinical Mycobacterium avium subspecies paratuberculosis strains. BMC Microbiol (2006) 6:10. doi:10.1186/1471-2180-6-10

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

20. Stabel JR. Host responses to Mycobacterium avium subsp. paratuberculosis: a complex arsenal. Anim Health Res Rev (2006) 7:61–70. doi:10.1017/s1466252307001168

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

21. Sweeney RW. Pathogenesis of paratuberculosis. Vet Clin North Am Food Anim Pract (2011) 27:537–546,v. doi:10.1016/j.cvfa.2011.07.001

CrossRef Full Text | Google Scholar

22. Olsen I, Sigurgardottir G, Djonne B. Paratuberculosis with special reference to cattle. A review. Vet Q (2002) 24:12–28. doi:10.1080/01652176.2002.9695120

CrossRef Full Text | Google Scholar

23. Coussens PM. Model for immune responses to Mycobacterium avium subspecies paratuberculosis in cattle. Infect Immun (2004) 72:3089–96. doi:10.1128/iai.72.6.3089-3096.2004

CrossRef Full Text | Google Scholar

24. Weiss DJ, Souza CD. Review paper: modulation of mononuclear phagocyte function by Mycobacterium avium subsp. paratuberculosis. Vet Pathol (2008) 45:829–41. doi:10.1354/vp.45-6-829

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

25. Kabara E, Coussens PM. Infection of primary bovine macrophages with Mycobacterium avium subspecies paratuberculosis suppresses host cell apoptosis. Front Microbiol (2012) 3:215. doi:10.3389/fmicb.2012.00215

CrossRef Full Text | Google Scholar

26. Souza C, Davis WC, Eckstein TM, Sreevatsan S, Weiss DJ. Mannosylated lipoarabinomannans from Mycobacterium avium subsp. paratuberculosis alters the inflammatory response by bovine macrophages and suppresses killing of Mycobacterium avium subsp. avium organisms. PLoS One (2013) 8:e75924. doi:10.1371/journal.pone.0075924

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

27. Dore E, Pare J, Cote G, Buczinski S, Labrecque O, Roy JP, et al. Risk factors associated with transmission of Mycobacterium avium subsp. paratuberculosis to calves within dairy herd: a systematic review. J Vet Intern Med (2012) 26:32–45. doi:10.1111/j.1939-1676.2011.00854.x

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

28. McGettigan PA. Transcriptomics in the RNA-seq era. Curr Opin Chem Biol (2013) 17:4–11. doi:10.1016/j.cbpa.2012.12.008

CrossRef Full Text | Google Scholar

29. Nalpas NC, Park SDE, Magee DA, Taraktsoglou M, Browne JA, Conlon KM, et al. Whole-transcriptome, high-throughput RNA sequence analysis of the bovine macrophage response to Mycobacterium bovis infection in vitro. BMC Genomics (2013) 14:230. doi:10.1186/1471-2164-14-230

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

30. MacHugh D, Taraktsoglou M, Killick K, Nalpas N, Browne J, Park SDE, et al. Pan-genomic analysis of bovine monocyte-derived macrophage gene expression in response to in vitro infection with Mycobacterium avium subspecies paratuberculosis. Vet Res (2012) 43:25. doi:10.1186/1297-9716-43-25

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

31. Taraktsoglou M, Szalabska U, Magee D, Browne J, Sweeney T, Gormley E, et al. Transcriptional profiling of immune genes in bovine monocyte-derived macrophages exposed to bacterial antigens. Vet Immunol Immunopathol (2011) 140:130–9. doi:10.1016/j.vetimm.2010.12.002

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

32. Magee DA, Taraktsoglou M, Killick KE, Nalpas NC, Browne JA, Park SDE, et al. Global gene expression and systems biology analysis of bovine monocyte-derived macrophages in response to in vitro challenge with Mycobacterium bovis. PLoS One (2012) 7:e32034. doi:10.1371/journal.pone.0032034

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

33. McLoughlin KE, Nalpas NC, Rue-Albrecht K, Browne JA, Magee DA, Killick KE, et al. RNA-seq transcriptional profiling of peripheral blood leukocytes from cattle infected with Mycobacterium bovis. Front Immunol (2014) 5:396. doi:10.3389/fimmu.2014.00396

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

34. Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics (2013) 29:15–21. doi:10.1093/bioinformatics/bts635

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

35. Zimin AV, Delcher AL, Florea L, Kelley DR, Schatz MC, Puiu D, et al. A whole-genome assembly of the domestic cow, Bos taurus. Genome Biol (2009) 10:R42. doi:10.1186/gb-2009-10-4-r42

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

36. Liao Y, Smyth GK, Shi W. The subread aligner: fast, accurate and scalable read mapping by seed-and-vote. Nucleic Acids Res (2013) 41:e108. doi:10.1093/nar/gkt214

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

37. Liao Y, Smyth GK, Shi W. featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics (2014) 30:923–30. doi:10.1093/bioinformatics/btt656

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

38. Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics (2010) 26:139–40. doi:10.1093/bioinformatics/btp616

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

39. Durinck S, Spellman PT, Birney E, Huber W. Mapping identifiers for the integration of genomic datasets with the R/bioconductor package biomaRt. Nat Protoc (2009) 4:1184–91. doi:10.1038/nprot.2009.97

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

40. Benjamini Y, Hochberg Y. Controlling the false discovery rate – a practical and powerful approach to multiple testing. J R Stat Soc Ser B Methodol (1995) 57:289–300.

Google Scholar

41. Young MD, Wakefield MJ, Smyth GK, Oshlack A. Gene ontology analysis for RNA-seq: accounting for selection bias. Genome Biol (2010) 11:R14. doi:10.1186/gb-2010-11-2-r14

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

42. Barrett T, Wilhite SE, Ledoux P, Evangelista C, Kim IF, Tomashevsky M, et al. NCBI GEO: archive for functional genomics data sets – update. Nucleic Acids Res (2013) 41:D991–5. doi:10.1093/nar/gks1193

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

43. Gentleman RC, Carey VJ, Bates DM, Bolstad B, Dettling M, Dudoit S, et al. Bioconductor: open software development for computational biology and bioinformatics. Genome Biol (2004) 5:R80. doi:10.1186/gb-2004-5-10-r80

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

44. Hochreiter S, Clevert DA, Obermayer K. A new summarization method for Affymetrix probe level data. Bioinformatics (2006) 22:943–9. doi:10.1093/bioinformatics/btl033

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

45. Talloen W, Clevert DA, Hochreiter S, Amaratunga D, Bijnens L, Kass S, et al. I/NI-calls for the exclusion of non-informative genes: a highly effective filtering tool for microarray data. Bioinformatics (2007) 23:2897–902. doi:10.1093/bioinformatics/btm478

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

46. Smyth GK. Linear models and empirical bayes methods for assessing differential expression in microarray experiments. Stat Appl Genet Mol Biol (2004) 3:Article3. doi:10.2202/1544-6115.1027

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

47. Rachman H, Kaufmann SH. Exploring functional genomics for the development of novel intervention strategies against tuberculosis. Int J Med Microbiol (2007) 297:559–67. doi:10.1016/j.ijmm.2007.03.003

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

48. MacHugh DE, Gormley E, Park SDE, Browne JA, Taraktsoglou M, O’Farrelly C, et al. Gene expression profiling of the host response to Mycobacterium bovis infection in cattle. Transbound Emerg Dis (2009) 56:204–14. doi:10.1111/j.1865-1682.2009.01082.x

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

49. Bannantine JP, Talaat AM. Genomic and transcriptomic studies in Mycobacterium avium subspecies paratuberculosis. Vet Immunol Immunopathol (2010) 138:303–11. doi:10.1016/j.vetimm.2010.10.008

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

50. Maertzdorf J, Weiner J III, Kaufmann SH. Enabling biomarkers for tuberculosis control. Int J Tuberc Lung Dis (2012) 16:1140–8. doi:10.5588/ijtld.12.0246

CrossRef Full Text | Google Scholar

51. Berry MP, Blankley S, Graham CM, Bloom CI, O’Garra A. Systems approaches to studying the immune response in tuberculosis. Curr Opin Immunol (2013) 25:579–87. doi:10.1016/j.coi.2013.08.003

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

52. Rohde KH, Veiga DF, Caldwell S, Balazsi G, Russell DG. Linking the transcriptional profiles and the physiological states of Mycobacterium tuberculosis during an extended intracellular infection. PLoS Pathog (2012) 8:e1002769. doi:10.1371/journal.ppat.1002769

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

53. Lamont EA, Xu WW, Sreevatsan S. Host-Mycobacterium avium subsp. paratuberculosis interactome reveals a novel iron assimilation mechanism linked to nitric oxide stress during early infection. BMC Genomics (2013) 14:694. doi:10.1186/1471-2164-14-694

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

54. Guerra-Laso JM, Raposo-Garcia S, Garcia-Garcia S, Diez-Tascon C, Rivero-Lezcano OM. Microarray analysis of Mycobacterium tuberculosis-infected monocytes reveals IL-26 as a new candidate gene for tuberculosis susceptibility. Immunology (2015) 144:291–301. doi:10.1111/imm.12371

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

55. Killick KE, Magee DA, Park SDE, Taraktsoglou M, Browne JA, Conlon KM, et al. Key hub and bottleneck genes differentiate the macrophage response to virulent and attenuated Mycobacterium bovis. Front Immunol (2014) 5:422. doi:10.3389/fimmu.2014.00422

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

56. Magee DA, Conlon KM, Nalpas NC, Browne JA, Pirson C, Healy C, et al. Innate cytokine profiling of bovine alveolar macrophages reveals commonalities and divergence in the response to Mycobacterium bovis and Mycobacterium tuberculosis infection. Tuberculosis (Edinb) (2014) 94:441–50. doi:10.1016/j.tube.2014.04.004

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

57. Rue-Albrecht K, Magee DA, Killick KE, Nalpas NC, Gordon SV, MacHugh DE. Comparative functional genomics and the bovine macrophage response to strains of the Mycobacterium genus. Front Immunol (2014) 5:536. doi:10.3389/fimmu.2014.00536

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

58. Mortazavi A, Williams BA, McCue K, Schaeffer L, Wold B. Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat Methods (2008) 5:621–8. doi:10.1038/nmeth.1226

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

59. Nagalakshmi U, Wang Z, Waern K, Shou C, Raha D, Gerstein M, et al. The transcriptional landscape of the yeast genome defined by RNA sequencing. Science (2008) 320:1344–9. doi:10.1126/science.1158441

CrossRef Full Text | Google Scholar

60. Sultan M, Schulz MH, Richard H, Magen A, Klingenhoff A, Scherf M, et al. A global view of gene activity and alternative splicing by deep sequencing of the human transcriptome. Science (2008) 321:956–60. doi:10.1126/science.1160342

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

61. Wilhelm BT, Marguerat S, Watt S, Schubert F, Wood V, Goodhead I, et al. Dynamic repertoire of a eukaryotic transcriptome surveyed at single-nucleotide resolution. Nature (2008) 453:1239–43. doi:10.1038/nature07002

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

62. Ozsolak F, Milos PM. RNA sequencing: advances, challenges and opportunities. Nat Rev Genet (2011) 12:87–98. doi:10.1038/nrg2934

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

63. Mutz KO, Heilkenbrinker A, Lonne M, Walter JG, Stahl F. Transcriptome analysis using next-generation sequencing. Curr Opin Biotechnol (2013) 24:22–30. doi:10.1016/j.copbio.2012.09.004

CrossRef Full Text | Google Scholar

64. Coussens PM, Colvin CJ, Rosa GJ, Perez Laspiur J, Elftman MD. Evidence for a novel gene expression program in peripheral blood mononuclear cells from Mycobacterium avium subsp. paratuberculosis-infected cattle. Infect Immun (2003) 71:6487–98. doi:10.1128/IAI.71.11.6487-6498.2003

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

65. Kabara E, Kloss CC, Wilson M, Tempelman RJ, Sreevatsan S, Janagama H, et al. A large-scale study of differential gene expression in monocyte-derived macrophages infected with several strains of Mycobacterium avium subspecies paratuberculosis. Brief Funct Genomics (2010) 9:220–37. doi:10.1093/bfgp/elq009

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

66. Marfell BJ, O’Brien R, Griffin JF. Global gene expression profiling of monocyte-derived macrophages from red deer (Cervus elaphus) genotypically resistant or susceptible to Mycobacterium avium subspecies paratuberculosis infection. Dev Comp Immunol (2013) 40:210–7. doi:10.1016/j.dci.2013.02.004

CrossRef Full Text | Google Scholar

67. Lang T, Jahn R. Core proteins of the secretory machinery. Handb Exp Pharmacol (2008) 184:107–27. doi:10.1007/978-3-540-74805-2_5

CrossRef Full Text | Google Scholar

68. Hostetter J, Steadham E, Haynes J, Bailey T, Cheville N. Phagosomal maturation and intracellular survival of Mycobacterium avium subspecies paratuberculosis in J774 cells. Comp Immunol Microbiol Infect Dis (2003) 26:269–83. doi:10.1016/S0147-9571(02)00070-X

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

69. Molnar J, Fong KS, He QP, Hayashi K, Kim Y, Fong SF, et al. Structural and functional diversity of lysyl oxidase and the LOX-like proteins. Biochim Biophys Acta (2003) 1647:220–4. doi:10.1016/S1570-9639(03)00053-0

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

70. Iturbide A, Garcia de Herreros A, Peiro S. A new role for LOX and LOXL2 proteins in transcription regulation. FEBS J (2014). doi:10.1111/febs.12961

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

71. Shen H, Aeschlimann A, Reisch N, Gay RE, Simmen BR, Michel BA, et al. Kappa and delta opioid receptors are expressed but down-regulated in fibroblast-like synoviocytes of patients with rheumatoid arthritis and osteoarthritis. Arthritis Rheum (2005) 52:1402–10. doi:10.1002/art.21141

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

72. Sommer S, Pudrith C, Colvin C, Coussens P. Mycobacterium avium subspecies paratuberculosis suppresses expression of IL-12p40 and iNOS genes induced by signalling through CD40 in bovine monocyte-derived macrophages. Vet Immunol Immunopathol (2009) 128:44–52. doi:10.1016/j.vetimm.2008.10.294

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

73. Suttles J, Stout RD. Macrophage CD40 signaling: a pivotal regulator of disease protection and pathogenesis. Semin Immunol (2009) 21:257–64. doi:10.1016/j.smim.2009.05.011

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

74. Alleva DG, Kaser SB, Monroy MA, Fenton MJ, Beller DI. IL-15 functions as a potent autocrine regulator of macrophage proinflammatory cytokine production: evidence for differential receptor subunit utilization associated with stimulation or inhibition. J Immunol (1997) 159:2941–51.

Pubmed Abstract | Pubmed Full Text | Google Scholar

75. Coussens P, Lamont EA, Kabara E, Sreevatsan S. Host-pathogen interactions and intracellular survival of Mycobacterium avium subsp. paratuberculosis. In: Behr MA, Collins DM, editors. Paratuberculosis: Organism, Disease, Control. CABI: Wallingford, UK (2010).

Google Scholar

76. Arsenault RJ, Li Y, Bell K, Doig K, Potter A, Griebel PJ, et al. Mycobacterium avium subsp. paratuberculosis inhibits gamma interferon-induced signaling in bovine monocytes: insights into the cellular mechanisms of Johne’s disease. Infect Immun (2012) 80:3039–48. doi:10.1128/IAI.00406-12

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

77. Lamont EA, O’Grady SM, Davis WC, Eckstein T, Sreevatsan S. Infection with Mycobacterium avium subsp. paratuberculosis results in rapid interleukin-1beta release and macrophage transepithelial migration. Infect Immun (2012) 80:3225–35. doi:10.1128/IAI.06322-11

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

78. Fitzgerald KA. NLR-containing inflammasomes: central mediators of host defense and inflammation. Eur J Immunol (2010) 40:595–8. doi:10.1002/eji.201040331

CrossRef Full Text | Google Scholar

79. Briken V, Ahlbrand SE, Shah S. Mycobacterium tuberculosis and the host cell inflammasome: a complex relationship. Front Cell Infect Microbiol (2013) 3:62. doi:10.3389/fcimb.2013.00062

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

80. Schroder K, Tschopp J. The inflammasomes. Cell (2010) 140:821–32. doi:10.1016/j.cell.2010.01.040

CrossRef Full Text | Google Scholar

81. Ather JL, Ckless K, Martin R, Foley KL, Suratt BT, Boyson JE, et al. Serum amyloid A activates the NLRP3 inflammasome and promotes Th17 allergic asthma in mice. J Immunol (2011) 187:64–73. doi:10.4049/jimmunol.1100500

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

82. Sollberger G, Strittmatter GE, Kistowska M, French LE, Beer HD. Caspase-4 is required for activation of inflammasomes. J Immunol (2012) 188:1992–2000. doi:10.4049/jimmunol.1101620

CrossRef Full Text | Google Scholar

83. Kang TB, Yang SH, Toth B, Kovalenko A, Wallach D. Caspase-8 blocks kinase RIPK3-mediated activation of the NLRP3 inflammasome. Immunity (2013) 38:27–40. doi:10.1016/j.immuni.2012.09.015

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

84. Master SS, Rampini SK, Davis AS, Keller C, Ehlers S, Springer B, et al. Mycobacterium tuberculosis prevents inflammasome activation. Cell Host Microbe (2008) 3:224–32. doi:10.1016/j.chom.2008.03.003

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

85. Mishra BB, Rathinam VA, Martens GW, Martinot AJ, Kornfeld H, Fitzgerald KA, et al. Nitric oxide controls the immunopathology of tuberculosis by inhibiting NLRP3 inflammasome-dependent processing of IL-1beta. Nat Immunol (2013) 14:52–60. doi:10.1038/ni.2474

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

86. Weiss D, Evanson O, de Souza C, Abrahamsen M. A critical role of interleukin-10 in the response of bovine macrophages to infection by Mycobacterium avium subsp paratuberculosis. Am J Vet Res (2005) 66:721–6. doi:10.2460/ajvr.2005.66.721

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

87. Redford PS, Murray PJ, O’Garra A. The role of IL-10 in immune regulation during M. tuberculosis infection. Mucosal Immunol (2011) 4:261–70. doi:10.1038/mi.2011.7

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

88. Dorhoi A, Kaufmann SH. Perspectives on host adaptation in response to Mycobacterium tuberculosis: modulation of inflammation. Semin Immunol (2014) 26:533–42. doi:10.1016/j.smim.2014.10.002

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

89. Dorhoi A, Reece ST, Kaufmann SH. For better or for worse: the immune response against Mycobacterium tuberculosis balances pathology and protection. Immunol Rev (2011) 240:235–51. doi:10.1111/j.1600-065X.2010.00994.x

CrossRef Full Text | Google Scholar

90. Huynh KK, Joshi SA, Brown EJ. A delicate dance: host response to mycobacteria. Curr Opin Immunol (2011) 23:464–72. doi:10.1016/j.coi.2011.06.002

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

91. O’Garra A, Redford PS, McNab FW, Bloom CI, Wilkinson RJ, Berry MP. The immune response in tuberculosis. Annu Rev Immunol (2013) 31:475–527. doi:10.1146/annurev-immunol-032712-095939

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

92. Rajaram MV, Ni B, Dodd CE, Schlesinger LS. Macrophage immunoregulatory pathways in tuberculosis. Semin Immunol (2014) 26:471–85. doi:10.1016/j.smim.2014.09.010

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

93. Guida A, Lindstadt C, Maguire SL, Ding C, Higgins DG, Corton NJ, et al. Using RNA-seq to determine the transcriptional landscape and the hypoxic response of the pathogenic yeast Candida parapsilosis. BMC Genomics (2011) 12:628. doi:10.1186/1471-2164-12-628

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

94. Zhao S, Fung-Leung WP, Bittner A, Ngo K, Liu X. Comparison of RNA-Seq and microarray in transcriptome profiling of activated T cells. PLoS One (2014) 9:e78644. doi:10.1371/journal.pone.0078644

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

95. Auer PL, Doerge RW. Statistical design and analysis of RNA sequencing data. Genetics (2010) 185:405–16. doi:10.1534/genetics.110.114983

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

96. Tarazona S, Garcia-Alcalde F, Dopazo J, Ferrer A, Conesa A. Differential expression in RNA-seq: a matter of depth. Genome Res (2011) 21:2213–23. doi:10.1101/gr.124321.111

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

97. Marioni JC, Mason CE, Mane SM, Stephens M, Gilad Y. RNA-seq: an assessment of technical reproducibility and comparison with gene expression arrays. Genome Res (2008) 18:1509–17. doi:10.1101/gr.079558.108

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

98. Bradford JR, Hey Y, Yates T, Li Y, Pepper SD, Miller CJ. A comparison of massively parallel nucleotide sequencing with oligonucleotide microarrays for global transcription profiling. BMC Genomics (2010) 11:282. doi:10.1186/1471-2164-11-282

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

99. Roy NC, Altermann E, Park ZA, McNabb WC. A comparison of analog and next-generation transcriptomic tools for mammalian studies. Brief Funct Genomics (2011) 10:135–50. doi:10.1093/bfgp/elr005

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

100. Fang Z, Martin J, Wang Z. Statistical methods for identifying differentially expressed genes in RNA-Seq experiments. Cell Biosci (2012) 2:26. doi:10.1186/2045-3701-2-26

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar

101. Fasold M, Binder H. Estimating RNA-quality using GeneChip microarrays. BMC Genomics (2012) 13:186. doi:10.1186/1471-2164-13-186

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text | Google Scholar