# DEVELOPMENTS IN BOVINE IMMUNOLOGY – AN INTEGRATED VIEW

EDITED BY: Kieran G. Meade PUBLISHED IN: Frontiers in Immunology

#### *Frontiers Copyright Statement*

*© Copyright 2007-2015 Frontiers Media SA. All rights reserved. All content included on this site, such as text, graphics, logos, button icons, images, video/audio clips, downloads, data compilations and software, is the property of or is licensed to Frontiers Media SA ("Frontiers") or its licensees and/or subcontractors. The copyright in the text of individual articles is the property of their respective authors, subject to a license granted to Frontiers.*

*The compilation of articles constituting this e-book, wherever published, as well as the compilation of all other content on this site, is the exclusive property of Frontiers. For the conditions for downloading and copying of e-books from Frontiers' website, please see the Terms for Website Use. If purchasing Frontiers e-books from other websites or sources, the conditions of the website concerned apply.*

*Images and graphics not forming part of user-contributed materials may not be downloaded or copied without permission.*

*Individual articles may be downloaded and reproduced in accordance with the principles of the CC-BY licence subject to any copyright or other notices. They may not be re-sold as an e-book.*

*As author or other contributor you grant a CC-BY licence to others to reproduce your articles, including any graphics and third-party materials supplied by you, in accordance with the Conditions for Website Use and subject to any copyright notices which you include in connection with your articles and materials.*

*All copyright, and all rights therein, are protected by national and international copyright laws. The above represents a summary only. For the full conditions see the Conditions for Authors and the Conditions for Website Use.*

ISSN 1664-8714 ISBN 978-2-88919-632-6 DOI 10.3389/978-2-88919-632-6

### About Frontiers

Frontiers is more than just an open-access publisher of scholarly articles: it is a pioneering approach to the world of academia, radically improving the way scholarly research is managed. The grand vision of Frontiers is a world where all people have an equal opportunity to seek, share and generate knowledge. Frontiers provides immediate and permanent online open access to all its publications, but this alone is not enough to realize our grand goals.

### Frontiers Journal Series

The Frontiers Journal Series is a multi-tier and interdisciplinary set of open-access, online journals, promising a paradigm shift from the current review, selection and dissemination processes in academic publishing. All Frontiers journals are driven by researchers for researchers; therefore, they constitute a service to the scholarly community. At the same time, the Frontiers Journal Series operates on a revolutionary invention, the tiered publishing system, initially addressing specific communities of scholars, and gradually climbing up to broader public understanding, thus serving the interests of the lay society, too.

#### Dedication to quality

Each Frontiers article is a landmark of the highest quality, thanks to genuinely collaborative interactions between authors and review editors, who include some of the world's best academicians. Research must be certified by peers before entering a stream of knowledge that may eventually reach the public - and shape society; therefore, Frontiers only applies the most rigorous and unbiased reviews.

Frontiers revolutionizes research publishing by freely delivering the most outstanding research, evaluated with no bias from both the academic and social point of view. By applying the most advanced information technologies, Frontiers is catapulting scholarly publishing into a new generation.

### What are Frontiers Research Topics?

Frontiers Research Topics are very popular trademarks of the Frontiers Journals Series: they are collections of at least ten articles, all centered on a particular subject. With their unique mix of varied contributions from Original Research to Review Articles, Frontiers Research Topics unify the most influential researchers, the latest key findings and historical advances in a hot research area! Find out more on how to host your own Frontiers Research Topic or contribute to one as an author by contacting the Frontiers Editorial Office: researchtopics@frontiersin.org

## **DEVELOPMENTS IN BOVINE IMMUNOLOGY – AN INTEGRATED VIEW**

Topic Editors:

**Kieran G. Meade,** Animal & Grassland Research and Innovation Centre, Ireland

Image is Property of the Topic Editor Dr Kieran Meade.

The world's population is predicted to hit 9 Billion by 2050, and with it food demand is predicted to increase substantially. The World Bank estimates that cereal and meat production needs to increase by 50% and 85% respectively between 2000 and 2030 to meet demand, putting serious pressure on the global agricultural industry. Critical to meeting this demand for food are mechanisms to reduce the incidence of animal disease. With in excess of 1.3 billion cattle globally, the total cost of infectious diseases is difficult to estimate. However in North America alone, the cost is predicted to be \$18 billion annually. Non-infectious diseases also account for another major impediment to the production capacity and welfare of animals as well as the economic sustainability of farming. However animal diseases have implications that spread far beyond the farm gate. Infectious agents can also contaminate the food chain, and potentially affect human health.

Controlling diseases, through better preventative and treatment methods requires

a detailed understanding of the immune response in livestock species. Multiple studies have identified associations between variation in immune genes and disease susceptibility, which potentially opens up new avenues to select animals with superior disease resistance. Detailed understanding of immunity in cattle is leading to the design of more effective vaccines. Furthermore, appreciation of the significant differences between rodent and human immune responses has also led to bovine models being developed for some human diseases.

The publication of the bovine genome and the advent of next-generation sequencing technologies have facilitated a massive expansion in our knowledge of the immune response in cattle. As a result there has been an explosion of exciting research findings including in metagenomics and epigenetics. Recently, there has been a welcome move to integrate our emerging understanding of the immune response with detailed studies of other important physiological processes including nutrition and reproduction. The interactions between the reproductive system, nutrition and the immune system are of particular interest, since each places significant demands on the animal at various stages through the production cycle. The interplay between these morphologically diffuse systems involves widely distributed chemical signals in response to environmental input, and each system must interact for the normal functioning of the other. A comprehensive "systems" approach is improving our understanding of normal physiological interactions between these systems and furthermore, how dysregulation can lead to disease.

The successful translation of bovine immunological research into improved treatments for animal disease requires tight interaction between diverse scientific and clinical disciplines including immunology, microbiology, endocrinology, physiology, nutrition, reproduction and clinical veterinary medicine. With so much recent progress in the field, we believe that it is valuable and well-timed to review the broad variety of the relevant studies that attempt to increase our understanding through comprehensive collaboration between these disciplines.

We are looking forward to a wide and vivid discussion of developments in bovine immunology and related issues, and we expect that our readers profoundly benefit from new exciting insights and fruitful collaborations.

**Citation:** Kieran G. Meade, eds. (2015). Developments in bovine immunology – an integrated view. Lausanne: Frontiers Media. doi: 10.3389/978-2-88919-632-6

# Table of Contents

*05 Advances in bovine immunology – new tools and new insights to tackle old foes*

Kieran G. Meade

*07 RNA-seq transcriptional profiling of peripheral blood leukocytes from cattle infected with Mycobacterium bovis*

Kirsten E. McLoughlin, Nicolas C. Nalpas, Kévin Rue-Albrecht, John A. Browne, David A. Magee, Kate E. Killick, Stephen D. E. Park, Karsten Hokamp, Kieran G. Meade, Cliona O'Farrelly, Eamonn Gormley, Stephen V. Gordon and David E. MacHugh

*19 Analysis of the bovine monocyte-derived macrophage response to Mycobacterium avium subspecies paratuberculosis infection using RNA-seq*

Maura E. Casey, Kieran G. Meade, Nicolas C. Nalpas, Maria Taraktsoglou, John A. Browne, Kate E. Killick, Stephen D. E. Park, Eamonn Gormley, Karsten Hokamp, David A. Magee and David E. MacHugh

*33 Key hub and bottleneck genes differentiate the macrophage response to virulent and attenuated Mycobacterium bovis*

Kate E. Killick, David A. Magee, Stephen D. E. Park, Maria Taraktsoglou, John A. Browne, Kevin M. Conlon, Nicolas C. Nalpas, Eamonn Gormley, Stephen V. Gordon, David E. MacHugh and Karsten Hokamp

*46 Comparative functional genomics and the bovine macrophage response to strains of the Mycobacterium genus*

Kévin Rue-Albrecht, David A. Magee, Kate E. Killick, Nicolas C. Nalpas, Stephen V. Gordon and David E. MacHugh

*60 The single intradermal cervical comparative test interferes with Johne's disease ELISA diagnostics*

Aideen E. Kennedy, Ana T. Da Silva, Noel Byrne, Rodney Govender, John MacSharry, Jim O'Mahony and Riona G. Sayers

*68 Regulation of toll-like receptors-mediated inflammation by immunobiotics in bovine intestinal epitheliocytes: role of signaling pathways and negative regulators*

Julio Villena, Hisashi Aso and Haruki Kitazawa


Trudee Fair

*93 Bovine mastitis: frontiers in immunogenetics*

Kathleen Thompson-Crispi, Heba Atalla, Filippo Miglior and Bonnie A. Mallard

## Advances in bovine immunology – new tools and new insights to tackle old foes

#### **Kieran G. Meade\***

Animal and Bioscience Research Department, Animal and Grassland Research and Innovation Centre, Teagasc, Trim, Ireland \*Correspondence: kieran.meade@teagasc.ie

**Edited and reviewed by:** Claudia Kemper, King's College London, UK

**Keywords: mastits, M. bovis, MAP, one health, immunobiotics, bovine immunology, pregnancy, miRNAs, tuberculosis**

The United Nations predicts that the current world population of 7.2 billion is projected to reach 9.6 billion by 2050, which can only mean one thing for global agriculture, and that is increased pressure on production. One step to increase food supply is the abolition of the milk quota restrictions in the EU in 2015. It is inevitable with expansion and intensification of production that risks associated with disease will be exacerbated. This will have critical consequences for the sustainability of farming systems, for the safety of the food chain, and for human health. In this context, it is appropriate that we redouble every effort to understand the immune response, particularly to recalcitrant infectious diseases in cattle.

The publication of the bovine genome in 2009 and the advent of new high-throughput technologies have facilitated a massive expansion in our knowledge of the immune response in cattle. Genomic selection means we can now identify animals with superior genetics at birth and use them as parents of the next generation, thereby rapidly increasing the rate of genetic gain. While these methods are currently in use to breed cattle with superior genetics for production traits, they are not yet available to improve disease resistance. Mycobacterial and mammary gland (Mastitis) infections represent two diseases that severely impact global cattle production, with an annual estimated cost of \$3 billion for TB globally (1) and multiples of this value for mastitis, especially when the costs associated with subclinical infection are included (2, 3). Both these diseases, as well as complementary analyses on the regulation of bovine immunity, are addressed by cutting edge papers in this edition of *Frontiers*.

*Mycobacterium tuberculosis* causes TB in human beings and over one-third of the world's population are infected with this disease. Similarly, in cattle, related mycobacterial species cause potentially zoonotic infections, which are endemic in many parts of the world, despite stringent global surveillance and control programs. The advent of Next-Generation Sequencing (NGS) technologies holds significant promise to overcome limitations in current generation test sensitivity and specificity and as discussed by McLoughlin et al., transcript biomarker signatures have been identified, which discriminate between TB patients with active and latent disease. Indeed, McLoughlin et al. used multi-dimensional scaling to unambiguously classify peripheral blood leukocyte transcriptomic profiles from *Mycobacterium bovis* infected and healthy cattle, thereby uncovering potential biomarkers for *M. bovis* infection (4). This technology is a powerful tool for unraveling the complexities of host immune response and provides new layers of information, which deepens our understanding of host–pathogen interactions that underlie Mycobacterial disease pathogenesis. The macrophage is recognized as the key effector cell driving antimycobacterial immunity but which can be hijacked by mycobacteria, thereby contributing to suboptimal bacterial clearance and disease chronicity. Using the macrophage as a model, RNA-seq was performed after challenge with *Mycobacterium avium* subspecies *paratuberculosis* (MAP), which has uncovered novel genes that have not previously been associated with the host response to MAP infection (5). One of the key challenges with NGS technologies is the bioinformatic analysis to extract meaningful and biologically relevant findings from the wealth of data generated. Sophisticated system-biology tools have been employed by Killick et al. to generate biological interaction networks that usefully identify key hub and bottleneck genes that are central to the immune response and thereby potential targets for immunomodulation – either naturally by pathogens in their eternal quest for survival, or therapeutically with the development of new intervention strategies (6). Furthermore, a comparative analysis of the macrophage response to various strains of mycobacteria has been reviewed; *M. bovis* induced a distinct transcriptional profile in monocytederived macrophages compared to the more similar profiles of both *M. bovis* BCG and MAP (7). The authors describe how differential expression of type-I interferon genes were specific to the virulent *M. bovis* strain supporting a role for these genes in the establishment of active tuberculosis in cattle.

The identification of the genes and pathways involved in orchestration of the immune response is not merely of fundamental importance but differentiating between the immune responses to closely related bacterial species can have very real implications for the success of current generation diagnostics. The study by Kennedy et al., examines the effects of annual mandatory testing for *M. bovis* infection on the ELISA performance routinely used for the diagnosis of MAP (8). In one herd, prior to testing for bovine tuberculosis (BTB), 7.9% of cattle serum samples and 5.8% of milk samples were positive for MAP antibodies. Shortly after the BTB test, the MAP ELISA positive rate increased to 39% in both sample types, clearly showing BTB test interference in MAP ELISA performance. Exploiting differences in host immunity induced in response to these closely related strains of mycobacteria could yield significant dividends in terms of increased specificity of diagnostics.

An optimal and effective host immune response must overcome pathogen-mediated efforts to subvert it while minimizing tissue damage and, therefore, the regulation of the immune response is critical. Furthermore, shedding light on the differential dialog between the host and pathogenic or comensal bacterial strains is a very exciting area of research. It is of interest, therefore, that Villena et al. review how immunobiotics can dampen TLRmediated inflammation in an intestinal cell line, and postulate how immunoregulatory feeds could be developed in the future (9). Inflammation is a feature of most diseases, and detailed knowledge on the pathways regulating it is critical toward developing effective immunoregulatory approaches in cattle. MiRNAs have also been shown to be powerful regulators of immunity and as reviewed by Lawless et al., 793 miRNAs have been identified to date in the bovine genome (10). Their rapid induction in response to challenge and the tissue-specific expression pattern of some has led to speculation on their potential use as diagnostics. Work by the same group has identified miRNA signatures of infection in mammary epithelial cells in response to a common mastitiscausing pathogen. The review highlights relevant studies on how miRNAs regulate the production of IFN-y and TNF, key cytokines in the immune response to TB.

In a comprehensive review, Thompson-Crispi et al. integrates multiple studies on the genetic regulation of the bovine immune response, particularly in relation to mastitis (11). Interestingly, the review discusses earlier work by the same group in which their High Immune Response (HIR) technology was used to identify Immunity<sup>+</sup> sires, daughters of which showed a 44% reduction in mastitis as well as reduced susceptibility to other diseases. The potential consequences of selection for a specific immune phenotype are discussed, and the review signposts how integration of complementary genetic, genomic, and epigenetic data – supplemented with accurate disease and health phenotype information – will enhance our ability to breed cattle with superior disease resistance in the future.

Early fetal mortality is a major contributory factor to poor reproductive outcomes and increased costs, especially in high producing dairy cows. In that regard, the review by Fair is a relevant one. Although immunological analyses in the cow during pregnancy are growing, attempts to evaluate the interaction between the cow and the developing fetus are few. While comparative immunology can shed some light, Fair argues that basic understanding in the bovine are required to more comprehensively understand the complex regulation of local and systemic immunity in the pregnant cow and thereby yield novel solutions to fertility problems (12). Understanding the immune shifts that occur during pregnancy is also critical to understanding the windows of susceptibility that may exist through which susceptibility to infectious diseases could be increased.

Multi-factorial challenges, for example achieving and sustaining excellent animal health, require multifaceted solutions that can only be achieved through intensive integration of knowledge and expertise from a diverse spectrum of research efforts (13). As the physicist Richard J. Feynman once wrote "In order to make progress, one must leave the door to the unknown ajar." There is a lot yet to learn in relation to bovine immunity and, therefore, this e-book is a timely integration of the most current research and scientific thinking on these critical issues and will contribute to the direction of future research in these areas. When dealing with infectious diseases, the old truism is perfectly apt – *ipsa scientia potestas est*. New knowledge also prepares us for the unforeseen challenges of the future.

#### **REFERENCES**


**Conflict of Interest Statement:** The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Received: 26 January 2015; accepted: 04 February 2015; published online: 17 February 2015.*

*Citation: Meade KG (2015) Advances in bovine immunology – new tools and new insights to tackle old foes. Front. Immunol. 6:71. doi: 10.3389/fimmu.2015.00071 This article was submitted to Molecular Innate Immunity, a section of the journal Frontiers in Immunology.*

*Copyright © 2015 Meade. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

### RNA-seq transcriptional profiling of peripheral blood leukocytes from cattle infected with Mycobacterium bovis

**Kirsten E. McLoughlin<sup>1</sup> , Nicolas C. Nalpas <sup>1</sup> , Kévin Rue-Albrecht <sup>1</sup> , John A. Browne<sup>1</sup> , David A. Magee<sup>1</sup> , Kate E. Killick <sup>1</sup> , Stephen D. E. Park <sup>1</sup>† , Karsten Hokamp<sup>2</sup> , Kieran G. Meade<sup>3</sup> , Cliona O'Farrelly <sup>4</sup> , Eamonn Gormley <sup>5</sup> , Stephen V. Gordon6,7 and David E. MacHugh1,7\***

<sup>1</sup> Animal Genomics Laboratory, UCD School of Agriculture and Food Science, University College Dublin, Dublin, Ireland

<sup>2</sup> Smurfit Institute of Genetics, Trinity College Dublin, Dublin, Ireland

<sup>3</sup> Animal and Bioscience Research Department, Animal and Grassland Research and Innovation Centre, Dunsany, Ireland

<sup>4</sup> Comparative Immunology Group, School of Biochemistry and Immunology, Trinity Biosciences Institute, Trinity College Dublin, Dublin, Ireland

<sup>5</sup> Tuberculosis Diagnostics and Immunology Research Centre, UCD School of Veterinary Medicine, University College Dublin, Dublin, Ireland

<sup>6</sup> UCD School of Veterinary Medicine, University College Dublin, Dublin, Ireland

<sup>7</sup> UCD Conway Institute of Biomolecular and Biomedical Research, University College Dublin, Dublin, Ireland

#### **Edited by:**

Uday Kishore, Brunel University, UK

#### **Reviewed by:**

Divyendu Singh, Indiana University Bloomington, USA Jaya Talreja, Wayne State University, USA

#### **\*Correspondence:**

David E. MacHugh, Animal Genomics Laboratory, Veterinary Sciences Centre, University College Dublin, Belfield, Dublin 4, Ireland e-mail: david.machugh@ucd.ie

#### **†Present address:**

Stephen D. E. Park, IdentiGEN Ltd., Dublin, Ireland

Bovine tuberculosis, caused by infection with Mycobacterium bovis, is a major endemic disease affecting cattle populations worldwide, despite the implementation of stringent surveillance and control programs in many countries. The development of highthroughput functional genomics technologies, including gene expression microarrays and RNA-sequencing (RNA-seq), has enabled detailed analysis of the host transcriptome to M. bovis infection, particularly at the macrophage and peripheral blood level. In the present study, we have analyzed the peripheral blood leukocyte (PBL) transcriptome of eight natural M. bovis-infected and eight age- and sex-matched non-infected control Holstein-Friesian animals using RNA-seq. In addition, we compared gene expression profiles generated using RNA-seq with those previously generated using the high-density Affymetrix® GeneChip® Bovine Genome Array platform from the same PBL-extracted RNA. A total of 3,250 differentially expressed (DE) annotated genes were detected in the M. bovisinfected samples relative to the controls (adjusted P-value ≤0.05), with the number of genes displaying decreased relative expression (1,671) exceeding those with increased relative expression (1,579). Ingenuity® Systems Pathway Analysis (IPA) of all DE genes revealed enrichment for genes with immune function. Notably, transcriptional suppression was observed among several of the top-ranking canonical pathways including Leukocyte Extravasation Signaling. Comparative platform analysis demonstrated that RNA-seq detected a larger number of annotated DE genes (3,250) relative to the microarray (1,398), of which 917 genes were common to both technologies and displayed the same direction of expression. Finally, we show that RNA-seq had an increased dynamic range compared to the microarray for estimating differential gene expression.

**Keywords: Mycobacterium bovis, tuberculosis, RNA-seq, biomarker, cattle, microarray, peripheral blood**

#### **INTRODUCTION**

Bovine tuberculosis (BTB) is caused by infection with *Mycobacterium bovis*, a member of the *Mycobacterium tuberculosis* complex (MTC). The mycobacterial species and strains that constitute the MTC can cause tuberculosis in a wide range of mammals and display 99.9% similarity at the nucleotide level (1). Econometric analyses incorporating agricultural production and human health indices have placed BTB as the fourth most important disease of cattle, costing an estimated \$3 billion on a global scale annually (2). Furthermore, as a zoonotic agent,*M. bovis* infection has important implications for human health (3, 4).

*Mycobacterium bovis* is normally transmitted via the inhalation of infectious bacilli, whereby infection is established in the lung. Evasion of host immune defenses enables the pathogen to survive and replicate within phagocytic macrophages, the primary innate immune cell that mediates the response to infection, and can result in dissemination of infection via the lymph system leading to disease progression and pathology (5). Subsequent transmission of the pathogen to susceptible hosts maintains the cycle of infection. The immune response to BTB is complex and is largely characterized by macrophage-mediated development of protective TH1-type responses following initial exposure to the pathogen. It has been reported that the development of disease involves a transition from TH1 to non-protective TH2-type responses (6). The progression of infection may also be due to the modulation and suppression of specific immune mechanisms by the pathogen (5).

In many developed countries, control and eradication programs have been put in place to facilitate early detection and removal of infected animals. In Ireland, the test and slaughter policy was introduced in the early 1950s as part of the national BTB eradication scheme (7). This policy includes mandatory screening of animals in the national herd using the single intradermal comparative tuberculin test (SICTT), alone or in combination with *in vitro* ELISA-based interferon-gamma (IFN-γ) release assays in infected herds, to increase the sensitivity of diagnosis (8). However, due to limitations in both of these tests and also the presence of wildlife reservoirs (including the Eurasian badger, *Meles meles*) *M. bovis* has remained recalcitrant to eradication (6). Consequently, there is a pressing need to develop novel diagnostics for early and reliable detection of *M. bovis* infection in cattle herds.

The availability of a well-annotated bovine genome sequence with concomitant advances in high-throughput genomics technologies offers novel approaches to interrogate and better understand the immune response to *M. bovis* infection. Many transcriptomics studies of the host response to *M. bovis* have involved the analysis of blood-derived RNA from naturally or experimentally infected animals, as previous work has shown that for BTB, host immune responses occurring in peripheral blood reflect those at the primary site of disease (9). Microarray analysis of peripheral blood-derived RNA has shown that transcriptional profiling can unambiguously differentiate animals by disease status and can identify immunomodulatory mechanisms associated with pathology (10–13).

The advent of high-throughput sequencing technologies has given rise to new methods for gene expression analysis based on RNA-sequencing (RNA-seq). This RNA-seq approach has a number of important advantages compared to microarray analysis, including unbiased whole-transcriptome profiling; the characterization and analysis of both sense and antisense transcription and novel transcripts; the identification of mRNA isoforms; increased precision and sensitivity for the quantification of lowly expressed transcripts; and the detection of expressed coding and regulatory DNA sequence variants that can influence phenotype [e.g., disease resistance and susceptibility] (14–16).

To gain a deeper knowledge of the host transcriptional response to *M. bovis* infection, we have used RNA-seq to compare the peripheral blood leukocyte (PBL) transcriptomes of eight animals naturally infected with *M. bovis* and eight non-infected control animals. Differentially expressed (DE) genes identified from this analysis were further investigated using the Ingenuity® Systems Pathway Analysis (IPA) Knowledgebase to detect overrepresented cellular pathways in response to *M. bovis* infection. We also compared the gene expression profiles generated from RNA-seq with data from the Affymetrix® GeneChip® Bovine Genome Array using the same PBL-extracted RNA samples (12).

#### **MATERIALS AND METHODS**

#### **ANIMALS**

The 16 age-matched female Holstein-Friesian animals used in this study have been previously described (12). The eight *M. bovis*infected cattle were selected from a panel of naturally infected animals identified during routine disease surveillance by the Irish Department of Agriculture, Food and the Marine. These animals had a positive result for both SICTT and whole blood IFN-γbased BOVIGAM® assay tests (Prionics AG, Zurich, Switzerland). In addition, *M. bovis* infection was confirmed following detailed post-mortem pathological examination and culture. Non-infected control animals were selected from a herd with no recent history of *M. bovis* infection and were shown to be negative for both the SICTT and IFN-γ tests. All animal procedures detailed were performed according to the provisions of the Cruelty to Animals Act (licenses issued by the Irish Department of Health and Children) and ethics approval for the study was obtained from the UCD Animal Ethics Committee.

#### **BLOOD COLLECTION AND RNA EXTRACTION**

The materials and methods used to isolate and purify PBL-derived RNA from all 16 animals have been described by us previously. Briefly, whole blood was collected from each animal in 8 ml heparin vacutainers® (Becton-Dickinson Ltd., Dublin, Ireland) and RNA extraction was performed within 2 h of blood collection. The complete methods used for blood collection, PBL isolation, and total RNA extraction and purification have been described by us previously (12). RNA quantity and quality checking was performed using the NanoDrop™ 1000 spectrophotometer (Thermo Fisher Scientific,Waltham, MA, USA) and the Agilent 2100 Bioanalyzer using an RNA 6000 Nano LabChip kit (Agilent Technologies, Cork, Ireland).All samples displayed a 260/280 ratio >1.8 and RNA integrity numbers (RIN) >8.0.

#### **RNA-seq LIBRARY PREPARATION**

The laboratory method used to generate RNA-seq libraries was adapted from a protocol previously described by our group (17). In total, 16 strand-specific RNA-seq Illumina® libraries were prepared (i.e., eight libraries each for the infected and control groups) using 1.2µg of total RNA. Total RNA was heated at 65°C for 5 min to disrupt any secondary structure and purification of poly(A) RNA was performed using a Dynabeads® mRNA DIRECT® Micro Kit according to the manufacturer's instructions (Invitrogen™). 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). 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., 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., 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.

To facilitate Illumina® GA adaptor ligation, a single "A" base was added to the 3<sup>0</sup> ends of the blunt-end-repaired cDNA samples. Thirty-two microliters of purified phosphorylated bluntend-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<sup>0</sup> to-5<sup>0</sup> 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<sup>0</sup> -dATP overhang with 1× Quick DNA ligase buffer (New England Biolabs®); 30 nM custom indexed singleread adaptors (see Table S1 in Supplementary Material for barcode index sequences); 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× SYBR® Safe DNA gel stain (Invitrogen™). Gels were electrophoresed at 100V 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 (25µl) were performed and contained 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 consisted of an initial denaturation step of 98°C for 30 s, 18 cycles of 98°C for 10 s, 65°C for 30 s, and 72°C for 30 s, followed by a final

extension step of 72°C for 5 min. PCR products were visualized following electrophoresis on a 2% agarose gel stained with 0.25× SYBR® Safe DNA gel stain (Invitrogen™) and purified to remove PCR-generated adaptor-dimer 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 (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 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.

The libraries were subsequently validated using conventional Sanger sequencing of individual library clones. Library fragments from two libraries were cloned using the Zero Blunt® TOPO® PCR Cloning system according to the manufacturer's instructions (Invitrogen™). Conventional Sanger sequencing of 10 plasmid inserts from each of the 2 libraries confirmed that the RNA-seq libraries contained inserts derived from bovine mRNA. Plasmid sequencing was outsourced (Source Bioscience Ltd., Dublin, Ireland) and sequences generated were validated using BLAST-searching of the DNA sequence database (18).

Cluster generation and sequencing of the pooled RNA-seq library was performed on an Illumina® Cluster Station and Genome Analyzer IIx sequencer according to the manufacturer's instructions. The pooled library was sequenced as single-end read 84-mers using Illumina® version 4.0 sequencing kits and the standard Illumina® Genome Analyzer IIx 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 GSE60265.

#### **BIOINFORMATICS AND STATISTICAL ANALYSIS OF RNA-SEQ AND MICROARRAY DATA**

All the bioinformatics pipeline bash, Perl, and R scripts used for computational analyses were deposited in a GitHub repository at https://github.com/kmcloughlin1 and these analyses were performed on a 32-node Compute Server running Linux Ubuntu (version 12.04.2).

An initial quality check was performed on each of the raw read data files using the FastQC software (version 0.10.1)<sup>1</sup> to determine the best sequence read quality filtering strategy. Subsequently, a custom perl script was used to: (1), deconvolute the pooled libraries into individual libraries of sample sequence reads based on the unique index barcode (allowing up to one mismatch as long as the barcode sequence can be associated to a single unique index barcode); (2) filter out single-end reads containing adaptor sequence (allowing up to three mismatches); and (3) remove

<sup>1</sup>http://www.bioinformatics.babraham.ac.uk/projects/fastqc

single-end reads of poor quality (i.e., reads containing 25% of bases with a Phred quality score below 20). Filtered individual libraries file were checked again with the FastQC software package to confirm sequence read quality. Single-end reads, from each filtered individual sample library, were aligned to the *Bos taurus* reference genome [UMD3.1.73; (19)] using the STAR aligner software package [version 2.3.0] (20).

For each library, raw counts for each annotated gene were obtained using the featureCounts software from the Subread package [version 1.3.5-p4] (21). The featureCounts parameters were set to unambiguously assign uniquely aligned single-end reads in a stranded manner to the exons of genes within the *B. taurus* reference genome annotation (UMD3.1.73 genome annotation).

Differential gene expression analysis was performed using the gene raw counts, obtained from featureCounts, within the Bioconductor edgeR package (22). The differential gene expression pipeline within the edgeR package was customized to (1) filter out all bovine rRNA genes; (2) filter out genes displaying expression levels below the minimally set threshold of one count per million (CPM) in at least eight individual libraries (i.e., equivalent to one group of biological replicates); (3) calculate normalization factors for each library using the trimmed mean of *M*-values method (14); (4) generate the density of counts per gene and multidimensional scaling (MDS) plots based on data from each individually barcoded library (using the Euclidean distance metric); (5) estimate the dispersion parameter for each library using the Cox-Reid method; (6) identify DE genes between infected versus non-infected control samples (i.e., unpaired-sample statistical model) using a negative binomial generalized linear model; and (7) adjust the *P*-value for multiple testing using the Benjamini– Hochberg correction (23) with a false discovery rate (FDR) ≤0.05. Mean fold-changes in gene expression are reported in the main body text as geometric mean values in the *M. bovis*-infected group relative to the control group; for genes displaying reduced relative expression, the negative reciprocal geometric mean fold-changes are given.

The IPA® Knowledgebase<sup>2</sup> was used to identify cellular pathways and gene ontology categories that were overrepresented based on the list of DE genes (*P-*value ≤0.05).

The raw microarray data generated from the same 16 PBLextracted RNA samples were retrieved from the NCBI GEO repository (24) with the accession number GSE33359 (12). The Affymetrix® GeneChip® Bovine Genome Array used to generate these data contains 24,072 probe sets representing more than 23,000 gene transcripts. To compare gene expression profiles from these samples using RNA-seq and the microarray, we first re-analyzed the microarray data using a series of Bioconductor packages (25) and the most recent build of the bovine genome [UMD3.1.73; (19)]. Normalization of raw data was performed using the Factor Analysis for Robust Microarray Summarization (FARMS) algorithm (26). The FARMS algorithm uses only perfect match (PM) probes and a quantile normalization procedure, providing both *P*-values and signal intensities. Normalized data were then further subjected to filtering for informative probes

sets using the informative/non-informative (I/NI) calls unsupervised feature selection criterion implemented in FARMS (27). This defines a probe set as being informative when many of its probes reflect the same change in mRNA concentration across arrays. To compare and contrast the two gene expression technologies, microarray probe sets were first annotated with the corresponding Ensembl ID using the Bioconductor biomaRt package (28). Genes displaying differential expression between control and infected groups were identified using the linear models for microarray data (LIMMA) bioconductor package (29). Following this, a Benjamini–Hochberg multiple-testing correction of ≤0.05 was applied to all DE genes (23) and the Euclidean distance was used as the distance metric for MDS plotting.

#### **RESULTS**

#### **SUMMARY STATISTICS FOR THE RNA-seq DATA**

All 16 RNA-seq libraries were sequenced across one full Illumina® GAIIX flow cell. Deconvolution and filtering of sequence reads to remove adaptor-dimer contamination yielded a mean of 13.2 million reads per individual barcoded RNA-seq sample library. Alignment of the filtered reads to the *B. taurus* UMD3.1.73 genome build yielded a mean of 11.8 million reads (90%) that aligned to unique locations in the bovine genome for each RNAseq library; a mean of 906,679 reads (7%) for each library that aligned to multiple locations in the genome; and a mean of 381,991 reads (3%) for each library that did not align to any genome location (Table S1 in Supplementary Material). Further analysis of the mean 11.8 million reads mapping to unique genome locations demonstrated that 52% of these mean reads were assigned to annotated regions of the genome, which were used to calculate raw counts for each sense gene and subsequently used for downstream bioinformatics and systems analyses; while 48% were not assigned to any annotated genome location or were assigned to overlapping (therefore ambiguous) annotated genomic regions.

#### **GENE EXPRESSION AND IPA ANALYSIS OF SENSE STRAND TRANSCRIPTION**

Analysis of the gene coverage based exclusively on sense strand sequence information, revealed that of the 24,616 annotated *B. taurus* genes in Ensembl (release 73), 17,792 genes (72.3%) had at least one sequence read count (i.e., one mapped read) in at least 1 of the 16 individual sample RNA-seq libraries. The 17,792 detectable genes were further filtered by removing lowly expressed genes, whereby only genes displaying more than 1 CPM reads in 8 or more individual libraries were used for subsequent analyses. This yielded 12,294 genes (49.9% of annotated *B. taurus* genes) that were suitable for downstream analyses.

Prior to differential gene expression analysis, the 12,294 filtered genes were used to generate an MDS plot to visualize gene expression and infection status for the 16 animal PBL samples (**Figure 1A**). This plot shows that samples were clearly differentiated according to infection status along dimension 1 and dimension 2 highlights one *M. bovis*-infected sample as a possible outlier (*Infected 34* – animal ID).

Statistical analysis of all 12,294 genes that passed the filtering process identified a total of 3,250 DE genes (FDR ≤0.05), of which 1,579 and 1,671 displayed increased and decreased expression,

<sup>2</sup>http://www.ingenuity.com

respectively, in the *M. bovis*-infected samples relative to the noninfected control samples (Table S2 in Supplementary Material). Among the DE genes showing the greatest mean fold-change increase in expression were *CXCL6* (+6.77), *IL8* (+5.39), and *CTLA4* (+3.61), while *CXCL10* (−3.27), *DEFB10* (−7.21), and *IL12* (−4.35) showed the greatest mean fold-change decrease in expression; all of these genes have been previously shown to be involved in the host response to mycobacterial infection (12, 30–33).

from 5,082 informative microarray probe sets in 8 M. bovis-infected and 8

Functional categorization of the DE genes using IPA revealed an overrepresentation of genes with roles in inflammatory response, immunological disease and infectious disease. Of the 3,250 genes found to be DE, 2,785 mapped to the IPA Knowledgebase. IPA analysis identified 201 statistically significant (*P-*value ≤0.05) enriched pathways, many of which were associated with immune

function (Table S3 in Supplementary Material). Based on the well-documented role of the T-cell response to mycobacterial infection, the top-ranking IPA canonical pathway (*T-cell receptor signaling* ) was overlaid with the RNA-seq gene expression data, which indicated activation of this pathway (**Figure 2**). Further inspection of the IPA results showed *Leukocyte extravasation signaling* to be the second-ranked canonical pathway. Transcriptional suppression was observed for this pathway with several genes required for migration of leukocytes to the site of infection displaying reduced relative expression (**Figure 3**). These include genes encoding leukocyte ligands required for endothelial adhesion such as *ITGB2* (−1.29-fold), *ITGAL* (−1.23-fold), and *SPN* (also known as *CD43*; −1.42-fold). Reduced relative expression of the gene encoding the LFA-1 protein was also indirectly observed in these BTB-infected animals. LFA-1 is a complex formed from an α-chain encoded by *ITGAL* and a β-chain encoded by *ITGB2*; both *ITGAL* and *ITGB2* exhibited decreased relative expression as discussed above. Decreased relative expression of *PECAM1* (−1.69-fold), which encodes a protein involved in the transmigration of leukocytes through or between endothelial layers into tissues during extravasation, was also detected in the present study (34). Furthermore, a reduction in relative expression was also observed for *MMP9* (−1.79), which degrades the extracellular matrix facilitating the transmigration of leukocytes (35, 36).

#### **COMPARISON OF THE NUMBER OF DIFFERENTIALLY EXPRESSED GENES IDENTIFIED FROM RNA-seq AND MICROARRAY PLATFORMS**

The RNA extracted from these animals has previously been analyzed using theAffymetrix® Bovine GenomeArray (12). To directly compare gene expression profiles generated by the two platforms, we re-analyzed all microarray data and of the 24,072 probe sets represented on the microarray, 5,082 of these passed the filtering process and were designated as informative. An MDS plot generated from these 5,082 informative probe sets shows samples clustered according to their infection status (**Figure 1B**); this pattern was also observed for these samples by Killick and colleagues using hierarchical clustering (12).

Analysis of all informative microarray probe sets identified 2,808 DE transcripts (FDR ≤0.05) and mapping of these transcripts to the *B. taurus* UMD3.1.73 genome build yielded 1,398 DE genes with Ensembl IDs (Table S4 in Supplementary Material). This is substantially lower than the 2,757 DE genes previously reported by us for the same RNA samples (12). This discrepancy may be explained by differences in the versions of the *B. taurus* reference genome used to annotate the microarray probe sets: for the previous study, we used the Btau4.0 genome assembly, while in the current study the UMD3.1.73 genome assembly was used. Discrepancies between the two genome annotations as detailed by Zimin et al. (37) most readily account for the reduced number of DE genes with Ensembl IDs observed for the present study.

Of the 1,398 DE genes obtained from re-analysis of the Affymetrix® Bovine Genome Array in infected animals relative to the control animals, 630 and 768 exhibited increased and decreased expression, respectively. Consequently, it is noteworthy that the number of Ensembl-annotated DE genes obtained

control PBL samples.

with the RNA-seq analysis (3,250 DE genes; 1,579 increased relative expression and 1,671 decreased relative expression) was markedly higher than the number of DE genes detected using the microarray. Further examination of the results showed that 917 DE Ensembl genes had the same direction of expression (i.e., increased or decreased relative expression) on both platforms; 2,331 DE Ensembl genes were unique to the RNA-seq results (i.e., DE using RNA-seq but not DE on the microarray platform); and 479 DE Ensembl genes were unique to the microarray (i.e., DE on the microarray platform but not DE on the RNA-seq platform). Finally, two of the genes were found to have conflicting patterns of expression, i.e., genes that displayed increased relative expression on one platform and decreased relative expression on the other platform. *CHTOP* showed decreased relative expression on the RNA-seq platform but increased relative expression on the

microarray platform and*GPR89* showed increased relative expression on the RNA-seq platform but decreased relative expression on the microarray platform (**Figure 4**).

Finally, analysis and comparison of the number of IPAidentified statistically significant canonical pathways (*P-*value ≤0.05) for the DE genes from both platforms revealed 101 pathways that were common to both, 100 that were unique to RNA-seq, and 36 that were unique to the microarray. The larger number of IPA-identified canonical pathways from the RNA-seq results reflects the greater number of DE genes detected using this platform. Despite this, however, there was a notable level of overlap in the number of pathways identified using both platforms. The percentage of overlapping pathways between the two platforms was estimated as 50.2% for RNA-seq and 73.7% for the microarray (Tables S3 and S5 in Supplementary Material).

**identified using IPA, the Leukocyte extravasation pathway**. Red shading indicates increased expression in M. bovis-infected animals

#### **CORRELATION OF THE LOG FOLD-CHANGE IN GENE EXPRESSION BETWEEN THE TWO PLATFORMS USING ALL GENES THAT PASSED FILTERING**

We analyzed the correlation between the log<sup>2</sup> fold-changes (infected versus control) on both gene expression platforms. We hypothesized that genes with high differential fold-changes in decreased expression in M. bovis-infected animals relative to the non-infected control group.

expression based on RNA-seq analysis should also show high differential fold-changes based on the microarray data. Accordingly, these data would be expected to yield a significant, positive correlation. For this, all 12,294 genes and 5,082 probe sets that passed the filtering criteria across both sample groups for both RNA-seq and microarray analysis, respectively, were considered. Next, we identified all the genes, irrespective of significance, that were common to both data sets based on Ensembl gene ID, this involved matching the Ensembl IDs of the filtered genes in the RNA-seq data set with the Ensembl IDs of the filtered microarray probe sets. In total, 2,265 Ensembl IDs were identified in common between the two platforms. The log<sup>2</sup> fold-change values for these 2,265 genes (infected versus control) for both the RNA-seq and microarray data, irrespective of the significance of differential expression generated a Spearman correlation coefficient of 0.88 (*P* ≥ 0.001), which underlines the reproducibility and robustness of both platforms for investigations of differential gene expression.

#### **COMPARISON OF THE FOLD-CHANGE AND DIFFERENCE IN EXPRESSION VALUES FOR RNA-seq AND MICROARRAY DATA**

We next investigated the correlation between fold-change in expression and gene expression levels. This enabled us to determine if the highest fold-changes in gene expression were observed for genes with low levels of expression. For this, we compared the log<sup>2</sup> expression fold-changes with the log<sup>2</sup> differences in CPM between infected and control groups for the RNA-seq platform. We also compared the log<sup>2</sup> expression fold-changes with the log<sup>2</sup> differences in hybridization intensity between infected and control groups for the microarray platform. We hypothesized that if transcripts with the lowest expression levels gave the highest foldchange values, a negative correlation would be observed between log<sup>2</sup> expression fold-changes and log<sup>2</sup> differences for genes displaying increased relative expression in the infected group relative to the control group. Reciprocally, a positive correlation would be expected between log<sup>2</sup> expression fold-changes and log<sup>2</sup> differences for genes displaying decreased relative expression for the same group contrast.

For the RNA-seq platform, all 12,294 genes that passed filtering were used. Of these, 6,243 displayed an increase in relative expression and 6,051 displayed a decrease in relative expression, irrespective of significance. Spearman rank correlation coefficients of 0.35 and −0.43 were observed for genes displaying increased and decreased relative expression, respectively (*P* ≤ 0.001). For the microarray platform, we used all 5,082 probe sets that passed filtering of which 2,551 displayed an increase in expression and 2,531 displayed a decrease in expression, irrespective of significance. Spearman rank correlation coefficients of 0.61 and −0.64 were observed for genes displaying increased and decreased relative expression, respectively (*P* ≤ 0.001). The correlation coefficients for both the RNA-seq and microarray platforms for genes with increased and decreased relative expression do not support the hypothesis detailed above; therefore, we conclude that there is no obvious relationship between gene expression level and fold-change in expression.

#### **DYNAMIC RANGE OF RNA-seq AND MICROARRAY DATA**

To investigate the dynamic range of the RNA-seq and microarray platforms, the log<sup>2</sup> reads per kilobase per million (RPKM) from the RNA-seq data and the log<sup>2</sup> intensities from the microarray data were used; only genes and probe sets that passed the filtering criteria (12,294 and 5,082, respectively) were considered for this analysis. The lowest expression value was subtracted from the highest expression value for each platform. For the RNA-seq platform, the gene displaying the lowest expression level was *MUC5B* (log<sup>2</sup> RPKM of −6.49), while the gene with the highest value was *COX1* (log<sup>2</sup> RPKM of 13.46); this yielded a log<sup>2</sup> dynamic range of 19.96. It is important to note that as RPKM values are proportions, values <1.0 will yield a negative result when log-transformed. For the microarray platform, the probe set with the lowest log<sup>2</sup> intensity was Bt.29403.1.s1\_at (4.89) and the probeset with the highest log<sup>2</sup> intensity was AFFX-Crex-5\_at (13.90), yielding a log<sup>2</sup> dynamic range of 9.01. Therefore, the dynamic range of RNA-seq for the current study is almost 2,000-fold greater than the microarray.

#### **DISCUSSION**

Whole-genome transcriptional profiling has been successfully used to study human and BTB and has facilitated high-resolution analysis of the host genes and cellular pathways that are activated and perturbed in response to mycobacterial pathogens (11–13, 17, 38–42). The target tissue for studies of the host response has generally been peripheral blood collected from infected and noninfected individuals or animals; peripheral blood provides an easily accessible biological sample that reflects the host immunological and pathological changes induced at the site of infection (9). For

example, Berry and colleagues identified a microarray-derived 393 transcript biomarker signature, characterized by type 1 interferoninducible genes, which discriminated active human TB patients from latently infected and healthy control individuals (38). Furthermore, microarray-based comparative analysis of the peripheral blood transcriptome of active human tuberculosis cases and sarcoidosis patients (an analogous granulomatous disease of the respiratory tract of unknown etiology) also revealed a transcriptional signature that differentiated these two pathologically similar diseases (39). Similarly for BTB, pan-genomic and immunospecific microarray analysis of peripheral blood has demonstrated that *M. bovis*-infected and non-infected animals can be unambiguously differentiated by disease status. Moreover, downstream analysis of the DE genes provided functional genomics evidence that active BTB is associated with the suppression of host innate immune responses and impairment of T-cell signaling (10, 12).

Notwithstanding the remarkable progress in functional genomics studies of the immunobiology and host–pathogen interactions, microarray studies of the host response to mycobacterial infections are not without limitations. For example, microarrays are limited to analysis of genes/transcripts for which probes can be generated from functionally annotated genome resources. In addition, quantification of gene expression indirectly from hybridization signal intensities constrains the dynamic range for lowly and highly expressed gene transcripts. Also, the measurement of gene expression can be hampered by probes that differ in their hybridization affinities with the target mRNA and by background non-specific hybridization, particularly for lowly expressed genes (15, 43, 44).

In contrast, RNA-seq technologies, which are based on highthroughput sequencing and subsequent counting of all expressed RNA transcripts present in a biological sample, have several advantages for quantifying RNA abundance and unraveling the complexity of the host transcriptome following mycobacterial infection. These include unbiased global gene expression analysis (the entire transcriptome is normally surveyed), detection of allele-specific expression, and cataloging of novel transcripts, RNA classes (e.g., long non-coding RNA transcripts), and splice variants that are rarely quantifiable using microarray technologies (15, 44). Also, RNA-seq analysis is based on digital counts of reads that map to annotated genes within a reference genome, thus offering a more precise and sensitive method to identify and quantify DE genes than analog microarray hybridization intensities. Consequently, for the present study we have compared the peripheral blood transcriptomes of non-infected control animals and animals naturally infected with *M. bovis* (eight samples per group) using RNA-seq and compared these results to parallel data obtained using the pan-genomic Affymetrix® Bovine Genome Array.

#### **FUNCTIONAL BIOLOGY OF RNA-seq RESULTS**

A mean of 13.2 million 78-mer reads per individually barcoded RNA-seq library was obtained and this yielded a mean of 1.03 Gb of sequence data per library. A mean of 11.8 million reads (89% of the mean number of library reads) per library mapped to unique locations in the bovine UMD3.1 reference genome. This is higher than previously reported comparable studies: Nalpas and colleagues observed that 63.6% of reads mapped to unique UMD3.1 regions and Churbanov and colleagues observed that 71.2% of reads mapping to the Btau 4.0 reference genome (17, 45). Analysis of the gene coverage from these uniquely mapped reads revealed that 17,792 genes from a total of 24,616 annotated bovine genes gave at least 1 sequence read in at least 1 of the 16 individual libraries. As accurate quantification of gene expression is reliant on sequencing depth, whereby low sequencing coverage can lead to the generation of false-positive DE genes (type I errors) (46, 47), we have used a stringent sequencing-depth filtering criterion to remove lowly expressed genes. This filtering involved the retention of genes displaying more than one CPM in eight or more individual libraries; 12,294 genes (49.9% of annotated *B. taurus* genes) were retained for downstream differential expression analyses. The number of biological replicates (*n* = 8) and mean sequencing depth per library achieved in the present study is sufficient for the accurate quantification and analysis of DE genes with a corresponding reduction in the number of type I errors due to lowly expressed genes (46–48).

Previous work by our group and others has demonstrated that a functional genomics approach can highlight novel aspects of the complex etiology of *M. bovis* infection in cattle. These studies have confirmed the role of TH1-type cytokines and chemokines and innate immune receptors (e.g., TLR genes) in mediating the response to *M. bovis* infection. Moreover, these investigations support the hypothesis that the immunoevasive mechanisms used by the pathogen during infection are reflected in the host transcriptome at the peripheral blood level; in particular, the suppression of innate immune signaling, which leads to an inferior adaptive immune response (10, 12).

In the current study, RNA-seq was used to identify a total of 3,250 DE genes in the infected group relative to the control group; of these, 1,579 genes displayed increased relative expression and 1,671 genes showed reduced relative expression. The number of DE genes identified here through RNA-seq analysis exceeds the number of DE gene previously reported by Killick et al. (12) for the same RNA samples (2,757 DE genes; 1,281 and 1,476 genes displaying increased and decreased relative expression, respectively). This finding emphasizes the increased sensitivity of RNA-seq compared to microarrays for studies of differential gene expression (17, 49).

Gene ontology analysis revealed enrichment for genes involved in inflammation and immunity. TH1-type cytokines and chemokines, such as *CXCL6* and *IL8*, were among the top-ranking DE genes based on fold-change in expression; additional innate immune genes, such as *CCL4*, *CXCR4*, *CXCR7*, *IL1A*, *IL8*, *IL10*, and*TLR4* also shown increased expression in the *M. bovis*-infected group relative to the controls. Several innate immune genes also displayed reduced relative expression including *CXCL10*, *DEFB10*, *IL12*, *IL18*, and *IL27*. These results suggest that although innate immune genes play a role in mediating the host response to *M. bovis* infection, these genes may also serve as targets for immunomodulation by the pathogen to facilitate survival in the host. For example, *IL12*, *IL18*, and *IL27* encode cytokines that have all been shown to play key roles in initiating and controlling the adaptive immune response to mycobacterial infection (31, 50); suppression of these genes may result in the development of an inferior cellular response to infection leading to disease progression. The increased relative expression of the *IL10*

gene, which encodes an immunosuppressive cytokine, may also result in the suppression of host innate immune responses to infection resulting in mycobacterial persistence within the host (51). Collectively, these findings support our previous work,which hypothesized that the suppression of innate immune expression and signaling limits the initiation and maintenance of an appropriate adaptive immune response, contributing to the progression of BTB disease (10, 12).

Further analysis of the DE genes using the IPA Knowledgebase identified additional cellular mechanisms within several of the top-ranking canonical pathways, which may be subject to immunomodulation by the pathogen, including the *Leukocyte Extravasation Signaling* and *Tec Kinase Signaling* pathways. The *Leukocyte Extravasation Signaling* pathway exhibited a decrease in expression for many genes encoding positive modulators of this pathway, including *SPN* (also known as *CD43*), *ITGAL*, *ITGB2*, and *PECAM1* (52–55). Leukocyte extravasation refers to the transendothelial migration of activated leukocytes from the blood into infected tissue and is vital for immune surveillance and defense (56). This process, which requires the adherence of leukocytes to the endothelial surface of blood vessels followed by transmigration through the endothelial blood vessel cell layer into the infected tissue, is mediated by chemokines and several cell surface proteins and adhesion molecules including selectins and integrins (57). Within the IPA-identified *Leukocyte Extravasation Signaling* pathway, transcriptional suppression was observed for several leukocyte ligands required for endothelial adhesion during extravasation. *SPN*, *ITGAL*, and *ITGB2* encode the CD43 and the CD11b and CD18 leukocyte cell surface ligands, respectively [the latter of which can complex with different protein partners to form different integrins such as LFA-1 and MAC1 (58)], which are required for leukocyte adhesion to the endothelial cells and subsequent transmigration of the leukocytes into infected tissue (54, 59, 60). *PECAM1* encodes a selectin protein found at intercellular endothelial junctions and is also required for transmigration across these barriers (52, 61). Lower expression of these ligands may result in reduced leukocyte recruitment to the site of infection, leading to an impaired adaptive immune response to contain or eradicate *M. bovis* infection in the host, ultimately leading to disease progression. In addition, these findings lend further support to the hypothesis that the immunoevasion mechanisms used by the pathogen are reflected in the host transcriptome.

Hematological analysis of blood samples taken from the animals analyzed in the current study showed a significant increase in the mean number of lymphocytes (*P* = 0.001) and a significant decrease in the mean number of monocytes (*P* = 0.002) for the infected animals relative to the control group (12). Conversely, no significant differences were observed in the mean number of neutrophils between the two sample groups (*P* ≥ 0.05). It is therefore likely that many of the gene expression changes observed in the current study reflect differences in white blood cell populations between the sample groups.

#### **WHOLE-GENOME EXPRESSION PROFILING: RNA-seq VERSUS MICROARRAYS**

The PBL-extracted RNA analyzed in the present study using RNAseq had also been previously analyzed by us using the Affymetrix® GeneChip® Bovine Genome Array (12); this enabled a direct technical comparison between the two platforms. RNA-seq analysis detected 2.3-fold (3,250/1,398) more DE genes with Ensembl IDs compared to the microarray platform. The percentage of overlapping DE genes with Ensembl IDs between the two platforms was estimated at 28.2% for RNA-seq (917/3,250) and 65.6% (917/1,398) for the microarray. Similarly, the concordance rate based on IPA-identified canonical pathways for the two platforms was estimated at 73.7% for the microarray and 50.2% for RNA-seq. These results demonstrate that the majority of DE genes detected using the microarray are also detected by RNA-seq. This finding also highlights the greater number of DE genes uniquely identified by RNA-seq compared to microarrays as previously reported by us and others (17, 62, 63).

In the current study, the greater number of DE genes detected using RNA-seq can be largely attributed to the greater dynamic range of RNA-seq, which enables sensitive detection of lowly, but DE genes between the infected and control groups (46, 47). Notably, the concordance rate for the microarray (65.6%) is higher than that previously reported by our group for monocyte-derived macrophages (MDM) infected *in vitro* with *M. bovis* (17) and by (64), who examined the transcriptome of anti-CD3- and anti-CD28-stimulated human CD4<sup>+</sup> T cells. The increased microarray concordance rate (based on the number of DE genes) observed for the current study compared to these previous studies is likely due to increased sequencing depth and a greater number of biological replicates (46–48).

Interestingly, two DE genes (0.06% of all RNA-seq DE genes and 0.14% of all microarray DE genes) displayed opposite directions of expression on the two platforms. These discordance rates are lower than that previously reported for RNA-seq and microarray analysis of human cancer cell transcriptomes (65) and may be explained by several technical factors including random error, differences in the transcript isoforms detected by both platforms, and the susceptibility of microarray probes to cross-hybridize with non-specific gene transcripts (66, 67).

#### **WHOLE-GENOME TRANSCRIPTOMICS: BIOMARKER DEVELOPMENT FOR M. BOVIS INFECTION**

Multidimensional scaling analysis using all RNA-seq genes that passed the filtering criteria unambiguously differentiated animals on the basis of their disease status (**Figure 1A**). This result is also supported by the microarray data generated (12) and re-analyzed here. These findings suggest that genome-wide expression profiling of peripheral blood from *M. bovis*-infected animals can be used to identify transcriptional biomarkers for the detection of infected animals within herds and thereby augment surveillance strategies in countries where BTB control programs have been implemented. In addition, recent work has demonstrated that circulating serum or plasma microRNAs may serve as a complementary source of robust biomarkers for tuberculosis and other infectious diseases (68–72). Notwithstanding this, further work using large PBL sample panels from additional animals infected with *M. bovis* and other microbial pathogens will be required to identify and validate robust *M. bovis*-specific transcriptional signatures of infection.

### **ACKNOWLEDGMENTS**

This work was supported by Investigator Grants from Science Foundation Ireland (Nos: SFI/01/F.1/B028 and SFI/08/IN.1/B2038), a Research Stimulus Grant from the Department of Agriculture, Food and the Marine (No: RSF 06 405), a European Union Framework 7 Project Grant (No: KBBE-211602-MACROSYS), and the UCDWellcome Trust Computational Infection Biology Ph.D. Programme (No: 097429/Z/11/Z). We would like to thank Anthony Duignan (DAFM) and staff at the Department of Agriculture, Food and the Marine Backweston Laboratory Campus, Celbridge, Co. Kildare for assistance with cattle blood sampling. We also thank Mairéad Doyle and Tara Fitzsimons of the UCD Tuberculosis Diagnostics and Immunology Research Centre for performing IFN-γ tests.

#### **SUPPLEMENTARY MATERIAL**

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

#### **REFERENCES**


bovine tuberculosis and the identification of potential diagnostic biomarkers. *PLoS One* (2012) **7**:e30626. doi:10.1371/journal.pone.0030626


counter-receptor in human pre-B-cell leukemia NALL-1. *Cancer Res* (2008) **68**:790–9. doi:10.1158/0008-5472.CAN-07-1459


transcription profiling. *BMC Genomics* (2010) **11**:282. doi:10.1186/1471-2164- 11-282


**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.

*Received: 02 July 2014; accepted: 04 August 2014; published online: 26 August 2014. Citation: McLoughlin KE, Nalpas NC, Rue-Albrecht K, Browne JA, Magee DA, Killick KE, Park SDE, Hokamp K, Meade KG, O'Farrelly C, Gormley E, Gordon SV and MacHugh DE (2014) RNA-seq transcriptional profiling of peripheral blood leukocytes from cattle infected with Mycobacterium bovis. Front. Immunol. 5:396. doi: 10.3389/fimmu.2014.00396*

*This article was submitted to Molecular Innate Immunity, a section of the journal Frontiers in Immunology.*

*Copyright © 2014 McLoughlin, Nalpas, Rue-Albrecht, Browne, Magee, Killick, Park, Hokamp, Meade, O'Farrelly, Gormley, Gordon and MacHugh. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

## Analysis of the bovine monocyte-derived macrophage response to Mycobacterium avium subspecies paratuberculosis infection using RNA-seq

#### **Maura E. Casey 1,2, Kieran G. Meade<sup>2</sup> , Nicolas C. Nalpas <sup>1</sup> , Maria Taraktsoglou<sup>3</sup> , John A. Browne<sup>1</sup> , Kate E. Killick 1,4, Stephen D. E. Park <sup>1</sup>† , Eamonn Gormley <sup>5</sup> , Karsten Hokamp<sup>6</sup> , David A. Magee<sup>1</sup>† and David E. MacHugh1,7\***

<sup>1</sup> Animal Genomics Laboratory, UCD School of Agriculture and Food Science, University College Dublin, Dublin, Ireland

<sup>2</sup> Animal and Bioscience Research Department, Animal and Grassland Research and Innovation Centre, Teagasc, Dunsany, Ireland

<sup>3</sup> Biological Agents Unit, Health and Safety Executive, Leeds, UK

<sup>4</sup> Systems Biology Ireland, UCD Conway Institute of Biomolecular and Biomedical Research, University College Dublin, Dublin, Ireland

<sup>5</sup> Tuberculosis Diagnostics and Immunology Research Centre, UCD School of Veterinary Medicine, University College Dublin, Dublin, Ireland

<sup>6</sup> Smurfit Institute of Genetics, Trinity College Dublin, Dublin, Ireland

<sup>7</sup> UCD Conway Institute of Biomolecular and Biomedical Research, University College Dublin, Dublin, Ireland

#### **Edited by:** Uday Kishore, Brunel University, UK

**Reviewed by:**

Anthony George Tsolaki, Brunel University, UK Annapurna Nayak, Brunel University, UK

#### **\*Correspondence:**

David E. MacHugh, Animal Genomics Laboratory, UCD Veterinary Sciences Centre, University College Dublin, Dublin 4, Ireland e-mail: david.machugh@ucd.ie

#### **†Present address:**

Stephen D. E. Park, IdentiGEN Ltd., Dublin, Ireland; David A. Magee, Department of Animal Science, University of Connecticut, Storrs, CT, USA

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 ageand 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.

**Keywords: cattle, immune response, Johne's disease, macrophage, microarray, Mycobacterium avium subspecies paratuberculosis, RNA-sequencing, transcriptome**

#### **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 NODlike receptors (NLRs) (16, 17); indeed, it has been demonstrated that TLR2, TLR4, and NOD2 can independently recognize MAP cellular components (18). Infected macrophages secrete proinflammatory cytokines, such as IL-1B and TNF, which activate an early protective TH1 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 monocytederived 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). Thirtyfive 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 H2O.

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<sup>0</sup> -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)<sup>1</sup> to determine the most appropriate read quality filtering methodology. Consequently, a custom perl script was used to deconvolute sequence

<sup>1</sup>www.bioinformatics.babraham.ac.uk/projects/fastqc

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<sup>0</sup> -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<sup>0</sup> -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 postinfection were used for downstream systems analysis to identify important cellular pathways with the Ingenuity® Systems Pathway Analysis Knowledgebase (IPA<sup>2</sup> ;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.

<sup>2</sup>http://www.ingenuity.com

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.

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 MAPinfected and control MDM samples at 6 hpi ranked by fold-change and with FDR-adjusted *P* values ≤0.05.

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.


**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**.


**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**.

#### **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 postinfection time points and are, for the most part, associated with immunobiology.

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).

#### **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

expression in MAP-infected MDM relative to the non-infected control MDM.

non-expression and non-differential expression, respectively.

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 noninfected 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<sup>2</sup> reads per kilobase per million mapped reads (RPKM) from the RNA-seq data and the log<sup>2</sup> 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<sup>2</sup> dynamic range of 25.31 was estimated based on the *FAT3* gene (ENSBTAG00000004081, log<sup>2</sup> RPKM = −9.15) and the *FTH1* gene (ENSBTAG00000011184, log<sup>2</sup> RPKM = 16.16). For the microarray platform, this calculation yielded an estimated log<sup>2</sup> dynamic range of 13.56 based on the *ZCCHC8* gene (ENSBTAG00000006114, log<sup>2</sup> intensity = 2.03) and the *B2M* gene (ENSBTAG00000012330, log<sup>2</sup> 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<sup>2</sup> EXPRESSION VALUES AND log<sup>2</sup> FOLD-CHANGES BETWEEN THE RNA-seq AND MICROARRAY PLATFORMS**

We next examined the correlation between the log<sup>2</sup> expression values generated using the RNA-seq (log<sup>2</sup> RPKM values) and microarray (log<sup>2</sup> 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−15) and 0.67 (*P* < 1.0 × 10−15) 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−15) for the MAPinfected sample group and 0.66 (*P* < 1.0 × 10−15) for the control sample group.

Following this, log<sup>2</sup> 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<sup>2</sup> expression fold-change values for RNA-seq and the microarray platform at both 2 hpi (0.46; *P* < 1.0 × 10−15) and 6 hpi (0.52; *P* < 1.0 × 10−15).

The results of these analyses, using both log<sup>2</sup> expression values and log<sup>2</sup> 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 preexisting 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 MAPinfected 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<sup>2</sup> 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<sup>2</sup> foldchange = +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 MAPinfected 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<sup>2</sup> 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<sup>2</sup> foldchange = +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<sup>2</sup> 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 proinflammatory 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<sup>2</sup> fold-change = +4.9; *IL1B* log<sup>2</sup> 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<sup>2</sup> fold-change at 2 hpi = +2.7) and *IRAK2* (log<sup>2</sup> fold-change at 2 hpi = +2.2) are important components of the NLRP3-inflammasome complex. Indeed, *Nlrp3*−/<sup>−</sup> 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 wellcharacterized 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<sup>2</sup> fold-change = +1.1 (2 hpi), +1.9 (6 hpi); *IL10* log<sup>2</sup> 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 proinflammatory 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 RNAseq (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 RNAseq-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<sup>0</sup> end of genes and are therefore 3<sup>0</sup> 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 RNAseq 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.

#### **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

#### **REFERENCES**


monocyte-derived macrophages in response to *in vitro* challenge with*Mycobacterium bovis*. *PLoS One* (2012) **7**:e32034. doi:10.1371/journal.pone.0032034


**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.

*Received: 15 September 2014; accepted: 10 January 2015; published online: 04 February 2015.*

*Citation: Casey ME, Meade KG, Nalpas NC, Taraktsoglou M, Browne JA, Killick KE, Park SDE, Gormley E, Hokamp K, Magee DA and MacHugh DE (2015) Analysis of the bovine monocyte-derived macrophage response to Mycobacterium avium subspecies paratuberculosis infection using RNA-seq. Front. Immunol. 6:23. doi: 10.3389/fimmu.2015.00023*

*This article was submitted to Molecular Innate Immunity, a section of the journal Frontiers in Immunology.*

*Copyright © 2015 Casey, Meade, Nalpas, Taraktsoglou, Browne, Killick, Park, Gormley, Hokamp, Magee and MacHugh. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

## Key hub and bottleneck genes differentiate the macrophage response to virulent and attenuated Mycobacterium bovis

#### **Kate E. Killick 1,2\*, David A. Magee<sup>1</sup>† , Stephen D. E. Park 1,3, Maria Taraktsoglou1,4, John A. Browne<sup>1</sup> , Kevin M. Conlon5,6, Nicolas C. Nalpas <sup>1</sup> , Eamonn Gormley <sup>7</sup> , Stephen V. Gordon5,8, David E. MacHugh1,8 and Karsten Hokamp<sup>9</sup>**

<sup>1</sup> Animal Genomics Laboratory, UCD School of Agriculture and Food Science, University College Dublin, Dublin, Ireland

<sup>2</sup> Systems Biology Ireland, UCD Conway Institute of Biomolecular and Biomedical Research, University College Dublin, Dublin, Ireland

3 IdentiGEN Ltd., Dublin, Ireland


<sup>6</sup> Science Foundation Ireland (SFI), Dublin, Ireland

<sup>7</sup> Tuberculosis Diagnostics and Immunology Research Centre, UCD School of Veterinary Medicine, University College Dublin, Dublin, Ireland

<sup>8</sup> UCD Conway Institute of Biomolecular and Biomedical Research, University College Dublin, Dublin, Ireland

<sup>9</sup> Smurfit Institute of Genetics, Trinity College, Dublin, Ireland

#### **Edited by:**

Kieran G. Meade, Teagasc, Ireland

#### **Reviewed by:**

Geanncarlo Lugo-Villarino, Centre National de la Recherche Scientifique, France Silvia Bulfone-Paus, University of Manchester, UK Graham Stewart, University of Surrey, UK

#### **\*Correspondence:**

Kate E. Killick, Systems Biology Ireland, UCD Conway Institute of Biomolecular and Biomedical Research, University College Dublin, Dublin 4, Ireland e-mail: kate.killick@ucd.ie

**†Present address:**

David A. Magee, Department of Animal Science, University of Connecticut, Storrs, CT, USA

**INTRODUCTION**

*Mycobacterium bovis*, the causative agent of bovine tuberculosis (BTB) is a facultative pathogen that has a genome sequence 99.95% identical to that of *M. tuberculosis*, the etiological agent responsible for human tuberculosis (TB) (1). The pathogenesis of BTB is broadly similar to that of human TB and many of the characteristics of *M. tuberculosis* infection are thought to apply to *M. bovis* infection in cattle (2, 3). Both pathogens enter the host principally via inhalation of infected aerosol droplets, where they are phagocytosed by alveolar macrophages and dendritic cells. Pathogenic mycobacteria have evolved mechanisms for subverting the host immune response, including prevention of phagosome maturation and subsequent lysosomal delivery (4), enabling survival and replication within the host macrophage following phagocytosis (5). Infected alveolar macrophages can then remain within the central core of granulomas during latent infection or disseminate to draining lymph nodes and to other host organs during active infection (6, 7).

Mycobacterium bovis is an intracellular pathogen that causes tuberculosis in cattle. Following infection, the pathogen resides and persists inside host macrophages by subverting host immune responses via a diverse range of mechanisms. Here, a high-density bovine microarray platform was used to examine the bovine monocyte-derived macrophage transcriptome response to M. bovis infection relative to infection with the attenuated vaccine strain, M. bovis Bacille Calmette–Guérin. Differentially expressed genes were identified (adjusted P-value ≤0.01) and interaction networks generated across an infection time course of 2, 6, and 24 h. The largest number of biological interactions was observed in the 24-h network, which exhibited scale-free network properties. The 24-h network featured a small number of key hub and bottleneck gene nodes, including IKBKE, MYC, NFKB1, and EGR1 that differentiated the macrophage response to virulent and attenuated M. bovis strains, possibly via the modulation of host cell death mechanisms. These hub and bottleneck genes represent possible targets for immuno-modulation of host macrophages by virulent mycobacterial species that enable their survival within a hostile environment.

**Keywords: Mycobacterium bovis, tuberculosis, BCG, cattle, macrophage, network, gene interaction network**

*Mycobacterium bovis* Bacille Calmette–Guérin (BCG) is an attenuated strain of *M. bovis* that has been used as a live vaccine for TB control for almost a century. It is estimated that over 100 million people are immunized with BCG annually, making it the most widely used vaccine in human populations (8). All BCG strains available today are derived from a virulent *M. bovis* strain that was isolated at the turn of the twentieth century from the milk of an infected cow with tubercular mastitis (9). This strain, named "lait Nocard," was transferred to the Institut Pasteur in Lille in 1901, where it was cultured on glycerol-soaked potato slices supplemented with ox bile by Albert Calmette and Camille Guérin as part of their studies of TB pathogenesis. Serial passage of the lait Nocard strain on this culture medium resulted in the isolation of a strain, in 1908, that displayed reduced virulence in guinea pigs and calves (9). Over the next 13 years, Calmette and Guérin performed 230 *in vitro* serial passages of this partially attenuated *M. bovis* strain, eventually leading to the first successful vaccination of a human infant with a fully attenuated BCG isolate in July 1921 (10, 11).

All currently existing BCG substrains originate from the initial attenuated strain developed by Calmette and Guérin during this 13-year period. It was during this time that BCG underwent the first of two attenuation phases. Early studies identified a region of difference 1 (RD1) that is deleted from all BCG strains (12, 13) and was likely lost during the initial attenuation phase between 1908 and 1921. RD1 encodes, among other genes, the early secretory antigen target (ESAT-6) protein secretion system (ESX-1), which is present in virulent strains of *M. bovis* and *M. tuberculosis* (14, 15). A second attenuation phase occurred after 1924 due to dissemination of the BCG vaccine from the Institut Pasteur to other countries (13). During this time further genomic deletions occurred, including the RD2 deletion (8, 9).

Virulent *M. bovis* and attenuated BCG differ in their ability to generate disease. Infection with *M. bovis* can lead to subversion of the host immune response and formation of tuberculous granulomas with characteristic BTB pathology (16). In contrast, the BCG vaccine strain does not generally cause progressive infection in immune-competent hosts, but can induce a protective proinflammatory TH1 response via activation of interferon-gamma (IFN-γ) secreting CD4<sup>+</sup> T-cells (17); however, the precise mechanisms of how protection is afforded against virulent challenge remain unclear.

Despite decades of research, the host innate immune mechanisms elicited following infection with virulent *M. bovis* relative to attenuated BCG remain largely unknown. To our knowledge, a systems biology approach to understanding innate immune responses to virulent *M. bovis* and attenuated BCG has not previously been described. Systems analysis of transcriptomics data generated from host macrophages infected with the two strains enables reconstruction of dynamic molecular interaction networks that capture the differential innate immune responses to virulent *M. bovis* and attenuated BCG (18). Network features such as hub and bottleneck genes – a consequence of scale-free behavior – can also be defined, cataloged, and further explored to identify novel biological attributes relevant to the pathogenesis of mycobacterial infections.

Previous studies have demonstrated that the bovine monocytederived macrophage (MDM) is a useful model cell type for understanding the early host immunogenetic response to *M. bovis* and *M. bovis*-derived antigens (19–22). In the present study, the Affymetrix® GeneChip® Bovine Genome Array was used to examine the transcriptome of MDM from seven age-matched Holstein-Friesian female cattle infected *in vitro* with *M. bovis* and BCG across an infection time course of 2, 6, and 24 h. Differentially expressed (DE) genes were identified using established methods and a list of interactions were generated using the InnateDB resource (23) and visualized as interaction networks with the Cytoscape software package (24).

#### **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 (3-year-old) unrelated Holstein-Friesian females were used in this study. All animals were maintained under uniform housing conditions and nutritional regimens at the UCD Lyons Research Farm (Newcastle, County Kildare, Ireland). The animals were selected from an experimental herd without a recent history of BTB infection and all animals tested negative for the single intradermal tuberculin test (SICTT). These cattle were also negative for infection with *Brucella abortus*, *M. avium* subsp. *paratuberculosis*, *Salmonella Typhimurium*, bovine herpesvirus 1 (BHV-1), and bovine viral diarrhea (BVD) virus.

#### **MONOCYTE EXTRACTION AND CULTURE OF MDM**

Monocyte extraction and culture of bovine MDM was performed as described by us previously (20). Briefly, for monocyte isolation, 300 ml of whole blood was collected in acid citrate dextrose buffer (Sigma-Aldrich Ireland Ltd., Dublin, Ireland) in sterile bottles. Blood was layered onto Accuspin™ tubes containing Histopaque® 1077 (Sigma-Aldrich Ireland Ltd., Dublin, Ireland), and following density gradient centrifugation (500 g for 20 min) performed at room temperature, peripheral blood mononuclear cells (PBMC) were collected. Contaminating red blood cells (RBC) were removed following resuspension and subsequent incubation of the PBMC in RBC lysis buffer (10 mM KHCO3,150 mM NH4Cl, 0.1 mM EDTA pH 8.0) for 5 min at room temperature. After incubation, PBMC were washed twice with sterile phosphate-buffered saline (PBS; Invitrogen™, Life Technologies Corporation, Paisley, UK) before resuspending cells in PBS containing 1% bovine serum albumin (BSA; Sigma-Aldrich Ireland Ltd., Dublin, Ireland). Monocytes were then isolated using the MACS® protocol and MACSH MicroBeads conjugated to mouse anti-human CD14 antibodies (Miltenyi Biotec Ltd., Surrey, UK), which has been shown to be cross-reactive with bovine monocytes (25). The MACS® protocol was performed according to the manufacturer's instructions.

The identity and purity of monocytes was confirmed by flow cytometry using an anti-CD14 fluorescein-labeled antibody (data not shown). This method has been previously shown by us to yield a purity of CD14<sup>+</sup> cells ≥99% (19). Purified monocytes were seeded at 1 × 10<sup>6</sup> per well on 24-well tissue culture plates in RPMI 1640 medium (Invitrogen™, Life Technologies Corporation, Paisley, UK) containing 15% heat inactivated fetal calf serum (FCS; Sigma-Aldrich Ireland Ltd., Dublin, Ireland), 1% non-essential amino acids (NEAA; Sigma-Aldrich Ireland Ltd., Dublin, Ireland), gentamicin (5 mg/ml; Sigma-Aldrich Ireland Ltd., Dublin, Ireland), and incubated at 37°C, 5% CO2. Following 24 h incubation (day one), the media was replaced with 1 ml fresh antibiotic-containing media to remove any non-adhered cells. On day three, media was replaced with 1 ml antibiotic-free culture media (RPMI 1640 medium containing 15% heat inactivated FCS and 1% NEAA only). To ensure that the same number of MDM were subjected to mycobacterial challenge, cells were dissociated on day five using 1× non-enzymatic cell dissociation solution (Sigma-Aldrich Ireland Ltd., Dublin, Ireland), counted, and then re-seeded at 2 × 10<sup>5</sup> cells per well in 24-well tissue culture plates (Sarstedt Ltd., County Wexford, Ireland) using antibiotic-free culture media. By day eight, 80–100% confluent monolayers of MDM

were generated that displayed the characteristic macrophage morphology as confirmed by Giemsa staining (data not shown). On day eight, MDM were used for the *in vitro* challenge experiments with *M. bovis* and *M. bovis* BCG.

#### **CULTURE OF M. BOVIS AND BCG AND INFECTION OF BOVINE MDM**

The culturing of *M. bovis* (strain M2137; spoligotype SB0142) has been described by us previously (20). Culturing of BCG (Pasteur strain) was performed in a Biosafety Containment Level 3 (CL3) laboratory and conformed to the national guidelines on the use of Hazard Group 3 infectious organisms. BCG stocks for the *in vitro* MDM infection experiments were cultured in Middlebrook 7H9 medium (Difco™, Becton, Dickinson Ltd., Oxford, UK) containing 10% (vol/vol) BBL™ Middlebrook OADC Enrichment (Difco™, Becton, Dickinson Ltd., Oxford, UK), 0.05% Tween 80 (Sigma-Aldrich, Dublin, Ireland), and 0.50% (weight/vol) glycerol (Sigma-Aldrich Ltd., Dublin, Ireland) at 37°C. Bacterial cultures were grown to mid-logarithmic phase as determined by spectrophotometric analysis prior to the challenge experiments using the purified bovine MDM. Colonyforming unit (cfu) counting was performed using Middlebrook 7H11 medium (Difco™, Becton-Dickinson Ltd., Dublin, Ireland) containing 10% (vol/vol) Middlebrook (ADC) enrichment (Difco™, Becton-Dickinson Ltd., Dublin, Ireland) and 0.50% (vol/vol) glycerol (Sigma-Aldrich, Dublin, Ireland).

Monocyte-derived macrophage infections with BCG were performed as previously described for MDM infections with *M. bovis* (20). In brief, MDM (seeded at 2 × 10<sup>5</sup> cells per well) were challenged with BCG (4 × 10<sup>5</sup> cells per well) prepared in antibiotic-free culture media (RPMI 1640 medium containing 15% heat inactivated FCS and 1% NEAA only). BCG cell counts were performed using a Petroff Hausser chamber (Fisher Scientific Ltd., Dublin, Ireland) following sterile filtering to prevent clumping using a 5µm filter (Millipore Ireland Ltd., County Cork, Ireland). BCG cell numbers were then adjusted and 4 × 10<sup>5</sup> bacilli were added to the appropriate wells giving a multiplicity of infection (MOI) of 2:1. Subsequent cfu for BCG counting also yielded a mean MOI of 2:1. Once challenged, MDM were incubated at 37°C, 5% CO<sup>2</sup> for 2, 6, and 24 h.

#### **RNA EXTRACTION AND MICROARRAY ANALYSIS**

RNA extraction and microarray analysis has been described by us previously (20). Briefly, global gene expression was analyzed using the pan-genomic high-density Affymetrix® GeneChip® Bovine Genome Array (Affymetrix UK Ltd., High Wycombe, UK<sup>1</sup> ). This array contains 24,072 probe sets representing over 23,000 gene transcripts and includes approximately 19,000 Uni-Gene clusters. cDNA labeling, hybridization, and scanning for the microarray experiments were performed by Almac Diagnostics Ltd. (Craigavon, Co. Armagh, Northern Ireland) using a one-cycle amplification/labeling protocol.

#### **STATISTICAL ANALYSIS OF MICROARRAY DATA**

Affymetrix® GeneChip® Bovine Genome Array data were analyzed using BioConductor (26) 2 contained within the R statistical package<sup>3</sup> . All raw data were normalized using the factor analysis for robust microarray summarization (FARMS) package. The FARMS package uses only perfect match (PM) probes and a quantile normalization procedure, providing both *P*-values and signal intensities (27). Normalized data were then further subjected to filtering for informative probes sets using the I/NI-calls package in R (28). This defines a probe set as being informative when many of its probes reflect the same change in mRNA concentration across arrays. DE genes for a paired comparison between *M. bovis* versus BCG infections at each time point were extracted using a moderated paired *t*-test within the Linear Models for Microarray Data (LIMMA) R package (29) 4 . Genes displaying differential expression were annotated using the Affymetrix® bovine gene annotation. The Benjamini–Hochberg multiple-testing correction method (30) was applied to all DE genes to minimize the false discovery rate (FDR) and adjusted *P*-values for all DE genes were calculated. Only DE genes with an adjusted *P*-value of ≤0.01 were used for subsequent network reconstruction and analysis. All data are MIAME compliant and have been submitted to the NCBI gene expression omnibus (GEO) database with experiment series accession number GSE59774.

#### **GENERATION OF InnateDB INTERACTION LISTS**

Affymetrix® probe IDs for the DE gene lists at each time point were mapped to human Ensembl IDs, using the BioMart search tool on the Ensembl genome database<sup>5</sup> . Each of the three DE gene lists were uploaded separately to the InnateDB<sup>6</sup> curated, interaction database (23), where a list of interactions between the uploaded molecules were generated. At the time of this analysis, the InnateDB resource contained over 195,000 experimentally determined interactions, of which more than 18,000 have been manually curated, allowing for a detailed systems-level analysis of the innate immune response. The interactions, for each time point, were then viewed as a reconstructed network using the Cytoscape software (24) 7 . Small clusters of nodes that did not connect to the main network were removed from the analysis.

#### **24 h INTERACTION NETWORK ANALYSIS**

The 24-h interaction network was analyzed as an undirected, Markov, network. Statistical analysis of the network was performed using the *NetworkAnalyzer* plugin of Cytoscape (31). Hub nodes were identified by calculating the degree of connectivity (DOC) for each node in the network using the *Degree Sorted Circle Layout* in Cytoscape. DOC is a measure of the number of interactions a node has within the network. Nodes with the highest DOC are classed as hub nodes. Bottleneck nodes were identified by calculating the betweenness centrality index (BCI) for each node in the network, also using the *NetworkAnalyzer* plugin. BCI is a measure of how often a node appears on the shortest path between nodes in the network and is a reflection of the amount of control a node exerts over the interactions of other nodes in the

<sup>1</sup>www.affymetrix.com

<sup>2</sup>www.bioconductor.org

<sup>3</sup>www.r-project.org

<sup>4</sup>www.bioconductor.org/packages/release/bioc/html/limma.html

<sup>5</sup>www.ensembl.org/biomart/martview

<sup>6</sup>www.innatedb.com

<sup>7</sup>www.cytoscape.org

network (32). A Spearman rank correlation between the DOC and log<sup>2</sup> fold-change in gene expression (*M. bovis*-infected MDM relative to BCG-infected MDM) was performed in the SPSS statistics version 20 package (IBM, New York, NY, USA).

#### **Combination of interactions and differentially expressed gene lists**

Several Perl scripts were implemented to combine gene lists and interactions in various ways. One was used to calculate the numbers of genes that are DE and linked to a set of hub genes by a certain distance of edges/interactions, another one to extract expression values for interactors with bottleneck genes at all time points. All scripts are available through GitHub at: https: //github.com/kkillick/Systems-biology-paper.git.

#### **Gene ontology of 24 h network**

Ingenuity® systems pathway analysis (IPA) was used to identify gene ontology (GO) functional categories over-represented within the DE genes found in the 24-h interaction network. 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. For IPA GO analysis, the Ingenuity® knowledge base was used as a reference set and only genes within the 24-h interaction network were uploaded. IPA then performed an over-representation analysis that categorized the DE genes within the uploaded list into functional GO categories. Each GO category in IPA is ranked based on the number of DE genes falling into each functional group. Right-tailed Fisher's exact tests were used to calculate an overlap *P*-value for each of the biological functions assigned to the list of DE genes.

#### **Expression analysis of the interactors of the top hub and bottleneck genes**

The top hub and bottleneck genes (*IKBKE*, *MYC*, *NFKB1*, and *EGR1*) were specified in InnateDB to generate lists of all human interactions. DE genes for a paired comparison of *M. bovis* versus BCG infections at each time point, using an adjusted *P*-value of ≤0.01, were extracted from the complete interaction lists and visualized using the Circos software package<sup>8</sup> (33). Interaction nodes colored in red indicate that the gene was up-regulated within the dataset, while green indicate that the gene was down-regulated.

#### **Identification of scale-free and small-world properties within the 24-h network**

If a network has a node degree distribution that can be fitted with a power law distribution, it can be a sign that the network has a scale-free architecture. Node degree distribution of the 24 h network was calculated using the *NetworkAnalyzer* plugin of cytoscape and a power law distribution was fitted using Microsoft Excel.

To determine whether a network displays small-world properties it must show *L*>˜ *L*random but *C C*random. Where *L* is the characteristic path length, defined as the number of edges in the

<sup>8</sup>http://circos.ca

shortest path between two nodes, averaged over all pairs of nodes, and *C*, the clustering coefficient, is defined as follows: node *n* has *k<sup>n</sup>* max number of edges. *C<sup>n</sup>* denotes the fraction of these edges that actually exist, while *C* is defined at the average of *C<sup>n</sup>* over all *n* (34).

The R package igraph<sup>9</sup> was used to generate 1,000 random networks, each with 716 nodes and 1,785 edges. The mean characteristic path length and the mean clustering coefficient were calculated across the 1,000 networks also using the average path length and transitivity functions of the igraph R package. All R code used here is available through GitHub at: https://github.com/ kkillick/Systems-biology-paper.git.

#### **REAL TIME QUANTITATIVE REVERSE TRANSCRIPTION (qRT)-PCR VALIDATION OF MICROARRAY RESULTS**

A panel of 12 genes that were DE across all three time points based on the microarray data was selected for microarray validation using real time quantitative reverse transcription PCR (qRT-PCR) analysis. The laboratory and statistical methods used to perform this analysis are described in Section "Real Time Quantitative Reverse Transcription (qRT)-PCR Validation of Microarray Results" in Supplementary Material, and Table S1 in Supplementary Material details the qRT-PCR primers used in this study.

#### **RESULTS**

#### **INCREASING DIFFERENTIAL GENE EXPRESSION BETWEEN M. BOVIS-VERSUS BCG-INFECTED MDM OVER TIME**

Comparison of gene expression profiles between paired *M. bovis* versus BCG *in vitro* infections revealed 702 DE genes at 2 h and 1,000 and 2,674 DE genes at 6 and 24 h post-infection, respectively (adjusted *P*-value ≤0.01). Of the 702 DE genes found at the 2-h time point, 334 genes displayed increased expression in *M. bovis* infections relative to BCG infections (hereafter referred to as up-regulated), while 368 genes showed decreased expression in *M. bovis* infections relative to BCG infections (hereafter referred to as down-regulated). At the 6-h time point, 569 genes were upregulated, with 431 genes displaying down-regulation, while at the 24-h time point 1,152 and 1,522 genes were up- and downregulated, respectively. Tables S2–S4 in Supplementary Material contain DE gene results for the 2, 6, and 24 h post-infection time points, respectively. **Figure 1** depicts the number of DE genes and the relative change in expression across the infection time course.

#### **InnateDB INTERACTION ANALYSES**

InnateDB is an open-source database containing experimentally verified gene and protein-interactions involved in the human, mouse, and bovine innate immune response (23). As a larger number of curated interactions are available for human data, the DE genes detected here were mapped to corresponding human Ensembl gene IDs prior to analysis with InnateDB. In total, 362 DE genes at 2 h, 554 DE genes at 6 h, and 1,503 DE genes at 24 h mapped to human Ensembl gene IDs. These gene lists were examined separately using InnateDB to generate a catalog of interactions for each infection time point. At 2 h post-infection, 205

<sup>9</sup>http://igraph.sourceforge.net/doc/R/igraph.pdf

interactions were identified, while 413 and 1,785 interactions were detected for 6 and 24 h, respectively. As the 24-h time point generated the largest interaction network, it was chosen for detailed downstream analyses.

#### **24 h INTERACTION NETWORK RECONSTRUCTION AND ANALYSIS**

Removal of small groups of nodes not linked to the main cluster left 615 nodes for analysis in the 24-h interaction network. These were connected via 1,670 edges. Two hundred eight of these nodes had self-interaction loops. The network diameter was 10 and the average shortest path length was 4.315. The *IKBKE*, *MYC*, *NFKB1*, *HDAC5*, and *TRAF2* genes were identified as the top hub nodes within the network. A degree of connectivity (DOC) of 57 was found for *IKBKE*, 57 for *MYC*, 55 for *NFKB1*, 52 for *HDAC5*, and 50 for *TRAF2*. The DOC for these five nodes was markedly higher than the mean number of connections found for nodes within the network of 3.857, suggesting that they are hub genes within the 24-h network. Also, for the 24-h interaction network, *IKBKE*, *MYC*, and *EGR1* displayed the highest betweenness centrality index (BCI) of 0.1661, 0.1360, and 0.1297, respectively (Figure S1 in Supplementary Material). This finding suggests that *IKBKE*, *MYC*, and *EGR1* are also key bottleneck nodes within the 24-h interaction network.

No significant correlation between the DOC and the log<sup>2</sup> foldchange in gene expression for nodes was found, indicating that nodes with the greatest changes in expression between groups may not have the biggest influence on other nodes within the network (Spearman rho value 0.066,*P* = 0.108) (**Figure 2**). Figures S2 and S3 in Supplementary Material show the biological interaction networks for the 2 and 6-h infection time points, respectively.

#### **SCALE-FREE AND SMALL-WORLD PROPERTIES OF THE M. BOVIS VERSUS BCG 24 h INTERACTION NETWORK**

A power law distribution was fitted to the node degree distribution of the 24-h network, generating an *R* 2 value of 0.867 (Figure S4 in Supplementary Material). This result suggests that the 24-h interaction network has a scale-free architecture, whereby many nodes have a small number of connections yet a small number of hub nodes have many connections. Scale-free networks tend to exhibit the small-world network property of high clustering and short average path lengths, *L*>˜ *Lrandom* but *C Crandom*, (where *L* is the characteristic path length and *C* is the clustering coefficient) (34). This was tested on the 24-h interaction network by comparing it to 1,000 random networks generated with the igraph *R* package using the same number of nodes and edges. The random networks yielded a mean *C* = 0.007, with a maximum *C* = 0.012 and *L* = 4.251. Conversely, for the observed 24 h interaction network, the *C* and *L* values were 0.112 and 4.315, respectively. With a *C* value 15.9-fold higher than the mean *C* value for the 1,000 random networks, and the *L* value only marginally higher, the 24-h interaction network exhibits small-world properties.

Highlighting the gene neighborhoods at increasing distances from a main hub node (*MYC*) in the 24-h interaction network demonstrates the scale-free architecture of the network. **Figure 3** shows that almost all other nodes in the network (587 out of 615) are connected to *MYC* through less than five edges. The combined effect of the three main hub gene nodes (*IKBKE*, *MYC*, and *NFKB1*) is presented in **Figure 4A**. It shows that nearly all DE genes in the network (96%) are linked through a maximum of four edges to one of the three key hub genes. This is significantly higher than expected by chance (*P*-value < 0.01), based on

100,000 permutation tests. In addition, the intersection between the linked genes from all hubs rises quickly from less than 10% at a distance of two edges to nearly 90% at a distance of four edges (**Figure 4B**), indicating complex intertwinement and a capacity for rapid dissemination of signals throughout the network.

#### **GENE ONTOLOGY OF THE M. BOVIS VERSUS BCG 24 h INTERACTION NETWORK**

The total number of DE genes that could be mapped to molecules in the Ingenuity® Knowledge Base was 507, out of the 615 genes in the uploaded dataset. GO analysis facilitates discovery of molecular functions that are enriched within a dataset. IPA GO analysis of molecules within the 24-h interaction network revealed *Cell death* to be a major over-represented functional category. Two hundred thirty-seven of the molecules within the 24-h network were involved in *Cell death* and an adjusted *P*-value of 4.18 × 10−<sup>25</sup> was obtained for this category. Further inspection showed that the top 15 hub nodes and the top 15 bottleneck nodes were also involved in the *Cell death* GO category.

#### **EXPRESSION ANALYSIS OF THE INTERACTORS OF THE TOP HUB AND BOTTLENECK GENES**

Expression changes of the interactors with the top hub and bottleneck genes *IKBKE*, *NFKB1*, *MYC*, and *EGR1* were examined at each of the three time points. Using an adjusted *P*-value of ≤0.01, it was observed that the majority of the interactors were DE at 24 h for a paired comparison of *M. bovis* versus BCG infections. This is in contrast to the earlier time points where only small numbers of interactors exhibited significant expression changes (adjusted *P*-value ≤0.01). For example, the number of DE genes that interact with *IKBKE* increased from 6 to 12 between time points 2 and 6 h, but expanded to 68 at 24 h. These results suggest that differential regulatory activity at the *IKBKE* hub gene and its interactors increases over time in MDM infected with *M. bovis* relative to BCG. A similar outcome was observed for the *NFKB1*

and *MYC* hub genes and the *EGR1* hub/bottleneck gene, where the equivalent numbers of DE genes for the 2, 6, and 24-h time points were 19, 27, and 42; 15, 22, and 72; and 23, 26 and 54, respectively (**Figure 5**).

#### **REAL TIME QUANTITATIVE REVERSE TRANSCRIPTION PCR (qRT-PCR) VALIDATION OF MICROARRAY RESULTS**

The results from real time qRT-PCR validation of the microarray results are provided in Table S5 in Supplementary Material.

#### **DISCUSSION**

The ability of pathogenic mycobacteria to cause disease has been largely attributed to their capacity to survive and replicate in macrophages via a diverse range of mechanisms that subvert the host innate immune response (35). It is believed that mycobacterial virulence is, in part, governed by several secreted virulence factors, such as members of the ESAT-6 protein family encoded by the RD1 genome region. Notably, the RD1 region is present in all virulent *M. tuberculosis* and *M. bovis* strains and is one of several regions of difference deleted from the genome of all strains of attenuated BCG (36). The presence or absence of these virulence factors are thought to result in differential host innate immune responses to virulent and attenuated mycobacteria that can be identified through the analysis of host gene expression data (37–39). In the present study, we have analyzed the global differential bovine macrophage transcriptional response to virulent *M. bovis* and attenuated BCG across an *in vitro* infection time course of 2, 6, and 24 h. Furthermore, we have applied a systems biology approach to reconstruct the complex molecular interactions and networks that illustrate the differential macrophage responses to these mycobacterial strains. Further work will be necessary to provide insights into whether key genes identified in this study represent aspects of attenuation regulated by the RD1 region.

Comparison of the MDM response to infection with *M. bovis* relative to BCG revealed the largest number of DE genes at

24 h post-infection (2,674 DE genes), which was 3.8- and 2.7 fold higher than the number of DE genes at the observed at 2 h (702 DE genes) and 6 h time points (1,000 DE genes), respectively (**Figure 1**). Notably, comparison between *M. bovis*- and BCG-infected MDM at the 2-h time point did not reveal significant differences in the expression profiles of key immune-related genes known to be important during mycobacterial infection. For example, mRNA transcripts were detected for genes encoding macrophage pattern recognition receptors (PRRs) such as *TLR2* and *TLR4*, which are involved in the recognition of evolutionarily conserved pathogen-associated molecular patterns (PAMPs), and inflammatory cytokines including *IL1B*, *IL6*, and

*IL10* following infection; however, these were not DE between treatments (adjusted *P*-value of 0.01).

Systems analysis of the DE genes identified at the 24-h time point generated a biological interaction network that displayed both scale-free and small-world properties (**Figure 3**). In a scalefree network, nodes are highly clustered and can be connected to each other by relatively few steps due, in part, to the presence of a small number of highly connected hub nodes (34, 40). Hub nodes, identified as those nodes that have the highest DOC within a network, enable scale-free networks to be highly robust. The majority of nodes in the network are non-hub nodes with only a small

number of connections; therefore, scale-free networks can withstand the removal of non-hub nodes and still retain the capacity to transmit information across the whole network. As such, these networks are said to be "error tolerant." There are many examples of both small-world and scale-free networks, including metabolic networks (41–43), neural networks in the brain (44, 45), and the Internet (46). However, scale-free networks are vulnerable: should key hub nodes be removed, the network will be seriously compromised. Therefore, although scale-free and small-world networks are highly robust against random removal of non-hub nodes, they are vulnerable to targeted attacks on key hub nodes (47).

In the current study, *IKBKE*, *MYC*, and *NFKB1* were the top three ranking hub genes within the 24-h network (**Figure 2**). The importance of these genes is highlighted by the effect they have on the network based on the number of DE genes that are close neighbors (up to four edges away, **Figure 4**). There is not enough information available from the interaction databases to infer directionality for all interactions in the network, but the type of products associated with these hubs (inhibitors and activators, see below) suggest that most of their interactors lie downstream within a signaling cascade. The fact that nearly all of the DE genes in the network are linked by a maximum of four edges to one of the three main hub genes emphasizes the rapid and extensive response that can be initiated by change of expression in a very small number of genes. Furthermore, although the first- and second-degree DE neighbors of the hubs are mostly separate sub-groups, these sets of genes merge quickly into one large group after the next two extensions of the network. In other words, the overall number of DE genes (but not necessarily the type of expression change) could potentially be derived from changes at any one of the hubs, which points toward a single response captured by the network, compared to a group of separate responses. This correlates well with the detection of one common category (*Cell death*) found to be significantly over-represented in the GO analysis; however, observable differences in macrophage viability or morphology between BCG and *M. bovis* infections were not examined. The identification of three hubs instead of just one could be indicative of (i) redundancy, which has evolved for more robustness in the response to infections; or (ii) fine-tuning derived from interplay between different regulators, which could provide greater control of the response.

The products of the *IKBKE* and *NFKB1* genes are central to the NF-κB signaling pathway that regulates the expression of many immune-related genes following infection. *IKBKE* encodes a subunit of the inhibitor of κB (IκB) kinase (IKK) complex, which phosphorylates the NF-κB inhibitors (IKBs), resulting in the release of the NF-κB transcription factor complex, a subunit of which is encoded by the *NFKB1* gene (48). Notably, at the 24 h time point, both *IKBKE* and *NFKB1* were up-regulated in the *M. bovis*-infected MDM relative to BCG-infected MDM, suggesting that there is increased activation of NF-κB signaling following infection with the virulent mycobacterial strain. This is further supported by the increased relative expression of several NF-κBinducible genes known to play a role in the host response to *M. bovis*, such as *IL1A*, *IL1B*, *IL6*, and *TNF* (49–51).

The up-regulation of the *IKBKE* and *NFKB1* hub genes indicates that activation of NF-κB and consequent downstream signaling and interactors of these hub genes (**Figure 5**) are crucial for the differential success of virulent and attenuated *M. bovis* strains within the host macrophage. Indeed, several pathogens, including *Helicobacter pylori* and *Chlamydia pneumoniae*, have been shown to manipulate the host NF-κB signaling pathway to ensure their survival (52). Evidence of the modulation of monocyte NF-κB signaling following mycobacterial infection has also been demonstrated. Dhiman et al. (53) showed that human monocytes infected with virulent *M. tuberculosis* exhibited increased NF-κB activity that resulted in up-regulation of NF-κB-inducible anti-apoptotic genes and increased survival relative to monocytes infected with avirulent *M. tuberculosis*. In the current study, we also observed increased expression of several NF-κB-inducible anti-apoptotic genes including *BIRC3*, *CFLAR*, and *TRAF2*, which were among the downstream interactors of the *IKBKE* and *NFKB1* hub genes (**Figure 5**). Furthermore, GO analysis of all nodes in the 24-h network revealed enrichment for genes involved in apoptosis, which is increasingly regarded as a host innate immune mechanism that contains and limits mycobacterial growth following infection (54– 56). However, the induction of anti-apoptotic genes may represent an immuno-modulation strategy employed by the pathogen that inhibits apoptosis, enabling survival, and replication within the macrophage (55). Therefore, it is possible that virulent *M. bovis* targets key members of the NF-kB pathway causing the induction of anti-apoptotic mechanisms that result in prolonged viability within the bovine macrophage.

*MYC*, another highly connected hub gene identified within the 24-h network, displayed the same DOC as *IKBKE*. *MYC* encodes a multi-functional transcription factor that serves as a key regulator of cell proliferation and apoptosis and has been shown to suppress or activate the expression of its target genes (57, 58). A recent study has shown that *MYC* expression was induced in primary human blood-derived macrophages following *in vitro* infection with pathogenic and non-pathogenic mycobacterial species. Moreover, it was demonstrated that the MYC transcription factor mediated the suppression of intra-macrophage mycobacterial growth via the activation of cytokines including TNF and IL-6 (59). In the current study, *MYC* was down-regulated in *M. bovis*-infected MDM relative to BCG-infected MDM within the 24-h network. Furthermore, 51 of 73 (69.8%) of the downstream *MYC* interactors at the 24-h time point were also down-regulated (**Figure 5**). Previous work by our group using the same *M. bovis*-infected MDM described here demonstrated that *MYC* was down-regulated relative to non-infected MDM from the same animals (20, 22). It is also important to note that *MYC* was not DE in the BCG-infected MDM analyzed here relative to the same control MDM (unpublished results). Collectively, these results suggest that suppression of *MYC* following *M. bovis* infection results in the suppression of host pro-inflammatory responses leading to intracellular microbial survival and growth. Hence, *MYC* may serve as a key target for mycobacterial modulation of innate immune mechanisms, facilitating persistence within host macrophages.

Bottleneck nodes in the 24-h network were identified as nodes with a high BCI (Figure S1 in Supplementary Material). BCI has been shown to be associated with constrained evolutionary rates for gene or protein nodes in biological networks; studies involving eukaryote protein-interaction networks have demonstrated that bottleneck nodes evolve more slowly than non-bottleneck nodes (60, 61). Analyses of protein-interaction networks have also shown that bottleneck nodes are more likely to be essential for organismal survival than non-bottleneck nodes (62, 63).

In the current study, *EGR1*, together with *IKBKE* and *MYC*, was identified as a top bottleneck node within the 24-h interaction network (Figure S1 in Supplementary Material), with the majority of significantly DE genes interacting with *EGR1* showing changes at this time point (55 out of 85,**Figure 5**). *EGR1* encodes a zinc-finger transcription factor that regulates the expression of a large number of genes involved in cellular differentiation and mitogenesis. In addition, EGR-1 has also been shown to regulate the transcription of pro-inflammatory cytokines in murine macrophages in response to stimulation with bacterial antigens, including TNF and IL-6 (64, 65). To our knowledge, only one publication exists associating *EGR1* with mycobacteria–macrophage interactions (66), although it has recently been identified as part of a conserved macrophage core transcriptional response module that is DE in response to intracellular bacterial pathogens (67). Here, identification of *EGR1* as a key bottleneck gene suggests that this gene is a component of a differential bovine macrophage response module to virulent and attenuated mycobacterial strains.

It has been shown that many of the species and subspecies in the *M. tuberculosis* complex exhibit specific host association (68). It is believed that this host specificity is an evolutionary consequence of reciprocal adaptive genetic changes that have occurred due to strong selection pressures between the interacting pathogenic mycobacteria and their hosts (69, 70). Several examples of reciprocal traits involved in host-mycobacterial pathogen coevolution have been described, including the ability of the host immune system to clear infection versus the ability of virulent pathogens to subvert and evade host immune responses (71). It is likely that the small-world properties of the 24-h network reflect host–pathogen co-evolution. Small-world networks are characterized by the presence of hub and bottleneck genes that interact with all other nodes within the network via a small number of steps. This model predicts that the expression and products of these hub and bottleneck genes enable the host to regulate an appropriate and efficient immune response to virulent or attenuated mycobacterial strains. Alternatively, key hub and bottleneck genes may serve as targets for the subversion of host immune responses by virulent pathogens that underlie successful infection, such as the suppression of apoptosis-related mechanisms by host transcription factors as demonstrated by Dyer et al. (72).

At the individual gene level, from the results presented here, it is not possible to determine whether host macrophage gene expression changes are due to regulation of host gene networks by the pathogen, the host, or a combination of both. In this regard, previous studies have demonstrated that pathogen-encoded protein virulence factors and, more recently, pathogen-encoded small regulatory RNA molecules manipulate and modulate host macrophage responses to facilitate intracellular survival and dissemination of pathogenic mycobacteria (73–79). Consequently, to fully understand the complex transcriptional regulatory interplay underlying host–mycobacterium interactions, it will be necessary to perform parallel high-throughput transcriptional profiling of mRNA and microRNA in both host and microbial cells (80, 81) and complement these with experimental detection of host–pathogen molecular interactions (82, 83).

#### **ACKNOWLEDGMENTS**

This work was supported by Investigator Grants from Science Foundation Ireland (Nos: 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 PhD program10. We would like to thank Mr. Eamon Costello and staff at the Irish Department of Agriculture, Food and the Marine Backweston Laboratory Campus (Celbridge Co., Kildare) for providing the isolate of *M. bovis*. We also would like to thank the staff at the

UCD Lyons Research Farm for assistance with blood collection and we thank Prof. Joseph Keane (TCD) for valuable comments.

#### **SUPPLEMENTARY MATERIAL**

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

#### **REFERENCES**


<sup>10</sup>http://bioinfo-casl.ucd.ie/PhD

*Mycobacterium bovis* isolates. *Biomed Res Int* (2013) **2013**:458278. doi:10.1155/ 2013/458278


**Conflict of Interest Statement:** The Guest Associate Editor Kieran Meade declares that, despite having collaborated with the authors, Kate E. Killick, David A. Magee, Stephen D. E. Park, John A. Browne, Nicolas C. Nalpas, Eamonn Gormley, Stephen V. Gordon, David E. MacHugh, Karsten Hokamp the review process was handled objectively. 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.

*Received: 23 May 2014; paper pending published: 30 June 2014; accepted: 19 August 2014; published online: 01 October 2014.*

*Citation: Killick KE, Magee DA, Park SDE, Taraktsoglou M, Browne JA, Conlon KM, Nalpas NC, Gormley E, Gordon SV, MacHugh DE and Hokamp K (2014) Key hub and bottleneck genes differentiate the macrophage response to virulent and attenuated Mycobacterium bovis. Front. Immunol. 5:422. doi: 10.3389/fimmu.2014.00422*

*This article was submitted to Molecular Innate Immunity, a section of the journal Frontiers in Immunology.*

*Copyright © 2014 Killick, Magee, Park, Taraktsoglou, Browne, Conlon, Nalpas, Gormley, Gordon, MacHugh and Hokamp. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

## Comparative functional genomics and the bovine macrophage response to strains of the Mycobacterium genus

#### **Kévin Rue-Albrecht <sup>1</sup> , David A. Magee<sup>1</sup>† , Kate E. Killick 1,2, Nicolas C. Nalpas <sup>1</sup> , Stephen V. Gordon3,4 and David E. MacHugh1,4\***

<sup>1</sup> Animal Genomics Laboratory, UCD School of Agriculture and Food Science, University College Dublin, Dublin, Ireland

<sup>2</sup> Systems Biology Ireland, UCD Conway Institute of Biomolecular and Biomedical Research, University College Dublin, Dublin, Ireland

<sup>3</sup> UCD School of Veterinary Medicine, University College Dublin, Dublin, Ireland

<sup>4</sup> UCD Conway Institute of Biomolecular and Biomedical Research, University College Dublin, Dublin, Ireland

#### **Edited by:**

Kieran G. Meade, Teagasc, Ireland

#### **Reviewed by:**

Paul M. Coussens, Michigan State University, USA Annapurna Nayak, Brunel University, UK

#### **\*Correspondence:**

David E. MacHugh, Animal Genomics Laboratory, UCD School of Agriculture and Food Science, University College Dublin, Belfield, Dublin 4, Ireland e-mail: david.machugh@ucd.ie

#### **†Present address:**

David A. Magee, Department of Animal Science, University of Connecticut, Storrs, CT, USA

Mycobacterial infections are major causes of morbidity and mortality in cattle and are also potential zoonotic agents with implications for human health. Despite the implementation of comprehensive animal surveillance programs, many mycobacterial diseases have remained recalcitrant to eradication in several industrialized countries.Two major mycobacterial pathogens of cattle are Mycobacterium bovis and Mycobacterium avium subspecies paratuberculosis (MAP), the causative agents of bovine tuberculosis (BTB) and Johne's disease (JD), respectively. BTB is a chronic, granulomatous disease of the respiratory tract that is spread via aerosol transmission, while JD is a chronic granulomatous disease of the intestines that is transmitted via the fecal-oral route. Although these diseases exhibit differential tissue tropism and distinct complex etiologies, both M. bovis and MAP infect, reside, and replicate in host macrophages – the key host innate immune cell that encounters mycobacterial pathogens after initial exposure and mediates the subsequent immune response.The persistence of M. bovis and MAP in macrophages relies on a diverse series of immunomodulatory mechanisms, including the inhibition of phagosome maturation and apoptosis, generation of cytokine-induced necrosis enabling dissemination of infection through the host, local pathology, and ultimately shedding of the pathogen. Here, we review the bovine macrophage response to infection with M. bovis and MAP. In particular, we describe how recent advances in functional genomics are shedding light on the host macrophage–pathogen interactions that underlie different mycobacterial diseases. To illustrate this, we present new analyses of previously published bovine macrophage transcriptomics data following in vitro infection with virulent M. bovis, the attenuated vaccine strain M. bovis BCG, and MAP, and discuss our findings with respect to the differing etiologies of BTB and JD.

**Keywords: cattle, BCG, gene expression, Johne's disease, macrophage, Mycobacterium avium subspecies paratuberculosis, Mycobacterium bovis, tuberculosis**

#### **INTRODUCTION**

*Mycobacterium* is a Gram-positive genus of Actinobacteria that includes more than 120 species (1, 2). Although the majority of species in this genus are non-pathogenic environmental bacteria, a few species are highly successful intracellular pathogens of human beings and other mammals including *Mycobacterium tuberculosis* and *Mycobacterium bovis –* the causative agents of human being and bovine tuberculosis (BTB), respectively – and *Mycobacterium avium* subspecies *paratuberculosis* (MAP), the causative agent of Johne's disease (JD) in cattle (3, 4). The success of these pathogenic mycobacteria is partly due to their ability to infect, reside, and proliferate inside host macrophages, despite the antimicrobial properties of these cells. Macrophages serve as key effector innate immune cells that mediate the initial host response to infection via the activity of inflammatory cytokines and chemokines; this initial interaction leads to either the eradication of intracellular

bacilli or the formation of organized collections of immune cells, termed granulomas, which contain infection (5).

Infections with pathogenic mycobacteria can manifest as acute or chronic disease or involve lengthy subclinical phases of infection with the potential to reactivate later. It is also understood that the establishment of successful infection is underpinned by subversion and modulation of host macrophage antimicrobial mechanisms, including the prevention of macrophage phagosome– lysosome fusion, inhibition of macrophage apoptosis, and suppression of antigen presentation and signaling mechanisms within the macrophage (6–8). Furthermore, it has been proposed that virulent mycobacteria exploit host defense mechanisms, such as the induction of cytokine-induced necrosis, which results in immunopathology, the dissemination of infection through the host and ultimately pathology that leads to shedding of the pathogen from the host, thereby maintaining the cycle of infection (9). Consequently, investigating the complex interplay between mycobacterial pathogens and the host macrophage is critical to our understanding of the immuno-pathogenesis of mycobacterial diseases.

#### **THE Mycobacterium tuberculosis COMPLEX**

The genus *Mycobacterium* contains the *Mycobacterium tuberculosis* complex (MTBC) that includes seven major pathogenic species and subspecies that cause tuberculosis in a range of mammalian hosts, the most well-studied member of which is *M. tuberculosis* – the causative agent of human tuberculosis. Typically, the members of the MTBC display greater than 99.95% nucleotide sequence identity at the genome level, with little or no evidence for the exchange of genetic material between species and strains (10). Despite this high level of genome similarity, the members of the MTBC differ with respect to host range and pathogenicity: *M. tuberculosis* and *Mycobacterium africanum* are almost exclusively human pathogens; *Mycobacterium microti* causes disease in rodents including voles; *Mycobacterium pinnipedii* causes tuberculosis in marine mammals including seals and sea lions; and *Mycobacterium caprae* is very closely related to *M. bovis* and infects both goats and deer. The species with the largest host range is *M. bovis*, which is mainly isolated from cattle, but can also be responsible for outbreaks in wild animals. Furthermore, *M. bovis* can cause disease in human beings yet rarely transmits between immunocompetent hosts. A closely related mycobacterial species, *Mycobacterium canettii*, causes pathology in human beings, but differs from the other members of the MTBC in that it displays smooth colony morphology rather than the characteristic rough morphology of the other MTBC members (11, 12).

Phylogenetic analyses using insertion/deletion DNA sequence polymorphisms (indels), such as the variable regions of difference (RD – see below) and whole gene and genome sequences have revealed that the evolutionary history of the MTBC represents a pattern of genome downsizing characterized by chromosomal DNA sequence deletions and the inability of these species to repair deletions through recombinogenic processes (10, 13). These studies support a distinct phylogenetic position of *M. canettii* from all other MTBC members: *M. canettii* strains possess intact RD sequences that are absent from the other MTBC species together with one species-specific deletion (RDcan). *M. canettii* strains also have 26 additional spacer sequences that are not found in other MTBC species [**Figure 1**] (10). Indeed, it has been recently proposed that *M. canettii* and other smooth tubercle bacilli (STB) lineages diverged from the common ancestor of all tubercle bacilli prior to the clonal radiation of non-smooth MTBC lineages

and that non-smooth MTBC lineages evolved from an STB-like mycobacterial ancestor, sometimes referred to as *Mycobacterium prototuberculosis* (14).

Genomic comparisons across the MTBC revealed a number of "regions of difference" (RD loci), with the presence or absence of these loci capable of differentiating the constituent strains of the MTBC. Notably, deletion of chromosomal region RD9 distinguishes *M. tuberculosis* strains from animal-adapted lineages, including *M. bovis*, while RD1 [which encodes, among other genes, the ESX-1 secretion system plus key secreted effectors including early secretory antigen target 6 (ESAT-6) and culture filtrate protein 10 (CFP-10)] is deleted in *M. bovis* Bacille Calmette–Guérin (*M. bovis* BCG) vaccine strains. The use of deletions to differentiate mycobacterial strains led to the proposal of an evolutionary scenario positing *M. tuberculosis* as being closer to the common ancestor of the MTBC and rejection of the hypothesis that *M. tuberculosis* infection in human populations arose from an animal pathogen such as *M. bovis* in parallel with cattle domestication and husbandry (15, 16).

#### **THE Mycobacterium avium COMPLEX**

A second major evolutionarily distinct cluster of mycobacterial species is represented by the *Mycobacterium avium* complex [MAC] (**Figure 1**). This complex shares an estimated 40% nucleotide sequence similarity with members of the MTBC based on proteome-derived DNA sequences (17). The MAC comprises several closely related, slow-growing, pathogenic, and nonpathogenic species. Among the pathogenic species within the MAC are *M. avium* and its subspecies (MAP; *M. avium* subsp. *avium*; *M. avium* subsp. *silvaticum*; and *M. avium* subsp. *hominissuis*) and *Mycobacterium intracellulare*. Phylogenetic analyses based on partial gene, whole gene, and complete genome DNA sequences have revealed that subspecies of *M. avium* typically share over ≥95% nucleotide sequence identity, while nucleotide identity between *M. avium* subspecies and *M. intracellulare* was estimated between 80 and 94% (18–21).

Despite their genetic similarity, members of the MAC show differential host and tissue tropisms. For example, *M. avium* subsp. *avium* is the classical causative agent of tuberculosis in birds, while *M. avium* subsp. *silvaticum* has been shown to cause tuberculous lesions in wood pigeons. In this regard, *M. avium* subsp. *avium* and *M. avium* subsp. *silvaticum* together represent a distinct lineage of avian pathogens. *M. avium* subsp. *hominissuis* and *M. intracellulare* are opportunistic pathogens widely distributed in the environment and can cause disseminated tuberculosis and pulmonary disease in a range of mammalian hosts, including pigs, cattle, and human beings (22). From a human perspective, *M. avium* subsp. *hominissuis* is considered to be most clinically relevant member of the MAC member, where it has been previously shown to have caused disseminated infections among immunocompromised patients (23). MAP, in contrast, is an obligate intracellular pathogen of ruminants that causes JD characterized by chronic enteritis, with severe economic losses for the dairy industry in many countries (3). MAP can be differentiated from the other subspecies of *M. avium* by its very slow growth rate *in vitro* (between 8 and 24 weeks growth is required for visible colony formation). In addition, MAP is dependent on the siderophore

mycobactin J, an iron-chelating cell wall component, for growth in primary cultures (24).

#### **Mycobacterium bovis INFECTION AND BOVINE TUBERCULOSIS**

Bovine tuberculosis is caused by infection with *M. bovis*, which continues to pose a threat to livestock worldwide. Furthermore, as a zoonotic pathogen, *M. bovis* also has serious implications for human health (25). It has been estimated that BTB contributes losses of \$3 billion to global agriculture annually (26, 27), while comprehensive econometric analyses place BTB as the fourth most important livestock disease worldwide (28). The impacts of BTB infection are manifold, including significant economic and social effects due to the slaughter of infected animals, compensatory payments to producers, continual surveillance programs, and disruption to agricultural trade and productivity (29).

Bovine tuberculosis is predominantly a pulmonary disease, characterized largely by the formation of tuberculous lesions in the upper respiratory lymph nodes of the lung and thorax. In some cases, tuberculous lesions have also been detected in the cranial lymph nodes (30). The etiology and host immune response to *M. bovis* is similar to *M. tuberculosis* infections in human beings (31). Infection is normally caused by the inhalation of aerosolized respiratory secretions containing infectious bacilli, with the natural site of infection being the respiratory tract, presumably on the alveolar surface of the lung. Following inhalation, the bacilli are rapidly encountered by host alveolar macrophages and other phagocytic innate immune cells (such as dendritic cells), which serve as key innate immune effector cells that provide the first line of defense against the pathogen. At this stage, bacilli can be destroyed by the antimicrobial actions of the macrophage; however, bacilli that evade intracellular destruction can persist and multiply within infected macrophages. This results in the migration of infected macrophages to regional lymph nodes, where protective TH1 cell immunity is induced through the recruitment and interaction of additional innate and adaptive immune cells, culminating in the formation of granulomas – organized complexes of immune cells composed of lymphocytes, non-infected macrophages, and neutrophils that contain infected macrophages and prevent the dissemination of bacilli (5, 31, 32). However, in most cases, the pathogen is not eliminated by the host; rather, the pathogen persists in a dormant stage within the granuloma for prolonged periods of time, becoming metabolically and reproductively active following the breakdown of the granuloma and dysregulation of protective T-cell immunity. This results in the development of active tuberculosis, causing immunopathology in the host and enabling the transmission of infection (31).

#### **MAP INFECTION AND JOHNE'S DISEASE**

Johne's disease, caused by infection with MAP, is a chronic inflammatory disease that affects the gastrointestinal tract of cattle and other ruminants. Specifically, JD presents as a granulomatous inflammation of the intestinal tissue and regional lymph nodes due to a massive influx of monocytes and macrophages. This inflammation effectively prevents absorption of nutrients, and, therefore, during the later stages of disease cattle manifest significant weight loss and diarrhea, resulting in progressive physiological wasting

and death. Disease progression is generally classified into four stages: silent infection, subclinical, clinical, and advanced clinical disease; in particular, the subclinical phase can be extremely lengthy (between 2 and 5 years), with pathology largely restricted to the ileum, rendering diagnosis difficult (24).

Exposure to MAP in ruminants generally occurs within the first months of life, through either a fecal–oral route or by ingestion of contaminated colostrum or milk, although evidence suggests that some infections can occur *in utero* (33). Following ingestion, the bacilli colonize the mucosa-associated lymphoid tissues of the upper gastrointestinal tract and are subsequently endocytosed by the microfold cells (M-cells) that cover the ileal Peyer's patches. The bacilli are subsequently phagocytosed by subepithelial and intraepithelial intestinal macrophages, where they reside and multiply (34). The subsequent host cellular immune response leads to the development of granulomas involving adjacent lymph nodes. After years of latent infection, bacilli are assumed to reactivate and trigger a state of active cellular proliferation, leading to the corrugated intestinal epithelia and clinical manifestations with shedding of bacilli into the environment and grassland completing the infection cycle. The disease finally presents as a malnutrition syndrome that culminates in the death of the animal (3, 35, 36).

Johne's disease has major implications for domestic animal health worldwide causing significant economic loss in affected herds, which is largely due to decreased milk yields, reduced slaughter weight, premature culling of infected animals, and losses due to continued spread of infection (37). In cattle, JD results in estimated losses of \$250 million to the US dairy industry annually, while dairy herd prevalence of JD is estimated to be greater than 50% in certain US states and European countries (35, 38, 39). Furthermore, it has been hypothesized that MAP infection may trigger or exacerbate Crohn's disease, an inflammatory disease of the intestines in human beings with similar granulomatous pathology at the ileocecal valve; however, this proposed link between MAP infection and Crohn's disease remains contentious (40).

#### **THE ROLE OF THE MACROPHAGE DURING MYCOBACTERIAL INFECTION**

Although BTB and JD exhibit distinct complex etiologies, the causative agents of these diseases display a propensity to infect, reside, and replicate in host macrophages – the key host innate immune cell that mediates the immune response following infection. Macrophage recognition of mycobacteria occurs through the interaction of mycobacterial pathogen-associated molecular patterns (PAMPs) – such as lipopolysaccharide, and various lipoproteins and glycolipids (e.g., lipoarabinomannan) – with host pathogen recognition receptors (PRRs) proteins, such as Toll-like receptors (TLRs), which are expressed on the macrophage cell surface (41). Macrophage PRR activation induces signaling pathways resulting in the production of endogenous NF-κB-inducible cytokines and chemokines that promote a TH1 immune response characterized by the release of proinflammatory IFN-γ, primarily from CD4<sup>+</sup> T-cells, and the lysing of infected macrophages by cytotoxic CD8<sup>+</sup> T-cells. IFN-γ induces microbicidal activity in infected macrophages and enhances the expression of major histocompatibility complex (MHC) class I and II molecules necessary for the presentation of mycobacterial antigens on the

macrophage cell surface to CD8<sup>+</sup> and CD4<sup>+</sup> T-cells, respectively. These mechanisms can lead to either the immediate killing of the pathogen and clearing of infection, or the containment of infection through the formation of granulomas (42–45).

Pathogenic mycobacteria have evolved a diverse range of immunoevasive mechanisms that facilitate survival and replication within the host macrophage. These immunoevasive mechanisms include inhibition of phagosome maturation necessary for destruction of the pathogen and antigen presentation (46, 47); evasion of macrophage apoptosis and activation of macrophage necrosis, which facilitates release of bacilli from the macrophage and encourages dissemination of infection to other cells (7, 48); and the subversion of innate cell signaling, which is critical to the establishment of infection and progression to active disease (49, 50). It has also been recently demonstrated that virulent *M. tuberculosis* strains preferentially infect permissive macrophages and evade microbicidal macrophages through the masking of PAMPs with cell surface associated lipids (51).

Failure or subversion of an appropriate innate immune response is critical to the establishment of infection and progression to disease; central to this process is the macrophage response to infection (31). Consequently, analysis of the bovine macrophage response to *in vitro* infections with *M. bovis* and MAP may provide insights into the cellular mechanisms that underlie and govern the divergent immunopathology of BTB and JD (36, 52).

#### **FUNCTIONAL GENOMICS ANALYSIS OF THE BOVINE MACROPHAGE RESPONSE TO MYCOBACTERIA**

Early investigations of the bovine macrophage response to mycobacterial infection focused on the analysis of the expression of single or small numbers of immunological parameters. For example, the quantification of gene or protein expression using reverse transcription quantitative real-time PCR (RT-qPCR) and ELISA technologies; however, focused studies such as these, are unable to provide a high-resolution overview of the global macrophage response to infection. Pathogen-induced activation of host macrophages is characterized by large-scale changes in the expression profile of genes critical for the control and eradication of the pathogen, while modulation of host gene expression critical for pathogen survival is also expected to be reflected in the transcriptome of the macrophage (53, 54). Consequently, genomics technologies that assay pan-genomic changes in gene expression have been widely used to discern patterns of host-gene regulation during infection. In particular, the development of highthroughput gene expression technologies, such as microarrays and RNA-sequencing (RNA-seq), over the past decade, coupled with dramatic improvements in mammalian genome resources and increasingly sophisticated computational tools for the analysis of large-scale gene expression datasets are providing new opportunities for detection, cataloging, and analysis of the large numbers of host macrophage genes expressed in response to mycobacterial infection in cattle (55–64).

A primary goal of our research group is to use highthroughput functional genomics technologies to analyze the bovine macrophage transcriptome following infection with *M. bovis* and *M. paratuberculosis* to improve our understanding of host–pathogen interactions that characterize and underlie BTB and JD. Previously, we used a 24-h time course infection model to investigate the transcriptome response of bovine monocytederived macrophages (MDM) infected with *M. bovis* and MAP using data generated from the Affymetrix® GeneChip® Bovine Genome microarray platform (61, 62). These analyses revealed a large number of differentially expressed (DE) genes following *M. bovis* infection relative to non-infected control MDM such that the number of DE genes also increased across the time course, with the highest number observed 24 h post-infection [hpi] (62). This contrasted with results from MAP-infected MDM relative to the same control macrophages, which showed the highest number of DE genes at the 2 hpi time point with a decrease in the number of DE genes at the later time points post-infection (i.e., 6 and 24 hpi) (61). These findings suggest that *M. bovis* and MAP have differential survival strategies once internalized by macrophages, which in turn, may underlie the divergent immunopathology associated with BTB and JD.

#### **THE MONOCYTE-DERIVED MACROPHAGE INFECTION MODEL AND GENE EXPRESSION DATASETS USED FOR COMPARATIVE FUNCTIONAL ANALYSIS**

To further investigate the similarities and differences of the bovine MDM response to virulent and attenuated mycobacterial species/strains, we have reanalyzed and directly compared the Affymetrix® GeneChip® Bovine Genome microarray data from our earlier work (i.e., the non-infected control MDM and the *M. bovis*- and MAP-infected bovine MDM) together with corresponding microarray data from MDM infected with *M. bovis* BCG (65). All infected and control MDM used to generate these data were derived from the same seven age-matched Holstein-Friesian females, while a multiplicity of infection (MOI) of 2:1 (i.e., 2 bacilli:1 MDM) was used for all MDM infections (61, 62). Gene expression omnibus (GEO) data series accession numbers used for these re-analyses were GSE33309, GSE35185, and GSE59774 (66).

For this new comparative analysis, gene expression data from *M. bovis*-, MAP-, and *M. bovis* BCG-infected MDM together with data from the non-infected control MDM at time points 2, 6, and 24 hpi were used. Prior to differential gene expression analysis, all microarray data were quality checked using the arrayQualityMetrics package in Bioconductor (67). Raw data from two microarrays (one MAP- and one *M. bovis* BCG infected MDM sample) did not pass the QC thresholds set in the arrayQualityMetrics package and were removed from all further downstream analyses. Furthermore, all arrays generated from these two animals were excluded to ensure a balanced experimental design for comparative gene expression analysis (**Figure 2**). Next, all raw gene expression data

**functional genomics analysis described in the current study**. Previously, MDM from seven age-matched females were infected with M. bovis and MAP (MOI 2:1); control MDM received culture media only (61, 62). We also infected, in parallel, MDM from the same animals with M. bovis BCG (MOI 2:1). Infections were performed across duplicate tissue culture plate wells (shaded circles); total RNA from duplicate

Pan-genomic gene expression data for each RNA sample was generated using the Affymetrix® GeneChip® Bovine Genome microarray platform (61, 62). For the comparative functional genomics analysis, microarray data from only five of these animals were used to ensure a balanced experimental design following quality control assessment (see main body text for details).

were normalized and filtered for non-informative probe sets using the I/NI algorithm implemented in the FARMS Bioconductor software package (version 1.14.0) (68); this analysis yielded 11,842 informative probe sets for use in downstream analysis.

The 11,842 informative probe sets identified post-filtering were then used to generate a multi-dimensional scaling (MDS) plot summarizing the transcriptomic relationship between samples (**Figure 3**). Notably, samples clustered according to time on the first dimension, while the second dimension clustered samples according to animal ID. However, there was a noticeable clustering of RNA samples from *M. bovis* BCG- and MAP-infected samples within each group of samples corresponding to a particular animal at a given time. Investigation of the expressed genes with the greatest animal effect revealed that loci within the bovine MHC displayed the greatest differences in expression among animals with relatively small gene expression differences within each animal across time and treatment (Figure S1 in Supplementary Material). The magnitude of the expression differences among animals for these genes presumably contributes to the separation of the samples by animal on the second dimension in the MDS plot (**Figure 3**).

We propose two hypotheses to explain the observed expression differences among animals. First, expressed MHC loci, which are among the most polymorphic loci in mammals, generate mRNA transcripts that display considerable nucleotide differences between individuals from the same species (69–71). mRNA transcripts that vary appreciably from the reference transcript sequences used to produce the microarray probe sets may, therefore, have a reduced hybridization efficiency compared to mRNA transcripts that are identical to or differ only slightly from the reference transcript sequences. In turn, this may result in Type 1 errors for differential gene expression estimates between samples or sample groups (72–74). Consequently, the animal effect observed in these data may be due, in part, to the technical limitations of the microarray platform used. Second, it is possible that the observed inter-animal differential gene expression is due, in part, to real differences in mRNA abundance generated by genotypic differences at loci that regulate gene expression (71, 74, 75); indeed, genotypic differences at loci that regulate gene expression in response to mycobacterial infection may contribute to phenotypic differences in the ability of an animal to clear or succumb to infection. In addition, a combination of both these technical and biological factors could also explain the observed animal effect.

#### **COMPARATIVE FUNCTIONAL GENOMICS ANALYSIS REVEALS SIMILARITIES AND DIFFERENCES IN THE MACROPHAGE TRANSCRIPTOME RESPONSE TO M. bovis, MAP, AND M. bovis BCG**

**Figure 4** shows the results of differential gene expression analysis for the infected MDM (i.e.,*M. bovis*, MAP, and *M. bovis* BCG) relative to the non-infected control MDM [false-discovery rate (FDR) adjusted *P* value ≤0.05]. Notably, for each infected MDM/control, MDM contrast the number of DE genes varied with respect to time. *M. bovis*-infected MDM exhibited the greatest number of DE genes, with the number of DE genes increasing across the 24-h time course. Furthermore, for *M. bovis*, the number of downregulated genes exceeded the number of upregulated genes at each of the time points post-infection.

In contrast, the number of DE genes observed for both the MAP- and *M. bovis* BCG-infected MDM (relative to the

non-infected controls) was highest at the 6 hpi time point for both sample groups (**Figure 4**). These results indicate that for MAP and *M. bovis* BCG infection, MDM differential gene expression had largely abated at the 24 hpi time point and that the MDM transcriptome reverted to a transcriptional state similar to that of the control MDM. These observations support previous work that showed that differential gene expression changes in MAPinfected bovine MDM are transient and are largely undetected 24 hpi relative to non-infected control MDM (76, 77).

A comparison of the lists of DE genes obtained for all mycobacteria/control contrasts at each time point (**Figure 5**) identified a core set of DE genes at the 2 and 6 hpi time points common to all three mycobacterial treatments consisting of 170 and 236 DE genes, respectively. Among the DE genes common to all three mycobacterial infections were *IL1A*, *IL1B*, *TNF*, *NFKB1*, and *NFKB2*; all of these genes were upregulated at one or more time points in all types of infected MDM, suggesting a robust inflammatory reaction to all of the mycobacteria used in this study. Ingenuity® Systems Pathway Analysis (IPA; Ingenuity Systems, Redwood City, CA, USA; www.ingenuity.com) was used to identify canonical pathways within the list of DE genes that were common to all three mycobacterial treatments relative to the control group. The top ranking canonical cellular pathways (based on the lowest adjusted *P*-values; FDR ≤ 0.05) enriched for the 170 and 236 common DE genes identified at 2 and 6 hpi included *IL-10 signaling* (2 hpi); *dendritic cell maturation* (2 and 6 hpi); *IL-6 signaling* (2 hpi); *TNFR2 signaling* (2 hpi), *TWEAK signaling*

(2 hpi); *communication between innate and adaptive immune cells* (6 hpi), *NF-*κ*B signaling* (6 hpi); and *TREM1 signaling* (2 and 6 hpi; FDR ≤ 0.0001) (see Tables S1 and S2 in Supplementary Material).

We also observed large numbers of DE genes that were specific to *M. bovis*-infected MDM, which increased over time from 1,013 significant genes at 2 hpi to 2,167 at 24 hpi (**Figure 5**). The percentages of DE genes unique to *M. bovis* infection relative to the total number of DE genes detected following *M. bovis* infection were 80.3% (2 hpi), 78.8% (6 hpi), and 96.6% (24 hpi). IPA analysis of the DE genes unique to *M. bovis*-infected macrophages across all three time points (Tables S3–S5 in Supplementary Material) revealed enrichment for genes involved in cell signaling, including *IL-6 signaling* and *mitogen-activated protein kinase (MAPK) signaling*, the latter of which regulates the expression of several transcription factors, such as those encoded by *FOS* and *JUN* that are critical for the activation of immune cells (78). In contrast, the number of DE genes specific to MAP and *M. bovis* BCG infection across the time course was markedly lower; for example, no gene was specific to either MAP or *M. bovis* BCG infection 24 hpi. Indeed, for each post-infection time point, more than 98.3% of the DE genes induced by MAP and *M. bovis* BCG relative to the controls were among the list of DE genes induced by *M. bovis* relative to the controls.

We next analyzed differential gene expression directly between each pair of mycobacteria-infected sample groups (**Figure 6**). Large and increasing numbers of DE genes were found across the time course in *M. bovis*-infected MDM relative both MAP- and BCG-infected MDM, confirming the divergence between the two types of MDM transcriptional responses. Notably, no DE genes were detected between MAP- and BCG-infected MDM at 2 and 6 hpi, while only two DE genes were identified at 24 hpi; this contrasts with the larger numbers of DE genes identified through the indirect comparison of MAP- and BCG-infected MDM using the controls as a common reference (**Figure 5**). This discrepancy is most readily explained by differences in variances of gene expression within each sample group as illustrated in Figure S2 in Supplementary Material.

The overlap in DE genes across all three treatment groups (relative to the control groups) at the 2 and 6 hpi time points suggest that a "core" MDM transcriptional response is induced by all three mycobacterial species/strains during the early stages of infection. This core transcriptome is characterized by genes involved in innate cytokine signaling and production, which encode proteins that activate the adaptive immune response following mycobacterial infection. However, the large and increasing number of DE genes specific to *M. bovis* across the infection time course demonstrates that *M. bovis* is a more potent inducer of proinflammatory genes than *M. bovis* BCG or MAP and highlights the distinct MDM gene expression profile elicited in response to this pathogen. In addition, the enrichment of *M. bovis*-specific DE genes for roles in macrophage cell signaling, such as IL-6 and MAPK signaling, suggests that additional cellular pathways are triggered by this pathogen relative to *M. bovis* BCG or MAP.

Although proinflammatory cytokines and chemokines play a pivotal role in mediating the host immune response to control mycobacterial infection, several lines of evidence suggest

that these molecules and their associated pathways can be exploited by virulent mycobacterial pathogens to promote granuloma formation, which recruit new macrophages to the site of infection enabling persistence within the host (79). For example, non-regulated production of proinflammatory cytokines and chemokines can result in immunopathology, including destructive inflammation and necrosis, allowing dissemination of the pathogen from infected cells (9, 80, 81). Therefore, the immunopathology of BTB may be associated with the increased induction of innate immune genes following infection with *M. bovis*. Furthermore, the divergent transcriptomic profile observed in MDM infected with virulent *M. bovis* relative to *M. bovis* BCGinfected MDM is presumably governed, in part, by the presence of several secreted virulence factors, such as those encoded by the RD1 locus (57, 82, 83). This locus is present in all virulent strains of *M. bovis* (and *M. tuberculosis*) but is absent in attenuated strains of *M. bovis* BCG (84).

Conversely, the reduced number of DE genes detected between the MAP-infected and the non-infected control and *M. bovis* BCGinfected MDM across the time course suggests that MAP infection does not result in a major perturbation of the MDM transcriptome; rather, bovine MDM sense and respond to MAP in a similar manner to attenuated *M. bovis* BCG despite their markedly distinct evolutionary histories and different pathogenicities. These results support the hypothesis that MAP infection of host MDM is achieved via a capacity to appear "benign," which enables it to reside and replicate within the macrophage for prolonged periods of time, and may underlie the lengthy subclinical phase of infection characteristic of JD (85). Our results also suggest that the immunoevasive mechanisms used by MAP involve the suppression of the proinflammatory response, such that the transcriptome of an infected macrophage resembles that of a non-infected cell. In support of this, *in vivo* MAP infection models demonstrate that cattle initially develop an early proinflammatory and TH1-type response to infection, which gradually declines in animals that progress to active disease, favoring a TH2-type response that does not control infection (3, 86–88). Notably, the immunosuppression observed *in vivo* may originate at a cellular level, whereby MAP-infected macrophages fail to properly respond to host-derived immune activators such as CD40L and IFN-γ (89, 90).

#### **TRANSCRIPTIONAL EVIDENCE FOR TYPE I INTERFERON-MEDIATED REGULATION OF INTERLEUKIN I PRODUCTION IN M. bovis-INFECTED MDM: TOWARD A MECHANISM OF PATHOLOGY**

Examination of the lists of DE genes for each mycobacterial infection/control contrast shows that type I interferon-inducible genes such as *IFIT1*, *IFIT2*, *MX1*, *MX2*, and *IL27* and interferondependent *CXCL10* were not DE following MAP and *M. bovis* BCG infection at any of the post-infection time points. However, all of these genes were DE for at least one time point following *M. bovis* infection (**Table 1**). Similarly, the gene encoding type II interferon (i.e., *IFNG*) was also upregulated at 6 and 24 hpi following *M. bovis* infection, but was not DE at any post-infection time point following infection with MAP and *M. bovis* BCG. We further observed that, in general, the fold-change of upregulation of the type I and type II interferon-inducible genes increased over

the time course of infection in the *M. bovis*-infected MDM, with an accompanying decrease in the fold-change of upregulation of the interleukin-1 genes (*IL1A* and *IL1B*).

IL-1 and type I IFN-signaling pathways have been shown to play important, yet opposing, roles in determining the host response to infection with virulent members of the MTBC. Mice deficient in IL-1B display increased susceptibility to virulent *M. tuberculosis*, indicating that IL-1 signaling is required for the host control of infection (91, 92). Conversely, mice deficient in type I IFN signaling show reduced bacterial loads following infection, suggesting that type I IFN plays a contributory role in tuberculosis disease progression (93, 94). Studies have also shown that virulent *M. tuberculosis* and attenuated *M. bovis* BCG use distinct signaling pathways for regulating IL-1B production in human MDM. *M. tuberculosis* induced the expression of IFN-related genes, while induction of type I IFN-signaling inhibited IL-1B secretion (95, 96). Notably, infection of human MDM with *M. bovis* BCG did not induce significant differential expression of IFN-related genes or IL-1B secretion. These results suggest that type I IFN-mediated suppression of IL-1B production is a key mechanism for the intracellular survival of *M. tuberculosis* (95, 96). Comparable *in vivo* and *in vitro* studies of bovine ileal tissue and MDM have demonstrated increased transcript and protein levels of IL-1A and IL-1B (*in vitro* only) in response to MAP infections, with a concomitant increase in downstream expression of TRAF1 (97, 98). This increase in TRAF1 has been proposed to enhance survival of MAP

**Table 1 | Log<sup>2</sup> fold-change values of manually curated genes for several immune-related pathways**.


Cells containing values are statistically different for the contrasts listed in the column headings (FDR ≤ 0.05); empty cells are not statistically different (FDR ≥ 0.05). The intensity of the shading is related to relative increased fold-changes in gene expression.

in macrophages due to the anti-apoptotic properties of TRAF1 and its capacity to interfere with normal macrophage activation, particularly via CD40/CD40L (98, 99).

The increase in type I IFN gene expression and the concomitant decrease in *IL1* induction that we observed in *M. bovis*-infected MDM over the 24-h infection time course suggests that the intracellular survival strategy of virulent *M. bovis* may also involve type I IFN-mediated suppression of *IL1* production. In contrast to this, these results suggest that the same immunoevasive mechanism is not used by virulent MAP or the attenuated *M. bovis* BCG in bovine MDM. These results indicate that differential activation of macrophage immunoregulatory pathways is central to the differential intracellular survival mechanisms of these related, yet distinct, bovine mycobacterial pathogens. A proposed model of the differential response of bovine MDM to the mycobacteria examined in this comparative study is shown in **Figure 7**.

#### **DIFFERENTIAL MYCOBACTERIAL VIRULENCE FACTORS AND THEIR IMPACT ON MACROPHAGE PATHOGEN RECOGNITION**

The region of difference 1 (RD1) locus,which is present in*M. bovis*, is a major genetic difference between this species and MAP and *M. bovis* BCG, which both lack RD1 (100). Moreover, RD1 (which is also present in virulent *M. tuberculosis* strains) has attracted increasing interest over the last two decades, because it contains the ESX-1 type VII secretion system, responsible for the secretion of virulence factors, such as the dimer ESAT-6/CFP-10 (101). ESAT-6, which is proposed to facilitate escape from macrophage phagolysosomal degradation, binds to TLR2 receptors, and activates TLR signaling cascades within the macrophage that culminate in cytokine production (102–104). It has been recently reported that TLR2 receptors mediate enhanced interferon production through

reprograming of murine macrophages following infection with viral ligands (105). In support of these observations, our results showed that differential expression of both type I and type II interferon genes is unique to *M. bovis*-infected MDM, which is not present in MDM infected with MAP and *M. bovis* BCG. Consequently, we hypothesize that the virulence factors encoded in the RD1 region and secreted by *M. bovis* – but not MAP or *M. bovis* BCG – trigger an additional cascade of signaling events, such as those mediated by TLRs and MAPKs (as revealed through IPA analyses of unique *M. bovis*-induced genes), relative to attenuated *M. bovis* BCG and the lengthy subclinical MAP. In turn, the combined activation of additional immune pathways and canonical PRR-dependent pathways may lead to a sustained (i.e., chronic) inflammatory response in infected macrophages, as opposed to a more transient inflammation following *M. bovis* BCG or MAP infections.

#### **CONCLUDING REMARKS**

In the present study, we highlight markedly different transcriptional response of bovine MDM infected with *M. bovis* over a 24-h time course compared to the closely related but attenuated *M. bovis* BCG strain and to virulent MAP. We hypothesize that RD1-encoded virulence factors provide a mechanistic basis for this differential response, as RD1 was lost during the derivation of the *M. bovis* BCG vaccine strain from *M. bovis* and is absent from the MAP genome. We also identified a common MDM transcriptional response to both attenuated *M. bovis* BCG and MAP. We propose that the respective attenuated and lengthy subclinical phenotypes of *M. bovis* BCG and MAP may induce similar responses in infected macrophages, at least during the early stages of infection. Finally, we identified type I interferon-dependent

genes among the DE genes specific to virulent *M. bovis-*infected MDM, adding further evidence supporting a key role for type I interferon in the establishment of active tuberculosis in cattle as has previously reported for human tuberculosis (106).

While the comparative functional genomics analysis presented here is based on data generated from microarrays, the changing landscape of transcriptomics, as represented by the advent of highthroughput RNA-seq, offers unprecedented opportunities to study the host macrophage response to mycobacterial infection at the nucleotide level, including investigation of the effect of genotype on gene expression levels. High-throughput sequencing technologies are providing novel insights into the cellular mechanisms governing mycobacteria–macrophage interactions, enabling further understanding of how modulation of these pathways can result in pathology. In addition, identification of transcriptional biomarkers of infection may lead to the development of novel diagnostics for BTB and JD, providing new molecular tools for disease control and eradication.

#### **ACKNOWLEDGMENTS**

This work was supported by Investigator Grants from Science Foundation Ireland (Nos: SFI/01/F.1/B028 and SFI/08/IN.1/B2038), a Research Stimulus Grant from the Department of Agriculture, Food and the Marine (No: RSF 06 405), and a European Union Framework 7 Project Grant (No: KBBE-211602- MACROSYS). Kévin Rue-Albrecht is supported by the UCD Wellcome Trust Computational Infection Biology PhD program (No: 097429/Z/11/Z). Kate E. Killick was supported by the Irish Research Council (IRC) funded Bioinformatics and Systems Biology thematic PhD program (http://bioinfo-casl.ucd.ie/PhD). We would like to thank Mr. Eamon Costello and staff at the Irish Department of Agriculture, Food and the Marine Backweston Laboratory Campus (Celbridge, Co., Kildare) for providing the isolate of *M. bovis*. We also would like to thank the staff at the UCD Lyons Research Farm for assistance with blood collection and Dr. Eamonn Gormley (UCD) and Prof. Joseph Keane (TCD) for valuable comments.

#### **SUPPLEMENTARY MATERIAL**

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

#### **REFERENCES**


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


106. Berry MP, Graham CM, McNab FW, Xu Z, Bloch SA, Oni T, et al. An interferoninducible neutrophil-driven blood transcriptional signature in human tuberculosis. *Nature* (2010) **466**:973–7. doi:10.1038/nature09247

**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.The Guest Associate Editor Kieran G. Meade declares that, despite having collaborated with authors Kévin Rue-Albrecht, David A. Magee, Kate E. Killick, Nicolas C. Nalpas, Stephen V. Gordon and David E. MacHugh, the review process was handled objectively and no conflict of interest exists.

*Received: 21 July 2014; paper pending published: 01 September 2014; accepted: 10 October 2014; published online: 05 November 2014.*

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

*This article was submitted to Molecular Innate Immunity, a section of the journal Frontiers in Immunology.*

*Copyright © 2014 Rue-Albrecht , Magee, Killick, Nalpas, Gordon and MacHugh. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

## The single intradermal cervical comparative test interferes with Johne's disease ELISA diagnostics

#### **Aideen E. Kennedy 1,2, Ana T. Da Silva<sup>1</sup> , Noel Byrne<sup>1</sup> , Rodney Govender <sup>2</sup> , John MacSharry <sup>3</sup> , Jim O'Mahony <sup>2</sup> and Riona G. Sayers <sup>1</sup>\***

<sup>1</sup> Animal and Bioscience Research Department, Animal and Grassland Research and Innovation Centre, Teagasc, Cork, Ireland

<sup>2</sup> Department of Biological Sciences, Cork Institute of Technology, Cork, Ireland

<sup>3</sup> Alimentary Pharmabiotic Centre, Biosciences Institute, University College Cork, Cork, Ireland

#### **Edited by:**

Uday Kishore, Brunel University, UK

#### **Reviewed by:**

Srinand Sreevastan, University of Minnesota, USA Divyendu Singh, Indiana University Bloomington, USA

#### **\*Correspondence:**

Riona G. Sayers, Animal and Bioscience Research Department, Animal and Grassland Research and Innovation Centre, Teagasc, Moorepark, Fermoy, Co. Cork, Ireland e-mail: riona.sayers@teagasc.ie

Enzyme-linked immunosorbent assays (ELISA) of milk and serum samples are a routinely used method of screening herds for Mycobacterium avium subspecies paratuberculosis (MAP). Infection with MAP causes granulomatous enteritis of ruminants known as Johne's disease (JD). The sensitivity (Se) and specificity (Sp) of MAP ELISAs leads to difficulties in the identification of both infected and infectious animals. Interference with MAP ELISA Se and Sp has been reported in MAP seronegative cows following administration of purified protein derivative (PPD) as part of intradermal testing for bovine tuberculosis (bTB). The aim of this study is to examine the impact of the single intradermal cervical comparative test (SICCT) for bTB, on both serum and milk MAP ELISA tests, in a herd containing both seropositive and seronegative cows pre-SICCT. A secondary objective is to provide appropriate timing of JD ELISA tests in relation to the SICCT. A herd of 139 cows were serum and milk sampled pre- and post-SICCT administration. Prior to SICCT, 6% of the herd tested seropositive for MAP using milk ELISA, with 8% positive on serum. ID Screen Paratuberculosis Indirect Screening Test (ID Vet) was used to screen the herd. Within 14 days of PPD administration, a significant increase in the prevalence of seropositive cows was recorded. Identical prevalence's were recorded with both test matrices (39%). ELISA values remained significantly higher until day 43 post-SICCT in milk (P = 0.850), and day 71 in serum (P = 0.602). If the "new" positives detected post-bTB testing are deemed false positives due to generation of cross-reacting antibodies by administration of PPD, milk would appear a more suitable sample for JD ELISA testing within 2 months of SICCT. In summary, sampling for JD utilizing milk ELISA should be avoided in the 43-day period following PPD administration, with serum ELISA sampling avoided for an additional 28 days.

**Keywords: Mycobacteriacea, Johne's disease,TB test, ELISA, PPD**

#### **INTRODUCTION**

*Mycobacterium avium* subspecies *paratuberculosis* (MAP), a member of the Mycobacteriacea family, causes chronic granulomatous enteritis known as Johne's disease (JD) (1). Clinical JD is characterized by diarrhea and progressive cachexia, which ultimately results in death (2). Uncertainty exists regarding a potential causal link between MAP and Crohn's disease in humans (3, 4). The potential damage to the global dairy industry, should a link between Crohn's and MAP be fully substantiated (5), combined with impacts on animal health, has prompted the establishment of JD control programs in a number of countries (6–8).

Use of enzyme-linked immunosorbent assays (ELISA) to identify animals at risk of being infected with MAP is common in control programs internationally (8, 9), including Ireland (10). ELISA is favored as a screening test due to its relatively low cost compared to fecal culture or polymerase chain reaction (PCR) (11). ELISAs also provide timely results compared to culture methods (11). The sensitivity (Se) and specificity (Sp) of MAP ELISAs, however, leads to difficulties in the identification of both infected and infectious individuals (12).

*Mycobacterium bovis*, the causative agent of bovine tuberculosis (bTB), is an additional pathogenic and definitively zoonotic (13) member of the Mycobacteriaceae. To reduce the zoonotic risk posed by bTB, address public/animal health concerns, and limit trade restrictions, a compulsory national eradication program for bTB was established in Ireland in 1962 (14). This eradication program involves ante-mortem testing of all registered bovines annually using the single intradermal cervical comparative test (SICCT) and post-mortem carcass inspection. All SICCT positive animals (reactors) are slaughtered, the herd of origin is restricted, and additional bTB testing is applied to the herd. The comprehensive nature of the testing program can lead to some animals being tested up to five times in a single year (15).

The SICCT utilizes intradermal introduction of *M. bovis* and *M. avium* subsp. *avium* purified protein derivatives (bPPD and aPPD) at two different sites on the neck to elicit a delayed hypersensitivity response mediated by T cells (16). Comparative measurements at both injection sites, taken 72 h post-PPD administration, are used to assess infection status (16). Additional ante-mortem testing methods used internationally for detection of bTB include

the single intradermal test and the caudal fold test,both less specific than SICCT (17).

Members of the Mycobacteriaceae family share several antigens, which can lead to diagnostic difficulties due to antibody cross reaction (18). MAP infection can interfere with specificity of bTB diagnostics (19), and likewise *M. bovis* infection can affect MAP serological tests (20). Varges et al. (21), has also shown interference by both single and comparative intradermal bTB tests on MAP sero diagnostics in bTB negative animals. The primary purpose of this current study was to investigate the impact of SICCT on the prevalence of ELISA positive results (serum and milk) in an Irish herd containing both MAP ELISA seropositive and seronegative animals over a period of 6 months. Secondary objectives included comparing milk and serum ELISA readings and investigating whether serum samples could be taken at the 72 h bTB visit without interference from PPD administration.

### **MATERIALS AND METHODS**

#### **STUDY HERD**

A 139-cow spring-calving dairy herd (mean-calving date February 19th) was recruited. This herd was depopulated in 1997 following a confirmed case of bovine spongiform encephalopathy (BSE). The experimental herd, therefore, consisted of descendants of cows used to repopulate the farm in 1998 (22). Annual statutory bTB test results were sourced from 1998 to provide a bTB history for the herd. Veterinary records were obtained in order to record a JD history for the herd post-repopulation. Approximately 60% of the cows were Holstein Friesian (HF), the remaining 40% purebred Jersey (Je) or Je cross-breeds. The study was licensed by the Irish Department of Health and Children.

#### **SAMPLE COLLECTION**

Milk and serum samples were collected 10 and 13 days prior to administration of the compulsory annual SICCT herd test in May 2012 (pre-SICCT). The SICCT was administered by the Department of Agriculture, Food and the Marine (DAFM) approved private veterinary practitioner (PVP) responsible for the care of animals on this farm as is standard practice for the Irish national bTB eradication scheme. Milk and serum samples were collected every 14 days (approximately) for 2 months post-SICCT and on a monthly basis thereafter until the composition of the herd changed materially due to end of lactation culling (longitudinal data). Sampling dates for serum and milk samples are outlined in **Table 1**. A limit of a 7-day interval between serum and milk sampling was applied in order to consider samples as "matched." Milk samples were not available for all cows at every sampling time point which is reflected as small variations in sample sizes. Additionally, milk samples were not collected in September 2012 due to an un-related health issue on farm. Fecal samples were collected on a weekly basis from consistently ELISA positive cows from 90 days post-SICCT. These cows were also subjected to a veterinary clinical exam.

Serum and milk samples were tested using a commercial ISO17025 accredited laboratory (designated laboratory for Irish voluntary JD control program) using the ID Screen Paratuberculosis Indirect Screening Test (ID Vet, Montpellier, France). The test is an *M. phlei* absorbed ELISA detecting anti-MAP IgG. Status of the sample was evaluated by examining the **Table 1 |Timetable of serum and milk samples and dates of SICCT**.


sample to positive ratio (*S*/*P* ratio) calculated using the formula *S*/*P* Ratio = [(ODSample − ODPositive control) ÷ (ODPositive control − ODNegative control) × 100]. Fecal samples were tested by microbial culture and real-time PCR (rtRT-PCR) using "in-house" methodologies developed by Cork Institute of Technology as outlined by Douarre et al. (23). The target gene was IS900. Primer sequences for the amplification were 5'-GAAGGGTGTTCGGGGCCGTCG CTTAGG-3' and 5'-GGCGTTGAGGTCGATCGCCC ACGTGAC-3' (reverse primer).

#### **DATA ANALYSIS**

Descriptive analysis, dataset construction, and graphical representations were completed in Excel (MS Office 2010). Normality of datasets was examined visually using ladders of power histograms in Stata (version 12)**.** Additional statistical analyses including chisquared test, *t*-test, box plot construction, Spearman rank correlation, and generalized estimating equations (GEE) were completed using Stata (version 12).

For the purposes of reporting within-herd MAP prevalence, ELISA *S*/*P* ratio results were interpreted according to manufacture instructions, i.e., ≥70 *S*/*P* (serum) and ≥15 *S*/*P* (milk) were categorized as positive, with a single exception. Cows recording *S*/*P* ratios of 60 ≥ SP < 70 (normally classified as inconclusive), were also categorized as negative. The prevalence of positive cows within the herd was plotted vs. trial day. Box plots were constructed to highlight trends in ELISA *S*/*P* % readings pre- and post-SICCT

Longitudinal milk and serum ELISA results were used to create datasets for statistical analysis. ELISA results were recorded as both a categorical variable (positive, negative) and a continuous variable (ELISA *S*/*P* readings). Multivariable GEE was used to investigate differences between pre- and post-SICCT categorical and continuous variables (dependent variables). Independent variables included in the models were sampling time point (pre-SICCT, post-SICCT), breed (Friesian, Jersey), parity (parities 1, 2, 3, ≥4), and date of calving (January, February, March, April). Second level

interactions between independent variables were examined and included in the model at *P* ≤ 0.05. For categorical variable analysis, a binomial distribution was assumed and a logit link function used. For continuous variable analysis, a Gaussian distribution and an identity link function was used. An exchangeable correlation was applied to both analyses. To investigate the correlation between milk and serum ELISA results, Spearman correlation (*rs*) was performed on categorical data sets.

#### **RESULTS**

Results of statutory bTB testing for this farm over the past 8 years indicate minimal issues with bTB in this herd. Similarly, no bTB positive reactor was identified following SICCT in 2012. From herd repopulation in 1998 to commencement of this study, no clinical case of JD had been diagnosed on the study farm.

Prior to administration of the SICCT, a total of 11 of 139 cows (7.9%) tested MAP ELISA positive in serum, with 8 of 137 (5.8%) milk samples testing positive. Following administration of SICCT, a significant increase in the prevalence of ELISA positives was recorded on both test matrices (serum *P* < 0.001; milk *P* < 0.001). The highest recorded prevalence of positive results for both serum and milk samples was 39% (**Figure 1**). No statistically significant difference (*P* = 0.668) was recorded in the prevalence of serum positive results, pre- and 72 h post-SICCT. Similarly no statistically significant difference (*P* = 0. 197) was recorded in *S*/*P* ratios of serum ELISA results pre- and 72 h post-SICCT. Both box plots and GEE analysis highlight an increase in both serum and milk *S*/*P* ratio readings subsequent to the 72 h sampling (**Figures 2** and **3**; **Tables 2** and **3**, respectively).

Statistically significant differences between pre- and post-SICCT milk ELISAs were recorded until 43 days postadministration of PPD, examined as both a continuous and categorical variable (**Table 2**). The prevalence of ELISA serum positive samples was not statistically different from pre-SICCT levels by day 58, while serum ELISA *S*/*P* ratios remained significantly elevated for 71 days post-SICCT (**Table 3**). It should be noted that a significant elevation in *S*/*P* ratios post-SICCT was again noted in November (trial day 143) for both milk and serum samples (**Tables 2** and **3**). No significant second level interactions were identified between independent variables.

Spearman correlation analysis of matched serum and milk samples generated pre SICCT values of *r<sup>s</sup>* 0.73. Post SICCT values ranged from *r<sup>s</sup>* 0.55 to 0.79 with the highest levels recorded at post SICCT test 1 (*r<sup>s</sup>* 0.77) and post SICCT test 6 (*r<sup>s</sup>* 0.79).

Weekly fecal culture of consistently ELISA positive cows yielded negative results. A total of 10 animals yielded PCR positive results, 2 of which recorded positive results at each sampling time point. Veterinary examination did not yield any clinical signs of JD in these animals.

#### **DISCUSSION**

The Irish cattle population is subjected to a comprehensive and compulsory bTB eradication program, involving administration of the SICCT on at least an annual basis (15). The purpose of the current study was to investigate the impact of SICCT (i.e., administration of bPPD and aPPD) on both the within-herd prevalence of positive cows and ELISA *S*/*P* ratios in an Irish dairy herd. The results of the current study can provide useful guidance to farmers and veterinarians on the optimum period to conduct MAP ELISA testing in regions engaging in comprehensive testing for bTB using SICCT.

Two international studies, one conducted in Brazil (21), and the second in the UK (24), have previously shown that tests for bTB interfere with MAP ELISA diagnostics. Varges et al. (21) reported ELISA interference occurring between 30 and 90 days post-PPD administration in bTB and MAP negative cattle. Of the 63 animals included in that study, 5 were classified as MAP ELISA positive post-PPD administration using both SICCT and single intradermal tuberculin test. Although the current study highlights a similar trend, the timescale over which interference is recorded differs. The increase in the number of animals detected ELISA positive post-SICCT and subsequent decrease to pre-SICCT prevalence occurred approximately 2 weeks earlier than the period of interference outlined by Varges et al. (21). The herd included in the

To improve visualization of interquartile ranges, only S/P values <150 shown.


**Table 2 | Multivariable GEE analysis of milk ELISA as a continuous (S/P % ELISA readings) and categorical (milk ELISA MAP positive/negative) dependent variable and independent variables**.

<sup>a</sup>May is the ELISA sample taken pre SICCT. <sup>b</sup>Parity 1: 1st lactation. No significant interactions identified with other independent variables. C.I., confidence interval. Coefficient, difference across the sample population. Statistically significant P values highlighted in bold.

#### **Table 3 | Multivariable GEE analysis of serum ELISA as a continuous (S/P % ELISA readings) and categorical (serum ELISA MAP positive/negative) dependent variable and independent variables**.


<sup>a</sup>May is the ELISA sample taken pre SICCT test.

<sup>b</sup>Parity 1 – 1st lactation. No significant interactions identified between independent variables. C.I., confidence Interval. Coefficient, difference across the sample population. Statistically significant P values highlighted in bold.

current study had a history of recording serum MAP ELISA positive individuals (within-herd prevalence of 8%). This contrasts with the Brazilian study where cattle were confirmed MAP fecal culture negative prior to inclusion in the trial. It is possible, therefore, that cows used in the current study had been pre-sensitized to MAP or additional mycobacterial-related antigens. This being the case, it would be expected that a more rapid immune response would result, i.e., a secondary humoral memory response (25). The longer duration taken to record an IgG response and the lower proportion of ELISA positive cows identified post-PPD administration by Varges et al. (21) may be indicative of a slower primary immune response (**Figure 4**). As mentioned previously, Irish cattle are tested annually using SICCT from the age of 6-weeks, which may also account for the suggested memory response in Irish cattle in contrast to their Brazilian counterparts. Additionally, the studies differed in the ELISA kits used for MAP antibody detection and used limited sample populations. More extensive studies are, therefore, required to compare the performance of all commercially available MAP ELISA kits with regard to administration of both aPPD and bPPD and the need for development of more specific antigens to improve the specificity and sensitivity of currently available MAP ELISAs has been clearly highlighted. The inclusion of a greater number, and diversity, of animals and herds would also strengthen findings, as would continuation of a study over a number of years incorporating multiple TB tests.

Varges et al. (21) examined both the single intradermal and comparative bTB test. Interestingly, cross-reacting antibodies were detected using both SICCT and single intradermal test, while administration of aPPD alone did not elicit cross-reacting antibodies. This would suggest that bPPD may be responsible for generation of cross-reacting antibodies in the MAP ELISA kits examined in both studies. This is supported by a study by Olsen et al. (18), which highlighted reduced MAP ELISA specificity in animals experimentally infected with bTB. Interestingly, animals with natural bTB infection did not elicit cross-reacting antibodies (18), which may again suggest that the intradermal administration of bPPD, is indeed, the stimulant for generation of cross-reacting antibodies. Commercially available MAP ELISA kits incorporate an *M. pheli* absorption step to increase the specificity of the assay. Again, Olsen et al. (18) showed this to be an ineffective method of improving MAP ELISA specificity with regard to bTB. It may be that while pre-absorption with *M. pheli* is somewhat successful in reducing binding of *M. avium* antibodies, repeated administration of bPPD negates its effect in preventing non-specific binding. The potential for a cumulative effect of PPD administration (either avian or bovine) from multiple bTB tests over a number of years, therefore, requires thorough investigation to fully characterize the impact of SICCT on MAP ELISA testing.

In Ireland, herds restricted due to a positive bTB diagnosis (Directive 64/432/EEC), undergo two repeat tests at a 60-day interval. For herds operating under these restrictions, the results of the current study highlight that milk samples may be a more suitable test matrix than serum ELISA to avoid test interference. Similar to the results obtained by Lombard et al. (26), there was moderate agreement between serum and milk samples. Milk samples, however, took a shorter interval to return to pre-SICCT levels

than serum in the current study. This may reflect the difference in IgG sub-classes between serum and milk and a lower milk IgG response (26). The post-SICCT period of elevated milk *S*/*P* ratios, however, may reflect a period of increased IgG production or IgG secretion from plasma to the mammary gland post-SICCT. This manifests as increased test sensitivity, stronger correlations highlighted between milk and serum results post-PPD administration. May et al. (24) also recorded significantly higher milk ELISA readings 4.5 weeks post-PPD administration in a UK herd. The limited statistical analysis completed by May et al. (24) and use of only a single testing timepoint post-SICCT presents difficulties in allowing direct comparisons between both datasets. Additionally, a number of regions in the UK administer the SICCT on one occasion every 4 years (27), a much longer testing interval than experienced by Irish herds. The differences highlighted between Varges et al., (21), May et al., (24) and the current study highlight the usefulness of examining the impact of SICCT on MAP ELISA results in multiple jurisdictions in order to more fully elucidate the impact of bTB testing on MAP diagnosis by ELISA.

It has previously been reported that exposure to environmental mycobacteria may yield low level protection against *M. tuberculosis* (28, 29). Hope et al. (30) also reported protection against *M. bovis* following exposure to *M. avium,* and that pre-exposure to *M. avium* results in an imprinting of memory against avian antigens onto T-lymphocytes. An amnestic response to environmental mycobacterial infection combined with continuous boosting of T cells in response to administration of PPD may, therefore, have the potential to assist in control of MAP at the animal level. In that regard, Ireland records a relatively low prevalence of MAP compared to additional milk exporting nations (31). For example, a total of 232 clinical cases of JD were reported in Ireland from 1995 to 2002 (32), yielding an average annual rate of approximately 0.0005%, given a cattle population of six million cattle (33). Additionally, Good et al. (34) reported that 20% of Irish herds contain at least one ELISA positive animal, again a relatively low prevalence (31). Given that environmental conditions in Ireland are conducive to the growth of mycobacteria (35), and that Irish farmers engage in high risk management practices with regard to spread of JD, e.g., widespread pooling of colostrum and milk for calf-feeding (36) (Kennedy et al. unpublished data), a higher prevalence of clinical cases and MAP ELISA positives might be expected. Another Irish study (37) recorded no significant associations between MAP seropositivity and milk production parameters, again contrasting with international studies (38, 39). It is our hypothesis that repeated annual administration of aPPD and bPPD in Ireland may induce a protective effect against MAP thereby lessening the clinical manifestations of MAP infection and resultant production losses. To more thoroughly investigate this hypothesis, it is necessary to complete in depth investigations as to whether the increase in antibody levels recorded post-PPD administration in the current study equates to an increased T-cell response, which would be required to achieve such a protective effect (40).

An advantage of the current study was the use of a compact spring-calving herd. This ensures that all cows examined were at a similar stage of lactation and physiological status. This allowed trends in MAP *S*/*P* % ratios over the latter half of lactation in a homogenous population to be examined. In agreement with a previous study (26), cows in late lactation were more likely to yield a MAP ELISA positive result using milk samples. The declining milk yields in late lactation result in a lessening of the dilution effect on antibody levels thereby increasing antibody concentrations (41). Interestingly, an increase in the prevalence of serum ELISA positives was also recorded in late lactation. This finding is in agreement with a Danish study (42). The increase in prevalence of serum ELISA positives in the current study corresponds with housing, which may increase the likelihood of exposure to mycobacterial antigens by increasing the potential for fecal contact. Nielsen et al. (42) also showed parity 2 and greater to be more likely to test ELISA positive relative to parity 1 cows, which is also highlighted in the current study. Parity 3 and 4 animals, however, were in general less likely to test positive than parity 2. The majority of Irish farmers target compact calving seasons (43) and strict culling practices are often in place (33). These culling practices may lead to less ELISA positive animals remaining in the system post second lactation. Results from this study indicate that age of animal at sampling and timing of JD ELISA tests relative to stage of lactation and time of bTB testing are important considerations when interpreting ELISA results.

#### **CONCLUSION**

Administration of PPD as part of the bTB test corresponds to an increased prevalence of ELISA positives for JD. Diagnostic sampling for JD utilizing milk ELISA should be avoided in the 43-day period following the bTB test, with serum ELISA sampling not recommended for an additional 28 days. Based on the increase in antibody titers in MAP ELISA recorded post-PPD administration, it is our hypothesis that repeated annual administration of aPPD and bPPD may induce a protective effect helping to curtail the clinical manifestations of MAP infection.

#### **ACKNOWLEDGMENTS**

This Teagasc research was funded by the Irish Dairy Levy. The authors wish to acknowledge farm personnel.

#### **REFERENCES**


**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.

*Received: 30 July 2014; accepted: 22 October 2014; published online: 12 November 2014.*

*Citation: Kennedy AE, Da Silva AT, Byrne N, Govender R, MacSharry J, O'Mahony J and Sayers RG (2014) The single intradermal cervical comparative test interferes with Johne's disease ELISA diagnostics. Front. Immunol. 5:564. doi: 10.3389/fimmu.2014.00564*

*This article was submitted to Molecular Innate Immunity, a section of the journal Frontiers in Immunology.*

*Copyright © 2014 Kennedy, Da Silva, Byrne, Govender, MacSharry, O'Mahony and Sayers. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

## Regulation of toll-like receptors-mediated inflammation by immunobiotics in bovine intestinal epitheliocytes: role of signaling pathways and negative regulators

### **Julio Villena1,2, Hisashi Aso<sup>3</sup> and Haruki Kitazawa<sup>4</sup>\***

1 Immunobiotics Research Group, Tucuman, Argentina

<sup>2</sup> Laboratory of Immunobiotechnology, Reference Centre for Lactobacilli (CERELA-CONICET), Tucuman, Argentina

<sup>3</sup> Cell Biology Laboratory, Graduate School of Agricultural Science, Tohoku University, Sendai, Japan

<sup>4</sup> Food and Feed Immunology Group, Laboratory of Animal Products Chemistry, Graduate School of Agricultural Science, Tohoku University, Sendai, Japan

#### **Edited by:**

Kieran G. Meade, Teagasc – The Agriculture and Food Development Authority, Ireland

#### **Reviewed by:**

Dirk Werling, Royal Veterinary College, UK Takao Mukai, Kitasato University, Japan

#### **\*Correspondence:**

Haruki Kitazawa, Food and Feed Immunology Group, Laboratory of Animal Products Chemistry, Department of Science of Food Function and Health, Division of Bioscience and Biotechnology for Future Bioindustries, Graduate School of Agricultural Science, Tohoku University, 1-1 Tsutsumidori-Amamiyamachi, Aobaku, Sendai 981-8555, Japan e-mail: haruki@bios.tohoku.ac.jp

Intestinal epithelial cells (IECs) detect bacterial and viral associated molecular patterns via germline-encoded pattern-recognition receptors (PRRs) and are responsible for maintaining immune tolerance to the communities of resident commensal bacteria while being also capable to mount immune responses against pathogens.Toll-like receptors (TLRs) are a major class of PRRs expressed on IECs and immune cells, which are involved in the induction of both tolerance and inflammation. In the last decade, experimental and clinical evidence was generated to support the application of probiotics with immunoregulatory capacities (immunobiotics) for the prevention and treatment of several gastrointestinal inflammatory disorders in whichTLRs exert a significant role.The majority of these studies were performed in mouse and human cell lines, and despite the growing interest in the bovine immune system due to the economic importance of cattle as livestock, only few studies have been conducted on cattle. In this regard, our group has established a bovine intestinal epithelial (BIE) cell line originally derived from fetal bovine intestinal epitheliocytes and used this cell line to evaluate the impact of immunobiotics in TLR-mediated inflammation. This review aims to summarize the current knowledge of the beneficial effects of immunobiotics in the regulation of intestinal inflammation/infection in cattle. Especially, we discuss the role of TLRs and their negative regulators in both the inflammatory response and the beneficial effects of immunobiotics in bovine IECs. This review article emphasizes the cellular and molecular interactions of immunobiotics with BIE cells through TLRs and gives the scientific basis for the development of immunomodulatory feed for bovine healthy development.

**Keywords: immunobiotics,TLR4, intestinal immunity, inflammation, bovine intestinal epitheliocytes,TLR negative regulators, lactobacilli, bifidobacteria**

### **INTRODUCTION**

Intestinal epithelial cells (IECs) are structurally and functionally polarized. These cells have an apical surface facing the intestinal lumen and a basolateral surface facing the underlying basement membrane and the lamina propria. IECs provide a physical barrier that separates commensal bacteria in the lumen from the underlying lamina propria and deeper intestinal layers (1). In addition, IECs are a central component of the immune system of the gut. Over the last decades, great progress has been achieved in understanding IECs immunobiology (2). It was amply demonstrated that the cross-talk between the epithelium with gut microbes significantly influences the activities of immune cells in the mucosa (2). The detection of commensal bacteria, pathogens, or probiotics by IECs is achieved through the families of germlineencoded pattern-recognition receptors (PRRs) that recognize conserved molecular structures known as microbe-associated molecular patterns (MAMPs). MAMPs and PRRs interaction and the subsequent signaling in IECs is involved in several important

mechanisms that are crucial for maintaining a healthy epithelial barrier including maintenance of tight junctions strength, epithelial cell proliferation and renewal, expression of antimicrobial peptides, and modulation of mucosal immune responses (3).

In recent years, worldwide interest has rapidly and significantly increased in the therapeutic and preventive effects of "friendly bacteria." These microorganisms, recognized as probiotics, are generally selected from *Lactobacilli* or *Bifidobacteria* strains (4). Several studies in animal models as well as clinical trials support a unique role for probiotics by beneficially modulating the mucosal immune system. Thus, a new term was required to identify probiotic bacteria that promote health by regulating the mucosal immune system. Clancy suggested the new term "immunobiotics" as appropriate for fulfilling this need (5). The quest for a better understanding of how immunobiotics works have led to an enormous interest in the molecular processes underlying host– microbe interactions. As reviewed by Lebeer et al. (6), the final conclusion of works that have studied the molecular mechanism

of probiotic immunomodulatory activities is that: "*their effect depends on the combination of distinct MAMPs that interact with various PRRs and the associated co-receptors that fine tune signaling, as well as on the quantity and quality of these MAMPs. Therefore, host-immunobiotic interactions are not univocal but involve complex interactions among various microbial molecules, host receptors, and adaptor molecules*" (6).

This expanding knowledge about the cellular and molecular effects of beneficial bacteria in innate mucosal immune system has raised the possibility of new treatments for improving health not only in humans but also in animals. In this review, we describe the recent advances in the impact of immunobiotics on bovine intestinal epithelial (BIE) cells and possible novel therapeutic approaches to beneficially modulate bovine epithelial cell immunobiology. Especially, we discuss the role of toll-like receptors (TLRs), their signaling pathways, and their negative regulators in both the inflammatory-intestinal injury and the beneficial effects of immunobiotics. This article emphasizes the cellular and molecular interactions of immunobiotics with bovine epithelial cells through TLRs and gives the scientific basis for the future development of immunomodulatory feed for improving bovine health.

#### **BOVINE INTESTINAL EPITHELIAL CELLS**

The development of bovine intestinal cell cultures and their characterization with regard to their permissiveness for bacterial adhesion and invasion, and the ability to sense PAMPs through PRRs represents an important step forward toward the establishment of *in vitro* systems to study molecular interactions of pathogenic, commensal, and probiotic microorganisms with the bovine host (7).

For cattle, primary cultures of ileum or colon epithelial cells have been used for toxicological assays, the study of microbial virulence factors, the efficacy of antimicrobial compounds, and the evaluation of innate immune responses through PRRs signaling (7–12). A combination of the enzymatic digestion (dispase and collagenase) together with soft mechanical agitation proved to be a successful method for releasing intact, viable bovine colonic crypts from underlying mesenchymal tissue. However, a series of purification steps was required to eliminate the majority of contaminating non-epithelial cells (mostly fibroblasts) from the crypt suspension (10). Similarly, cultures of bovine colonocytes and jejunocytes were obtained by Rusu et al. (11), using a combination of enzymatic and mechanical disruption of the intestinal epithelium. The study showed that primary cultures of bovine enterocytes isolated from colon and jejunum presented characteristics of epithelial cells, such as a typical pavement-like aspect, the formation of domes and apical tight junctions, and microvilli in confluent cultures. Moreover, these bovine colonocytes and jejunocytes expressed epithelial cell markers such as brush border enzymes and the epithelium typical cytoskeleton proteins, such as cytokeratins.

Dibb-Fuller et al. (8) developed primary bovine cell lines from ileum, colon, and rectum and those bovine primary gastrointestinal epithelial-derived cells were successfully used to assess adherence and invasion of several intestinal pathogenic bacteria including enterohemorrhagic *Escherichia coli* (EHEC) and *Salmonella enterica* serotype typhimurium. In addition, bovine colonic crypts cells isolated and purified from the mucosa, proved to be useful *in vitro* tools to study virulence factors of EHEC, verotoxins, in particular (10). The work showed the expression of globotriaosylceramide Gb3 by bovine colonocytes, which directly contrasts with the absence of this receptor on human intestinal epithelium. This fact represents a fundamental difference that could have major significance in the different pathogenicity of EHEC in these hosts, and make bovine colonocytes an invaluable tool for studying EHEC infection (10).

Bridger et al. (7) showed that epithelioid cells from bovine colonic crypts formed a confluent monolayer on the surface of collagen-coated culture flasks. Those bovine colonocytes expressed epithelial cell-specific cytokeratin and cell membraneassociated tight junctional ZO-1 in the contact area between neighboring cells. Semi-quantitative RT-PCR demonstrated variable amounts of gene transcripts for different TLR genes. Notably, primary bovine colonocytes expressed TLR4 mRNA while transcripts for TLR1, TLR3, and TLR6 were also detectable in some cultures. Moreover,the study showed that colonocytes significantly up-regulated the expression of IL-8, MCP-1, and RANTES when challenged with pathogenic *E. coli* or lipopolysaccharide (LPS) (7). Therefore, short-term bovine colonocytes cultures proved to be suitable *in vitro* models to study pathogens-specific responses of the bovine colonic mucosa.

Those studies clearly demonstrated that bovine primary IEC cultures represent valuable tools to assess the molecular mechanisms involved in pathologies caused by infectious agents. However, the cellular and molecular interactions of commensal or probiotic bacteria with BIE cells have been less explored.

There is an increasing research in the use of immunobiotics to beneficially modulate the mucosal immune system in animals. Immunobiotic bacteria could be used for improving resistance against pathogens and decreasing intestinal inflammatorymediated tissue damage (13–17). In this regard, our laboratories and others' have conducted *in vitro* and *in vivo* studies utilizing different lactobacilli and bifidobacteria strains to evaluate the effect of immunobiotics against infections and inflammation in animals, however, the majority of these studies were performed in swine and only few in the cattle (13).

We hypothesized that a cell line of BIE cells could be a useful *in vitro* model for the evaluation of the molecular interactions between IECs and bovine pathogens. In addition, such system would allow the selection of immunobiotic microorganisms and the study of the mechanisms by which probiotic lactic acid bacteria (LAB) functionally modulate bovine IECs. Therefore, in order to: (a) understand the functional role of IECs in bovine mucosal host defense and (b) select potential immunobiotic LAB strains that may be used to beneficially modulate the inflammatory response in bovine IECs; we have recently established an immortalized BIE cell line (18). Monolayer cobblestone and epithelial-like morphology is assumed by BIE cells when cultured, with close contact between the cells. Moreover, scanning electron microscopy examination of BIE cell reveled that 3-days old cells have irregular and slender microvilli-like structures on their surface and that this structures increase in complexity as the cells grow (18). BIE cells undergoes over at least 40 passages with no detectable loss of epithelial properties.

Expression of TLRs mRNA in BIE cells was evaluated by RTqPCR and it was demonstrated that all TLRs genes were expressed in BIE cells (19). TLR1, 3, 4, and 6 were strongly expressed, followed by TLR5, 8, 9, 10, 2, and 7. We were particularly interested in expression of TLR2 and TLR4 in BIE cells as the main receptors detecting Gram(+) probiotic bacteria and Gram(−) pathogens, respectively. Therefore, to confirm real-time PCR findings, we further examined the expression of TLR2 and 4 proteins in BIE cells using anti-TLRs antibodies that are able to cross-react with bovine TLRs (19). Visualization of the immunofluorescence staining confirmed the protein expression of TLR2 and 4 in BIE cells (**Figure 1A**).

**FIGURE 1 |Toll-like receptor 4 (TLR4) expression in bovine intestinal epithelial epitheliocytes (BIE cells) (A)**. Upon recognition of lipopolysaccharide (LPS), TLR4 dimerizes and initiates a signaling cascades that include phosphatidylinositol 3-kinase (PI3K), mitogen-activated protein kinase (MAPK), and nuclear factor κB (NF-κB) pathways; and culminates in the production of inflammatory mediators by BIE cells **(B)**.

The inflammatory response triggered by BIE cells in the face of a challenge with heat-stable enterotoxigenic *E. coli* (ETEC) MAMPs was also evaluated. Upon pathogen binding to TLR4 complex, this receptor recruits, through its short intracellular toll-interleukin-1 receptor (TIR) domain, adaptor molecules, and kinases, thus initiating a downstream signaling cascade that culminates in the production and secretion of inflammatory mediators such as TNFα, IL-1β, IL-6, and IL-8 (**Figure 1B**). The ETEC 987P strain used in our experiments does not express the TLR5-ligand flagellin, and we demonstrated that the main molecule responsible for the inflammatory response triggered by this bacterium is the LPS (15, 16). Stimulation of BIE cells with heat-stable ETEC MAMPs from strain 987P enhanced the production of the pro-inflammatory cytokines IL-6, IL-8, IL-1β, and MCP-1 by activating mitogenactivated protein kinase (MAPK) and nuclear factor κB (NF-κB) pathways (19). These findings are in line with our previous reports demonstrating that the heat-killed ETEC 987P strain triggers a TLR4-mediated inflammatory response in porcine intestinal epithelial (PIE) cells through NF-κB and MAPK pathways (20). In addition, our results in BIE cells correlate with studies of the immune response against ETEC in IECs of different hosts' species. It was shown that both NF-κB and MAPK pathways are important mediators of ETEC and LPS activation in human (HT-29 and T84) and mouse (CMT93) IECs (15, 21).

Available lines of evidence indicate that bovine epithelial cells, including intestinal, mammary, bronchial, and nasopharynx epithelial cells respond to bacterial LPS and other microbial products by producing pro-inflammatory cytokines required to combat invading pathogens. Therefore, the pro-inflammatory mediators produced by BIE cells in response to ETEC may have an important protective role during the course of intestinal infections. The chemokine IL-8 stimulates a strong infiltration of neutrophils in the gut lamina propria, a fact that is consistently observed upon ETEC infection. After IL-8 induced recruitment of neutrophils, increase of IL-6 production is able to induce degranulation of these cells, thereby enhancing the inflammatory response (22). In addition, IECs are able to produce MCP-1 in response to ETEC challenge. This chemokine has potent monocyte-activating and attracting properties and plays a major role during intestinal inflammation (23). Therefore, BIE cells respond to the presences of ETEC and LPS by activating the TLR4-signaling pathway, which is necessary to initiate a robust defensive action against intruders.

Our studies indicate then that the BIE cell line could be a useful cell line for evaluating inflammatory responses via TLR4 *in vitro*. Furthermore, considering that inflammatory responses induced by intestinal pathogens can lead to dysregulation of IECs signaling, disruption of membrane barrier integrity, enhancement of pathogen translocation and disease (24), BIE cells could be also used to evaluate therapies designed for preventing inflammatory damage caused by bovine intestinal pathogens or their associated PAMPs or virulence factors.

#### **PROBIOTICS FOR THE BOVINE HOST**

Several studies on the pathogenesis of intestinal inflammation/infection both in man and experimental animals continue to show the importance of commensal bacteria in the gastrointestinal tract in stimulating and directing the immune system. Moreover, the ability of immunobiotic bacteria to beneficially modulate the response against intestinal pathogens in animals through the improvement of resistance and the reduction of inflammatorymediated tissue damage has been described by several reports (25–27). Before weaning, dairy calves are highly susceptible to several pathogens. For several years, antibiotics have been used to overcome these problems also to obtain economic benefits in terms of improved calves performance and reduced medication costs. However, the use of antibiotics in animal husbandry is in question because of antibiotic resistance of microorganisms. In an effort to replace antibiotics from bovine feeds, many additives have been proposed including the use of probiotics (28, 29). In fact, some few studies have shown that probiotic bacteria can be used as growth promoters in calves instead of antibiotics to counteract the negative effects of their widespread use (30). Early studies of Abe et al. (31) showed that oral administration of *Bifidobacterium pseudolongum* or *Lactobacillus acidophilus* to calves improved body weigh gain and decreased the frequency of diarrhea occurrence. Similarly, Mokhber-Dezfouli et al. (32) demonstrated that probiotic treatments have the ability to beneficially modulate body weight gain, body height, and general health condition of calves. Additionally, oral treatment with probiotic *E. coli* significantly reduced the pathogenicity and fecal shedding of EHEC in calves (33, 34). It was also reported that a mixture composed of *Lactobacillus casei* DSPV 318T, *Lactobacillus salivarius* DSPV 315T, and *Pediococcus acidilactici* DSPV 006T protected calves against *Salmonella* Dublin infection (35). These studies clearly show the potential of probiotic bacteria to beneficially modulate gastrointestinal hemostasis in the bovine host. However, the cellular and molecular mechanisms involved in the probiotic activities in cattle have not been studied in depth.

#### **REGULATION OF INFLAMMATION IN BIE CELLS BY IMMUNOBIOTIC LACTOBACILLI**

The first contact of immunobiotic bacteria with the intestinal mucosa is mediated by the single cell layer of IECs. As mentioned before, these IECs are of paramount importance in host– immunobiotic cross-talk. Then, we thus sought to determine whether an immunobiotic *Lactobacillus* strain could regulate the inflammatory response induced by heat-stable ETEC MAMPs in BIE cells. Our previous studies with the strain *Lactobacillus jensenii* TL2937 showed that this bacterium has remarkable immunomodulatory effects in porcine IECs and immune cells [for a review see Ref. (3)]. The TL2937 strain is able to functionally modulate porcine IECs by inhibiting excessive MAPK- and NF-κBinduced pro-inflammatory cytokine production (IL-6 and IL-8) in response to TLR4 activation (3, 15). Consequently, we first focused on *L. jensenii* TL2937 to evaluate its anti-inflammatory effect in BIE cells. Preincubation of BIE cells with *L. jensenii* TL2937 significantly decreased IL-6 and IL-8 expressions in 20 and 25% with respect to the control, respectively, after heat-stable ETEC MAMPs challenge (19). However, this effect was lower when compared with the anti-inflammatory activity of this strain in PIE cells (15). In porcine, IECs previously treated with the TL2937 strain, stimulation with heat-stable ETEC MAMPs reduced IL-6 and IL-8 expressions by 35 and 30% when compared to control cells, respectively (15). Although the effect of the *L. jensenii* TL2937 in BIE cells was lower than the one previously reported for porcine IECs, our first studies in BIE cells indicated that probiotic lactobacilli could be beneficial for attenuating inflammatory damage caused by TLR4 activation in bovine epithelial cells (19). Thus, we next aimed to screen and select the most effective immunoregulatory lactobacilli strains able to modulate TLR4-mediated pro-inflammatory response in BIE cells. Several lactobacilli strains were evaluated in our bovine IECs line and we found that some of these bacteria were capable to downregulate the expression of inflammatory cytokines. Among these strains, *L. casei* OLL2768 showed the most pronounced effect (19). Notably, the anti-inflammatory activity of the OLL2768 strain was more pronounced than that observed for *L. jensenii* TL2937 in BIE cells, while the effect of OLL2768 strain was lower in PIE cells (15). It is well known that probiotic activities are strain specific. In addition, our findings clearly indicated that is necessary to carefully evaluate different strains according to the specific host, because the effect of the same *Lactobacillus* may be different according to the host that consumes it. Then, our *in vitro* bovine system could be of great value to find potential immunobiotic strains suitable for the improvement of the bovine host health.

We also aimed to define the molecular mechanisms by which *L. casei* OLL2768 attenuated heat-stable ETEC MAMPs-induced pro-inflammatory response in BIE cells. Our data showed that the immunoregulatory effect was related to the capacity of OLL2768 strain to inhibit NF-κB and MAPK p38 signaling pathways in bovine IECs after TLR4 activation.

Nuclear factor κB is composed of several protein subunits regulating DNA transcription. Under non-stimulatory conditions, it is bound to the inhibitor molecule IkB in the cytoplasm. After TLR activation IkB is phosphorylated by IKK and once freed from IkB, NF-κB subunit p65 (RelA) migrates into the nucleus, where it binds to target promoters and activates transcription of effector genes including TNF-α, IL-8, and others (36, 37). Among many up-stream signaling proteins involved in NF-κB activation, TLR4 plays a critical role and it is well-documented that TLR4/NF-κB pathway has a pivotal role in the pathogenesis of several intestinal inflammatory diseases and infection-induced tissue damage (38). Some studies have reported the ability of probiotic lactobacilli to modulate TLR4/NF-κB pathway in the gut (39, 40). It was showed that *Lactobacillus suntoryeus*, a gut commensal, blocks inflammatory mediators (Cox2, TNF-α, IL-1, and IL-6) through suppression of TLR4-linked NF-κB activation in mice with 2,4,6 trinitrobenzene sulfonic acid-induced colitis (41). Liu et al. (42) reported that *Lactobacillus reuteri* strains DSM 17938 and ATC-CPTA4659 led to decrease intestinal protein levels of TLR4 and decreased pro-inflammatory cytokine levels in parallel with inhibition of TLR4-signaling via the NF-κB pathway in newborn rats with necrotizing enterocolitis. In addition, it was reported that some probiotics strains are able to suppress TNF- or *Salmonella typhimurium*-induced IL-8 gene expression and secretion by IECs in a NF-κB-dependent manner (39, 40). Our experiments also demonstrated that *L. casei* OLL2768 is able to inhibit TLR4/p38 signaling pathway since we demonstrated that in lactobacillitreated BIE cells the phosphorylation of p38 was reduced after challenge in heat-stable ETEC MAMPs (19). Regulation of MAPK p38 pathway by probiotics has been described before. *L. rhamnosus* GG was found to significantly down-regulate expression of p38 in human monocyte-derived DCs after the challenge with LPS (43). It was also showed that *Lactobacillus bulgaricus* LBG, inhibited the activation of the TLR4-signaling pathway and IL-8 production induced by *Helicobacter pylori* LPS in human gastric adenocarcinoma cells through blocking MAPK p38 (44). Then, our findings in BIE cells are reminiscent of other studies showing that probiotic *L. casei* OLL2768 is capable of modulating TLR4/NF-κB- and TLR4/p38-induced inflammation.

The JNK and p38 MAPK pathways share several up-stream regulators, and accordingly there are multiple stimuli that simultaneously activate both pathways. In this regard, it was showed that the conditioned media from probiotic *L. rhamnosus* GG induced the expression of cellular heat shock protein (Hsp72) in IECs in a p38- and JNK-dependent manner (45). The work showed that *L. rhamnosus* GG conditioned media treatment resulted in a clear activation of both p38 and JNK pathways in IECs. Moreover, exposure of IECs to inhibitors against p38 and JNK before conditioned media treatment resulted in blockade of Hsp72 expression, thus, confirming a likely role for both MAPK signaling pathways in the probiotic effect (45). Then, we expected that *L. casei* OLL2768 had an inhibitory effect on JNK pathway in BIE cells as observed for the p38 MAPK pathway. However, increased levels of p-JNK were detected in BIE cells stimulated with the OLL2768 strain. It was also reported that JNK and p38 MAPK pathways may induce opposite effects. In fact, there is evidence indicating that the p38 MAPK pathway can negatively regulate JNK activity in several contexts (46, 47). The first evidence of the interaction of these two pathways was the observation that inhibition of p38 strongly increased the activation of JNK (46). The work analyzed the effect of specific p38 MAPK inhibitors, SB202190 or SB203580, on JNK phosphorylation in A549 human lung alveolar epithelial cells, and found that inhibition of p38 MAPK could induce JNK activation due to a compensatory mechanism (46). In addition, it was showed that the p38 inhibitor SB203580 enhances the activation of JNK isoforms after the challenge of IECs with IL-1β or by LPS in macrophages (46). In line with those studies, the kinetic analysis of p38 and JNK phosphorylation in BIE cells showed an early upregulation of p–p38 between 5 and 10 min after heat-stable ETEC MAMPs challenge that was followed by a down-regulation of p-JNK between 10 and 20 min (19). Therefore, we can speculate that *L. casei* OLL2768 has a direct influence in p38 pathway while its effect in JNK is the result of the inhibition of p38 phosphorylation. Further research is needed to clarify completely the influence of *L. casei* OLL2768 in MAPK pathways in BIE cells.

Following TLR activation, there must be a checkpoint where TLR signaling is abolished and the system is returned to a normal physiological state to avoid a harmful response toward the host immune system. Regulatory factors able to modulate the duration and intensity of TLRs signals are therefore key components for the protection of the hosts (48). Several regulatory mechanisms have been described for TLRs including soluble and decoy factors, membrane-associated protein regulators, negative regulators of the adaptor complex, and microRNA (3). To assess the expression of these negative regulators of TLRs, we first cloned cDNAs corresponding to these proteins in the bovine (19). We demonstrated the expression of A20-binding inhibitor of nuclear factor kappa B

activation 3 (ABIN-3); B-cell lymphoma 3-encoded protein (Bcl-3); interleukin-1 receptor-associated kinase M (IRAK-M); single immunoglobulin IL-1-related receptor (SIGIRR); toll interacting protein (Tollip); and mitogen-activated protein kinase 1 (MKP-1) in BIE cells (19). Consequently, the effect of *L. casei* OLL2768 on the expression of these negative regulators of the TLRs signaling was next evaluated. *L. casei* OLL2768 is able to up-regulate Tollip and Bcl-3 in BIE cells, and in this way to negatively regulate TLR4 signaling (19). Tollip is able to suppress the activity of IL-1 receptor-associated kinase (IRAK), and inhibit TLR4-triggered NF-κB and MAPK signaling pathways (49, 50). It was showed that stimulation of IECs with a TLR ligand, such as LPS, induces a state of hyporesponsiveness through up-regulation of Tollip that limits pro-inflammatory responses triggered by a second challenge with the same or another TLR ligand (50). The expression of Tollip was reported in bovine mammary epithelial cells (51). The work showed that expression of Tollip was increased in response to LPS, suggesting that the bovine mammary epithelium possesses the necessary immune repertoires required to regulate TLR4 activation. Recently, it was demonstrated by Fu et al. (52) that Tollip is significantly up-regulated in bovine endometrial epithelial cells after the stimulation with LPS, and that this up-regulation of Tollip was necessary for the regulation of the overexpression of NF-κB and the protection against the inflammatory damage. On the other hand, the Bcl-3 is a nuclear protein and member of the NF-κB family. Bcl-3 is able to stabilize repressive NF-κB homodimers in a DNA-bound state, and in this way prevents the binding of transcriptionally active dimers. Therefore, Bcl-3 functions as an inhibitor of NF-κB activity. In recent years, a role of Bcl-3 has been revealed in LPS tolerance via its ability to stabilize the p50 homodimer, and thus, has been identified as a negative regulator of TLR4 signaling (53). Furthermore, by selectively affecting chromatin remodeling, Bcl-3 mediates repression of pro-inflammatory genes, and also facilitates the expression of the anti-inflammatory gene *IL-10* (54). Therefore, the induction of Bcl-3 and Tollip by *L. casei* OLL2768 in BIE cells is important in establishing tolerance against heat-stable ETEC MAMPs.

It is not possible to give a precise molecular mechanism for the anti-inflammatory action of *L. casei* OLL2768 on BIE cells at present. However, it could be hypothesized that interaction of *L. casei* OLL2768 with BIE cells through one or more PRRs induces the up-regulation of the negative regulators Bcl-3 and Tollip, which reduce the production of inflammatory mediators in response to heat-stable ETEC MAMPs (**Figure 2**). One of the possible PRR involved in *L. casei* OLL2768 immunoregulatory capacities could be TLR2. Studies with the TLR2 ligand Pam3CSK4 in BIE cells demonstrated that the treatment with the TLR2 agonist up-regulate the expression of Tollip and reduce activation of NF-κB and p38 MAPK pathways (19). However, further research is needed to resolve which PRR is activated by *L. casei* OLL2768 for the induction of negative regulators.

#### **REGULATION OF INFLAMMATION IN BIE CELLS BY IMMUNOBIOTIC BIFIDOBACTERIA**

Members of the genus *Bifidobacterium* are considered to be important constituents of the microbiota of animals, from insects to mammals. They are gut commensals extensively used by the

food industry as probiotic microorganisms, since some strains have been shown to have specific beneficial effects. Bifidobacteria are able to prevent or alleviate infectious diarrhea through their effects on the immune system and resistance to colonization by pathogens. In addition, some bifidobacteria strains have potent anti-inflammatory capacities that could be used to reduce inflammatory-intestinal damage. *Bifidobacterium animalis* strain AHC7 decrease NF-κB activation in mice infected with *S. typhimurium*. *B. animalis* AHC7 consumption in this mouse model was associated with protection against inflammatory damage through modulation of secreted IL-10 and IL-12p70 and enhancement of Foxp3 expression in naïve T cells (55). In line with these results, it was showed that *Bifidobacterium bifidum* W23 was able to induce a suppression of IL-8 synthesis by Caco-2 cells challenged with *S. enterica* serovar enteritidis, and that the protective role of this probiotic strain was mediated, at least in part, via Hsp70 expression (56). Moreover, it was recently reported that *Bifidobacterium adolescentis* FRP 61, *Bifidobacterium longum* FRP 68 and FRP 69, and *Bifidobacterium breve* FRP 334 significantly reduced IL-8 production by HT-29 cells challenged with *S. typhimurium* (57). Others studies evaluating the effect of bifidobacteria in intestinal Caco-2 cells showed that *B. animalis* MB5 avoid cytokine deregulation upon ETEC challenge by inducing upregulation of IL-1β and TNF-α, and the down-regulation of TGF-β expression (58, 59). Additionally, we demonstrated in porcine IECs cells that treatment with *B. breve* MCC-117 significantly reduced the expression of inflammatory cytokines in response to heat-stable ETEC MAMPs. Moreover, studies with porcine immune cells showed that *B. breve* MCC-117 was able to reduce the levels of IFNγ in CD4<sup>+</sup> and CD8<sup>+</sup> lymphocytes and improved IL-10 levels in CD4+CD25highFoxp3<sup>+</sup> lymphocytes (14). These are among several other studies that clearly showed that bifidobacteria are highly effective in regulating pathogenic inflammation in the gut. Therefore, we next aimed to select potential immunomodulatory bifidobacteria able to beneficially modulate the inflammatory response in BIE cells.

The potential use of bifidobacteria as probiotic for cattle is supported by some new reports indicating the presence of these bacteria in young calf intestines and the fact that their presence in high numbers is associated with good health status of the host (60). Therefore, some bifidobacteria strains, previously selected in our porcine systems, were used to evaluate their anti-inflammatory capacities in heat-stable ETEC MAMPs-challenged-BIE cells (61). Similarly to the effect of *L. casei* OLL2768, some bifidobacteria were able to reduce the production of inflammatory mediators triggered by TLR4 activation in BIE cells. Considering their ability to reduce the expression of IL-6 and MCP-1, bifidobacteria strains were divided in the following two groups: (1) strains able to reduce both IL-6 and MCP-1 (*B. adolescentis* MCC-75 and *B. breve* MCC-117) and (2) strains able to reduce only IL-6 (*B. longum* BB536, *B. adolescentis* ATCC15705 and *Bifidobacterium infantis* MCC-1021) (61).

As described for the immunoregulatory *L. casei* OLL2768 strain, we also aimed to evaluate signaling pathways and TLR negative regulators expression in BIE cells after the treatment with bifidobacteria belonging to the two functional groups defined by our studies. Then, we selected *B. adolescentis* MCC-75, *B. breve* MCC-117 (strains with high anti-inflammatory capacities), and *B. adolescentis* ATCC15705 (strain with moderate anti-inflammatory capacity) for further experiments. Activation of MAPK, NF-κB, and phosphatidylinositol 3-kinase (PI3K) pathways, and changes in the expression of TLR negative regulators in MCC-75-, MCC-117- and ATCC15705-treated BIE cells were then studied. We found that each bifidobacteria strain induces unique changes in TLR4 signaling in bovine IECs (61) (**Figure 2**).

As mentioned before, several negative regulatory mechanisms control TLRs-mediated inflammatory responses and restore immune system balance in the gut. Although the NF-κBdependent gene expression is critical to the induction of an efficient immune response, excessive, or prolonged NF-κB signaling can contribute to the development of several inflammatory diseases. Therefore, this signaling transduction pathway has to be tightly regulated by several intracellular proteins. The ubiquitinediting enzyme A20 is key regulator of the TLRs signaling. It was showed that A20 deficiency in IECs renders mice sensitive to TNFα-induced lethal inflammation (62, 63). Moreover, it was reported that A20 is an early response negative regulator of TLR4 and TLR5 signaling in IECs that functions during intestinal inflammation to control the innate immune system (64). In addition, the A20 binding inhibitor of NF-κB activation (ABIN) is LPS-inducible proteins that negatively regulate NF-κB activation in response to TNF-α and LPS (65). ABINs have been described as three different proteins (ABIN-1, -2, and -3) that bind A20. Overexpression of ABINs inhibits NF-κB activation by TNF-α and several other stimuli. Similar to A20, ABIN-3 expression is NF-κB dependent, implicating a potential role for the A20/ABIN complex in the negative feedback regulation of NF-κB activation (66). Therefore, the induction of A20/ABIN complex by bifidobacteria in BIE cells is important in establishing tolerance against heat-stable ETEC MAMPs. This is in line with our previous reports in porcine IECs. In our works in the porcine systems, we showed that the bifidobacteria strains with the highest capacity to downregulate the expression of inflammatory cytokines in response to heat-stable ETEC PAMPs were also able to up-regulate A20. In fact, the most potent anti-inflammatory bacteria evaluated in our laboratory, bifidobacteria strains BB536 and M-16V and *L. jensenii* TL2937, strongly up-regulated the ubiquitin-editing enzyme A20 [for a review see Ref. (3)].

Bcl-3 protein functions as an inhibitor of NF-κB activity as mentioned before. In addition, SIGIRR, Tollip, and IRAK-M are also known to be expressed at high levels in IECs, and to thereby contribute to the hyporesponsiveness of IECs to commensals (64, 67, 68). Therefore, induction of these five negative regulators by bifidobacteria in BIE cells may be important for establishing tolerance against heat-stable ETEC MAMPs (**Figure 2**). Moreover, the fold expression increase of the negative regulators of the TLRs signaling should be also important since the levels of ABIN-3, IRAK-M, and Bcl-3 were significantly higher in *B. breve* MCC-117- and *B. adolescentis* MCC-75-treated BIE cells when compared with BIE cells treated with the moderate anti-inflammatory strain *B. adolescentis* ATCC15705 (61).

On the other hand, the MAPK pathway is involved in the upregulation of several inflammatory genes, and MKP-1 plays a role in the inhibition of pro-inflammatory mRNA expression, because it can inactivate MAPK pathway (69). Therefore, we expected that bifidobacteria with high anti-inflammatory activity significantly up-regulate MKP-1 expression and reduce MAPK activation as we have observed with other anti-inflammatory immunobiotic strains (3). However, when ERK, p38, and JNK MAPK activation and MKP-1 expression were studied in BIE cells treated with bifidobacteria, we found that *B. adolescentis* MCC-75 and *B. breve* MCC-117 activated ERK MAPK pathway and only moderately up-regulated MKP-1. On the contrary, *B. adolescentis* ATCC15705 strongly increased expression of MKP-1 and inhibit p38 and JNK pathways (61). It is known that the ERK pathway play key regulatory functions in a diverse spectrum of biological processes such as cell proliferation, differentiation, survival, and motility (70). It was also reported that TGF-β induces ERK activity in IECs and this TGF-β/ERK interaction regulates genes that are crucial for cell growth, migration, and survival of IECs (71, 72). In fact, treatment with TGF-β prevents mucosal-injury, enhances p-ERK and β-catenin, induces enterocyte proliferation, inhibits enterocyte apoptosis, and improves intestinal recovery following methotrexate-induced intestinal-mucositis in rats (73). Moreover, TGF-β increases protein levels, collagen I, TGF-β of type-1 inhibitor of plasminogen activator, and the TGF-β-converting enzyme furin in various IEC lines via ERK (74) indicating an important immunoregulatory role of the ERK pathway in maintaining homeostasis in IECs. Therefore, the activation of the ERK pathway by *B. adolescentis* MCC-75 and *B. breve* MCC-117 during ETEC MAMPs-mediated inflammation could have an important protective role against inflammatory damage (61).

Our studies also demonstrated that both *B. adolescentis* MCC-75 and *B. breve* MCC-117 were able to inhibit PI3K pathway in heat-stable ETEC MAMPs-challenged-BIE cells. It is known that PI3K regulates TLR signaling in both positive and negative ways. By mutating specific tyrosine residues in the cytosolic domain of TLR2, it was showed that there is a loss in the capacity of p85 to associate with this receptor and in the ability of TLR2 to activate NF-κB pathway. Furthermore, inhibition of PI3K during TLR2 stimulation has been shown to reduce NF-κB activation (75). On the contrary, studies in PI3K or p85a deficient mice showed that PI3K negatively regulates TLR signaling (76, 77). Then, it is well established that PI3K could affect TLRs signaling pathways in different ways and that it effect depends on cell type and readout. In our studies, we showed that stimulation of BIE cells with heat-stable ETEC MAMPs activated PI3K pathway, indicating that PI3K is positively involved in TLR4 signaling in BIE cells (**Figure 1B**). Moreover, we demonstrated that bifidobacteria able to reduce activation PI3K pathway were the strains with the highest anti-inflammatory activity (61) (**Figure 2**).

As mentioned before, our works evaluating the immunoregulatory activity of immunobiotics demonstrated that the upregulation of some regulatory cytokines and down-regulation of inflammatory mediators is dependent of TLR2 activation (15, 16). Therefore, we also investigated the role of TLR2 in the immunoregulatory effects of *B. adolescentis* MCC-75, *B. breve* MCC-117, and *B. adolescentis* ATCC15705 by using anti-TLR2

blocking antibodies (61). It was showed that the reduction of IL-6 induced by bifidobacteria in ETEC MAMPs-challenged-BIE cells was abolished when anti-TLR2 antibodies were used. This is in line with other reports conferring a key role to TLR2 in the recognition of bifidobacteria, which possess anti-inflammatory activities (78– 80). It is known that stimulation of TLR2 is able to induce tolerance against a subsequent LPS challenge (81). Therefore, it is possible that bifidobacteria could induce this type of cross-tolerance in BIE cells through their interaction with TLR2. In addition, we showed that the reduction of MCP-1 levels after challenge of BIE cells was not abolished when anti-TLR2 antibodies were used. This finding indicates that additional PRRs may be involved in the anti-inflammatory effects of *B. adolescentis* MCC-75, and *B. breve* MCC-117 in BIE cells.

#### **CONCLUSION**

The knowledge of the cellular and molecular interactions of human IECs with commensal and probiotic bacteria is rapidly progressing. An exciting possibility is that similar systems developed for the bovine host could serve as a platform for medicine and research. However, to achieve this goal, bovine IECs cultures must be enhanced and improved to allow that functional assays can be performed. Our research work has demonstrated that the BIE cell line is a useful *in vitro* tool for the study of TLR4-induced inflammatory responses in the bovine intestinal epithelium. We have also demonstrated that BIE cells could be used for a rapid screening and selection of potential immunobiotic bacteria as well as for studying the molecular mechanisms involved in their beneficial protective activity.

Despite the unique effect of each lactobacilli or bifidobacteria strain, some general conclusions can be made when comparing the effect of the two different immunoregulatory groups: high anti-inflammatory activity (*L. casei* OLL2768, *B. adolescentis* MCC-75, and *B. breve* MCC-117) and moderate antiinflammatory activity (*B. adolescentis* ATCC15705) (**Figure 2**): (1) anti-inflammatory capacity in BIE cells is strain dependent, as demonstrated by the differential effect induced by each strain, even those of the same specie (*B. adolescentis* MCC-75 and *B. adolescentis* ATCC15705). (2) The upregulation of TLR negative regulators and the intensity of that upregulation would be related to the different immunomodulatory capacity of each immunobiotic strain. Notably, upregulation of Tollip and Bcl-3 seems to be related to a high anti-inflammatory capacity. (3) The inhibition of PI3K pathway would be related to the high anti-inflammatory effect of immunobiotics in BIE cells. (4) The balance between MAPK activation and MKP-1 upregulation would be related to the anti-inflammatory effect of bifidobacteria in BIE cells. (5) The anti-inflammatory effect of immunobiotics in BIE cells is partially dependent on TLR2. Further research is needed to resolve which other PRR is involved in the immunoregulatory effects. In addition, one general conclusion can be made when comparing the effect of the two different immunoregulatory groups of bifidobacteria (**Figure 2**). (6) The upregulation of IRAK-M and ABIN-3 and the intensity of that upregulation would be related to the different immunomodulatory capacity of each bifidobacteria strain.

We believe that studies in BIE cell would provide useful information that may help in the near future to develop new functional feeds able to beneficially modulate the mucosal immune system in the bovine host. In line with this, we showed in our studies that the immunoregulatory strains *L. casei* OLL2768,*B. adolescentis* MCC-75, and *B. breve* MCC-117 are able to functionally modulate BIE cells by attenuating TLR4-induced NF-κB, MAPK, and IP3K activation and inflammatory cytokines production. Then, OLL2768, MCC-75, and MCC-117 strains would be good candidates for *in vivo* evaluation of the protective effect of immunobiotics against inflammatory damage induced by bovine intestinal pathogens or their associated PAMPs. We also believe that the immunobiotic application in cattle could contribute to produce safety animal foods via improving bovine health.

#### **ACKNOWLEDGMENTS**

This study was partially supported by a Grant-in-Aid for Scientific Research (B) (No. 24380146) and Challenging Exploratory Research (No. 26660216) from the Japan Society for the Promotion of Science (JSPS), and the Japan Dairy Association (J-milk) to Dr. Haruki Kitazawa. We thank Leonardo Albarracin for his help with the design and development of figures.

#### **REFERENCES**


porcine intestinal epithelial cell culture system: finding new anti-inflammatory immunobiotics. *FEMS Immunol Med Microbiol* (2011) **63**:129–39. doi:10.1111/ j.1574-695X.2011.00837.x


fecal shedding by treatment with probiotic *Escherichia coli*. *J Food Prot* (2003) **66**:924–30.


**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.

*Received: 11 June 2014; accepted: 19 August 2014; published online: 02 September 2014. Citation: Villena J, Aso H and Kitazawa H (2014) Regulation of toll-like receptors-mediated inflammation by immunobiotics in bovine intestinal epitheliocytes: role of signaling pathways and negative regulators. Front. Immunol. 5:421. doi: 10.3389/fimmu.2014.00421*

*This article was submitted to Molecular Innate Immunity, a section of the journal Frontiers in Immunology.*

*Copyright © 2014 Villena, Aso and Kitazawa. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

REVIEW ARTICLE published: 27 November 2014 doi: 10.3389/fimmu.2014.00611

## The role of microRNAs in bovine infection and immunity

#### **Nathan Lawless 1,2†‡, Peter Vegh1,3‡ , Cliona O'Farrelly <sup>2</sup> and David J. Lynn4,5\***


#### **Edited by:**

Uday Kishore, Brunel University, UK

#### **Reviewed by:**

Anthony George Tsolaki, Brunel University, UK Robert Braidwood Sim, University of Leicester, UK

#### **\*Correspondence:**

David J. Lynn, South Australian Health and Medical Research Institute, PO Box 11060, Adelaide, SA 5001, Australia e-mail: david.lynn@sahmri.com

#### **†Present address:**

Nathan Lawless, Compton Laboratory, The Pirbright Institute, Compton, Berkshire, UK

‡Nathan Lawless and Peter Vegh have contributed equally to this work.

#### **INTRODUCTION**

MicroRNAs (miRNAs) are short, non-coding RNAs, which posttranscriptionally regulate gene expression (1). Since their initial discovery in 1993 (2),studies have convincingly demonstrated critical roles for miRNAs in the regulation of many cellular processes, such as differentiation and proliferation (3). There is also substantial evidence, primarily from human and mouse studies, which miRNAs regulate innate and adaptive immune mechanisms (4, 5). However, the regulatory potential of miRNAs in agriculturally important species, such as cattle, is poorly explored. Here, we aim to review the role of miRNAs in bovine immune and inflammatory systems.

#### **MicroRNA BIOGENESIS**

MicroRNAs are transcribed by RNA Polymerase II or III as primary miRNAs (pri-miRNAs) in the nucleus, and are then processed into pre-miRNAs by the Microprocessor multiprotein complex, and the co-factor DiGeorge Syndrome Critical Region 8 (DGCR8/Pasha) (6). Following export to the cytoplasm by exportin 5 (XPO5), a mature 22 nucleotide long duplex is formed by an RNAse type III enzyme, called Dicer (7). The duplex, together with Trans-Activation Responsive RNA-Binding Protein 2 (TARBP2) and Argonaute (AGO) family proteins, form a complex, which triggers the assembly of the RNA-Induced Silencing Complex (RISC). One miRNA strand is removed, and the other strand guides RISC to its target mRNA via base-pairing (6). The recognition of the target binding site, which can be found in the coding or the untranslated regions (UTR) of the mRNA, mostly depends on a part called the seed sequence (nucleotides 2–7 at the

MicroRNAs (miRNAs) are a class of small, non-coding RNAs that are recognized as critical regulators of immune gene expression during infection. Many immunologically significant human miRNAs have been found to be conserved in agriculturally important species, including cattle. Discovering how bovine miRNAs mediate the immune defense during infection is critical to understanding the etiology of the most prevalent bovine diseases. Here, we review current knowledge of miRNAs in the bovine genome, and discuss the advances in understanding of miRNAs as regulators of immune cell function, and bovine immune response activation, regulation, and resolution. Finally, we consider the future perspectives on miRNAs in bovine viral disease, their role as potential biomarkers and in therapy.

**Keywords: bovine, Bos taurus, microRNA, miRNA, immune system**

5 0 end) of the miRNA (7). Although both strands of the duplex can be functional, usually only one strand is used (8).

MicroRNAs are evolutionarily conserved and have been found in all eukaryotes from unicellular species to mammals (9). To date, 2588 have been identified in the human genome, 1915 in mouse, and 793 in bovine (miRBase version 21, http://www.mirbase.org) (10). Originally, miRNAs were thought to mainly regulate gene expression by inhibiting translation, however, transcriptional regulation has been shown recently to be the primary mechanism used by miRNAs to influence gene expression in mammals (11). Several important factors determine how miRNAs function, including the location and number of target sites within mRNA 3<sup>0</sup> UTRs (or other gene regions), RNA-binding proteins (RBPs), which interfere with RISC binding, RISC co-factors, and the modification of RISC-components (12). Transcripts with target sites for a common miRNA compete for recognition. These competing endogenous RNAs (ceRNA), such as other mRNAs, long non-coding RNAs (lncRNA), pseudogene transcripts, and independently transcribed UTRs, can reduce the effect of specific miRNAs *in vivo*. Binding to target sites also protects miRNAs from degradation, in a process called target mediated miRNA protection (TMMP) (12). Individual miRNAs can have many targets, explaining how relatively small numbers of alternatively expressed miRNAs can have a large impact on gene regulation and control several diverse biological processes.

#### **MicroRNAs IN IMMUNITY AND INFECTION**

Cell-type specific expression of miRNAs has been demonstrated in innate and adaptive immune cells (13, 14) and there is a growing body of evidence that miRNAs regulate the differentiation, development, and function of these cells (15–17). Hematopoietic stem cell differentiation into myeloid and lymphoid lineages, for example, has been shown to be under the influence of several miRNAs, including miR-125b, miR-126, and miR-196b (4, 18). Additionally, the deletion of the *dicer1* gene, which is critical for proper miRNA processing, results in impaired T cell development (19), while miR-17–92, miR-150, and miR-155 have been demonstrated to be critical for B cell development. Other roles for miRNAs in regulating adaptive immunity have also been shown, including the regulation of B and T lymphocyte functions, including antibody production, by miR-155 (20–22). Activation of the innate immune system is also regulated by miRNAs (23). The human miRNA, miR-146a, for example, has been shown to target tumor necrosis factor receptor-associated factor 6 (TRAF6) and interleukin-1 (IL1) receptor-associated kinase (IRAK1), key regulatory nodes, which control innate immune signaling in response to lipopolysaccharide (LPS) (24). Similarly, miR-19a has been shown to regulate expression of SOCS 3, an important suppressor of cytokine signaling (25).

MicroRNAs have also been clearly demonstrated to have important roles in regulating responses to infection (26). In particular, several miRNAs have been identified to have important functions in regulating immune responses to mycobacterial infection (27). Tumor necrosis factor (TNF) biosynthesis, for example, is inhibited by *Mycobacterium tuberculosis* – an intracellular mycobacterial pathogen that infects alveolar macrophages – by altering levels of human macrophage miRNAs, including miR-125b and miR-155, for its own benefit (28). Similarly, miR-29 and miR-99b regulate the production of multiple cytokines, including IFN-γ and TNF-α, which control *M. tuberculosis* growth (29, 30). miRNAs are frequently evolutionarily conserved and many of these miRNAs have orthologs in cattle, therefore data from human and mouse studies can provide a roadmap for revealing miRNAs likely to have important roles in bovine infectious diseases. Many miRNAs, however, exhibit pathogen or stimulus-specific response profiles and certain families of miRNAs are expanded or contracted in the bovine genome.

#### **MicroRNAs IN THE BOVINE GENOME**

The first studies demonstrating miRNA expression in bovine tissues were undertaken in 2007 (31, 32). Since then, 793 mature miRNAs, encoded on all 30 chromosomes, have been identified in the *Bos taurus* genome. These miRNAs account for approximately a quarter of all the 3825 non-coding RNAs predicted in the genome by Ensembl (v75) (33). Typically, miRNAs have been grouped into families based on shared sequence similarity of the miRNA seed region (2–8 nt), the mature sequence, or the precursor miRNA sequence (34). Often, miRNA families can be found clustered with target genes in specific genomic regions (35). Many human miRNAs, including some of the most extensively studied immune-related miRNAs, share significant functional and sequence similarities to their bovine counterparts indicating evolutionary conservation and, putatively, conservation of function. The human miRNA, hsa-miR-155, for example, is a perfect homolog to its bovine counterpart bta-miR-155. In humans, this miRNA acts as an anti-inflammatory agent targeting

the Toll-like receptor/Interleukin-1 receptor (TLR/IL1R) inflammatory pathway (36). Another miRNA with a conserved bovine ortholog, hsa-miR-146a-5p, is known to negatively regulate the retinoic acid-inducible gene 1 (RIG-I) pathway in humans by suppressing TRAF6 and IRAK1 during viral infection (37). There is also an exact seed sequence match between hsa-miR-146a-5p, bta-miR-146a, and mmu-miR-146a-5p.

While there is significant conservation of miRNAs between species, there are also notable differences that very likely have functional consequences. There are numerous cases, for example, of miRNAs found in the human genome that are apparently absent in bovine. Some of these differences may be due to better annotation of the human microRNAome but clearly there are real differences too. The human miRNA, hsa-miR-198, for example, has a role in human immunity and has no apparent homolog in the bovine genome. This miRNA targets the Cyclin T1 gene (CCNT1), which acts as a co-factor for HIV-1 (38).

In addition to single miRNA differences in the repertoire of human and bovine miRNAs, there are also several cases where entire families or clusters of miRNAs that are present in human have yet to be discovered or do not exist in the bovine genome. These include the majority of miRNAs numbered from miR-550 to miR-640; some 200 miRNAs, which include the hsa-miR-515 cluster (11 miRNAs), and interestingly, the miR-548 family. The miR-548 family comprises of over 70 miRNAs whose expression to date has only been described in simians. Members of this miRNA family have been shown to target interferon-λ1 (IFN-λ1), modulating the primate interferon response to viral infections (39).

There are also several miRNA families in the bovine genome that are apparently species-specific, at least when compared to available genomes. The bta-miR-2284 and bta-miR-2285 families, for example, encode more than 100 mature miRNAs in the bovine genome but do not appear to have homologs in either human or mouse. These miRNA families have been shown to be expressed in a number of bovine immune-relevant tissues including CD14<sup>+</sup> monocytes, mammary epithelial cells, and alveolar macrophages (40–42), however, the genes targeted by the miRNAs in this family are currently unknown.

#### **ROLE OF THE miRNAs IN THE BOVINE IMMUNE SYSTEM**

The roles that miRNAs play in regulating immune activation and resolution in response to infection is less well understood in cattle, compared to human and mouse. Investigations in cattle have primarily focused on characterizing miRNA expression in immune-related tissues and investigating whether miRNA expression is perturbed during bacterial/viral infections – but detailed mechanistic studies are, to date, largely lacking. One of the first studies to profile immune-relevant miRNA expression in cattle, characterized the expression of more than 100 bovine orthologs of known human miRNAs, as well as novel bovine miRNAs, in several immune-relevant tissues and provided an early bovine miRNA expression atlas for many later studies (31). More recently, Vegh et al. utilized a next generation sequencing (NGS) strategy to profile miRNA expression on a genome-wide scale in bovine alveolar macrophages, the primary host cell for *M. bovis*, an economically important pathogen (41). NGS has several advantages over classical sequencing technologies, which include the

ability to accurately measure expression of all miRNAs simultaneously, with high reliability, single-nucleotide resolution and across the broad dynamic range of expression (43). miRNA expression has now been demonstrated in 10 immune-related bovine tissues (bovine embryo, thymus, small intestine, mesenteric lymph node, abomasum lymph node, CD14<sup>+</sup> blood isolated monocytes, CD14<sup>+</sup> milk-isolated monocytes (MIMs), mammary epithelial cells, alveolar macrophages, mammary tissue, and in the MDBK cell line) and tissue-specific expression of miRNAs in these tissues is readily apparent (**Figure 1**) (31, 40–42, 44–47).

Several studies have also examined whether miRNA expression is altered in response to Gram-positive or Gram-negative infections associated with bovine mastitis, a disease of the bovine mammary gland with significant economic consequences to the dairy industry. Naeem et al. examined a panel of 14 miR-NAs from mammary tissue biopsied during an *in vivo* intramammary *Streptococcus uberis* infection (46). Four of the fourteen miRNAs were found to undergo differential expression (**Table 1**). Another study compared transcriptional changes of five inflammation-associated miRNAs, including ones with extensively studied human orthologs: bta-miR-155, bta-miR-146a, and btamiR-223, in bovine CD14<sup>+</sup> monocytes stimulated with either LPS or *Staphylococcus aureus* enterotoxin B (SEB) (45). Four miRNAs were differentially expressed, and LPS was found to have a greater effect than SEB at inducing miRNA differential expression.

More recent studies have employed high-throughput sequencing approaches to temporally profile genome-wide changes in miRNA expression in different cell-types in response to challenge with bovine mastitis-causing pathogens such as *Escherichia coli*, *S. aureus*, and *S. uberis*. Lawless et al. identified 21 miRNAs that were differentially expressed in bovine mammary epithelial (BMEs) cells challenged *in vitro* with live *S. uberis*, a Gram-positive bacterium (40). Strikingly, the 21 miRNAs differentially expressed in response to this Gram-positive bacterium was substantially different to LPS-responsive miRNAs. Furthermore, the predicted target genes of miRNAs that were down-regulated in BMEs following *S. uberis* infection (but not the targets of up-regulated miRNAs), were statistically enriched for roles in innate immunity, suggesting that the repression of these miRNAs transcriptionally releases the innate immune response to this infection. Subsequently, Jin et al. also reported notable differences in miRNA expression profiles in MAC-T cells challenged with either heat-killed *E. coli* or *S. aureus* (47). This was the first bovine study to directly

**FIGURE 1 | Multidimensional scaling (MDS) plot of miRNA profiles shows that miRNA expression is cell-type specific**. Similar cell-types have similar profiles, as evidenced by milk (+) and blood CD14<sup>+</sup> cells (o), and by bovine mammary epithelial cells [primary (\*) and cell line (M)], which group together. Data are from four recent RNA-seq publications. +, Blood CD14<sup>+</sup> monocytes; o, milk CD14<sup>+</sup> monocytes; ×, bovine alveolar macrophage; \*, primary bovine mammary epithelial cells (BME); M, bovine mammary epithelial cell line (MAC-T).


#### **Table 1 | Summary of bovine miRNA literature**.

compare global miRNA expression of two pathogens in the same cell-type.

In the first NGS-based study to temporally profile infectioninduced miRNA responses *in vivo*, Lawless et al. simultaneously profiled genome-wide mRNA and miRNA expression at multiple time-points in both milk and blood FACS-isolated CD14<sup>+</sup> monocytes from cattle infected with *S. uberis* (42). Twenty-six miRNAs and more than 3500 genes were identified as being significantly differentially expressed over the 48 h time-course. In MIMs, up-regulated protein-coding genes were significantly enriched for inflammatory and innate immune pathways, while down-regulated genes were enriched for non-glycolytic metabolic pathways. Monocyte transcriptional changes in the blood were more subtle. Pathway analysis revealed that predicted targets of MIM down-regulated miRNAs were highly enriched for roles in innate immunity, while up-regulated miRNAs preferentially targeted genes involved in metabolism; suggesting that during *S. uberis* infection miRNAs are key amplifiers of monocyte inflammatory response networks and repressors of several metabolic pathways.

To date, only one study has examined whether bovine miR-NAs are transcriptionally perturbed during viral infection. In this study, an adult bovine kidney cell line was challenged with *bovine herpes virus 1*. Sequencing 3 miRNA libraries, more than 300 miR-NAs were found to expressed, but none were described as being differentially expressed (44).

#### **VALIDATED BOVINE miRNA TARGET GENES**

Two approaches are commonly used to validate miRNA targets. One involves transfection of miRNA mimics (or inhibitors) into cells and confirmation that the expression of predicted target genes is altered as expected (50). This approach provides evidence that a given miRNA can alter the expression of a target gene but does not prove direct regulation (51). The miRNA could, for example, regulate a transcription factor, which subsequently regulates the putative target gene. A more direct approach is to use a reporter assay, where the 3<sup>0</sup> UTR of the predicted target gene (or at least the portion containing the predicted miRNA binding site) is cloned upstream of a luciferase or green fluorescent protein (GFP) reporter gene. If binding of the transfected miRNA mimic to the 3 <sup>0</sup> UTR reduces the level of reporter protein this demonstrates a direct silencing effect of the miRNA on the gene (52).

To date, very few studies have functionally validated bovine miRNA targets. One example where miRNA regulation has been validated is bovine High Mobility Group Box 1 (HMGB1), a nuclear protein, which transcriptionally regulates inflammation (48). Bovine HMGB1 has been shown to be targeted by bta-miR-223, a miRNA that is up-regulated during *S. aureus* infection (48). Similarly, using reporter assays, bta-miR-124 has been shown to regulate expression of Monocyte Chemotactic Protein 1 (MCP1) and Polypyrimidine Tract Binding Protein 1 (PTBP1) in bovine fibroblasts (53).

#### **OTHER AREAS OF BOVINE miRNA IMMUNE BIOLOGY**

In addition to their role as intracellular transcriptional modulators of gene expression miRNAs are also stably expressed in a host of extracellular body fluids including milk, saliva, semen, and plasma (54–57). Extracellular miRNAs can be transferred to distant recipient cells via exosome-mediated transfer and have been demonstrated, in mouse dendritic cells, to modulate recipient cell transcription (58). Exosome packaged miRNAs have been shown to be highly stable and are resistant to degradation by RNases, freeze-thaw, and low pH (56, 59). Exosome miRNAs – including a number of immune-relevant ones, such as bta-miR-223 and bta-miR-125b – have been found in both human breast milk and bovine milk (54–56). MicroRNA expression levels in milk have also been observed to vary during different lactation periods and are present in milk products as well as in raw milk (54). Interestingly, they are particularly abundant in colostrum. Further investigation of the role of exosome packaged miRNAs play in regulating mammalian immunity is urgently needed.

#### **FUTURE RESEARCH DIRECTIONS IN BOVINE miRNAs AND THEIR EFFECT ON IMMUNOLOGY, INFLAMMATION, AND INFECTION**

It is clear that miRNAs play a key role in regulating human and mouse immune responses. In cattle, studies to date have been mainly limited to demonstrating differential expression of miR-NAs in immune-relevant tissues or cells challenged *in vitro* with specific pathogens. Importantly, annotated miRNAs are much fewer in the bovine genome than in murine or human genomes, and a bovine miRNA expression atlas across bovine tissues and cells is needed to bridge this gap. Among its many uses, better annotation of non-coding RNAs would aid in the interpretation of bovine genome-wide association study (GWAS) data. A previously unannotated small non-coding RNA, for example, was recently identified as the only gene in a novel genome-wide significant QTL for somatic cell score, a mastitis indicator trait (60).

Aside from understanding the important basic biology of how miRNAs regulate bovine gene expression, miRNAs could also potentially be of significant utility as biomarkers of specific diseases in cattle. Indeed, miRNAs exhibit many properties that have made them of significant interest as non-invasive biomarkers in human clinical studies. miRNAs are abundantly expressed, in a stable form, in a range of extracellular fluids, are easily measured, and in many cases exhibit temporal and spatial specificity (54– 57). They also have high information content – small numbers of miRNAs can serve as accurate biomarkers. Several studies have investigated miRNA expression profiles associated with different mastitis-causing pathogens (**Table 1**), however, many of these studies were carried out *in vitro*, and the potential of miRNAs as biomarkers of bovine disease is currently limited by a lack of studies of *in vivo* comparison of miRNA profiles associated with multiple different pathogens in the same challenge model. Such *in vivo* studies are needed to identify sensitive and specific biomarkers for particular infections. This could be used for identifying infections, such as tuberculosis, using a simple miRNA-based biomarker, or for distinguishing between different infections, for example between *E. coli* and *S. uberis* driven mastitis, helping veterinarians to select more specific therapeutic strategies.

A further limitation is the fact that research investigating the role of miRNAs in regulating bovine immunity has, to date, focused almost exclusively on bacterial infections and little is known about the role miRNAs play during bovine viral infections.

Some of the most economically important and high profile bovine infectious diseases are of viral origin including Foot and Mouth Disease Virus (FMDV) (61), Bovine Viral Diarrheal Virus (BVDV) (62), and the recently described Schmallenberg virus (63). In other species, it is clear that host miRNAs have a direct role in modulating the host immune response to pathogenic viral infections (38). Additionally, some viruses encode their own repertoire of miRNAs to subvert the host immune response. Nearly 200 viral miRNAs have been described (64). Describing their precise effect on the immune response could further our understanding of both bovine miRNA immune biology and virology.

Aside from their utility as biomarkers, miRNAs also have significant potential as therapeutic targets or agents. MicroRNA function can be augmented either by over-expression approaches, using miRNA mimics or vector based over-expression, or by inhibition,using miRNA sponges or anti-miR oligonucleotides (65). In humans, several miRNAs are currently in preclinical and clinical trials as novel therapeutics in cancer, viral infections, and cardiovascular disease (65). Human miR-122, for example, is being investigated for its therapeutic potential to modulate cholesterol metabolism (66). Additionally, targeting miR-122 using the antimiR miravirsen induces antiviral activity against hepatitis C virus (HCV) (67). The potential clinical utility of miR-122 is being investigated in Phase II clinical trials. Similarly, human miR-208 has been shown to have an important role in modulating cardiac function and remodeling (68), and is currently in preclinical trials. Interestingly from a bovine perspective, this miRNA also has a big impact on metabolism. Treatment with anti-miR-208 prevented weight gain in aging mice, which was due to a reduction in fat weight (69).

All of the miRNAs mentioned above with therapeutic potential have orthologs in cattle and these examples clearly suggest that there is a potential for the application of miRNA-based therapeutic strategies to combat disease and regulate metabolism in cattle in order to influence important economic traits, such as growth, feed efficiency, or milk production. Although the cost of miRNA therapy and the large size of animals may prevent agricultural use, bovine research models still could be valuable, as the conserved nature of miRNAs facilitates translation of research to human application (65).

A current limitation to the translational potential of miRNA biology in cattle is the lack of validated targets for known miRNAs, as only a handful of studies to date have functionally validated predicted miRNA targets (40, 41). Computationally, hundreds if not thousands of putative miRNA targets can be predicted and experimental validation is a costly and labor-intensive procedure. Methods which integrate and correlate miRNA expression with mRNA expression in the same sample can refine computational predictions and increase the validation hit rate. Other more recent technological advancements such as crosslinking immuneprecipitation sequencing (CLIP-Seq) can directly identify miRNA targets on a genome-wide scale (70), but such approaches have yet to be implemented in cattle.

In conclusion, miRNAs undoubtedly play a key role in regulating bovine immunity and disease. Future studies are poised to reveal their true potential as novel biomarkers or therapeutic agents in a range of bovine diseases as well as providing further insight into the fundamental biology of how they regulate bovine immune gene expression, insight which is essential before their translational potential can be realized fully.

#### **ACKNOWLEDGMENTS**

Nathan Lawless and Peter Vegh are supported by the Teagasc Walsh Fellowship scheme. David Lynn is supported by EMBL Australia.

#### **REFERENCES**


**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.

*Received: 05 September 2014; accepted: 13 November 2014; published online: 27 November 2014.*

*Citation: Lawless N, Vegh P, O'Farrelly C and Lynn DJ (2014) The role of microRNAs in bovine infection and immunity. Front. Immunol. 5:611. doi: 10.3389/fimmu.2014.00611*

*This article was submitted to Molecular Innate Immunity, a section of the journal Frontiers in Immunology.*

*Copyright © 2014 Lawless, Vegh, O'Farrelly and Lynn. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

## The contribution of the maternal immune system to the establishment of pregnancy in cattle

#### **Trudee Fair \***

School of Agriculture and Food Sciences, University College Dublin, Dublin, Ireland

#### **Edited by:**

Kieran G. Meade, Teagasc, Ireland

#### **Reviewed by:**

Annapurna Nayak, Brunel University, UK Peter James Hansen, University of Florida, USA Robert Owen Gilbert, Cornell University, USA

#### **\*Correspondence:**

Trudee Fair, School of Agriculture and Food Sciences, University College Dublin, Room 239, Agricultural Sciences Building, Belfield, Dublin 4, Ireland e-mail: trudee.fair@ucd.ie

Immune cells play an integral role in affecting successful reproductive function. Indeed, disturbed or aberrant immune function has been identified as primary mechanisms behind infertility. In contrast to the extensive body of literature that exists for human and mouse, studies detailing the immunological interaction between the embryo and the maternal endometrium are quite few in cattle. Nevertheless, by reviewing the existing studies and extrapolating from sheep, pig, mouse, and human data, we can draw a reasonably comprehensive picture. Key contributions of immune cell populations include granulocyte involvement in follicle differentiation and gamete transfer, monocyte invasion of the periovulatory follicle and their subsequent role in corpus luteum formation and the pivotal roles of maternal macrophage and dendritic cells in key steps of the establishment of pregnancy, particularly, the maternal immune response to the embryo.These contributions are reviewed in detail below and key findings are discussed.

**Keywords: cow, fertility, macrophage, cytokine, immune function**

#### **BACKGROUND**

It is estimated that fetal viability is only achieved in about 55% of fertilizations in non-compromised cattle, indicating an embryonic/fetal mortality of about 35%. It is estimated that 70– 80% of the total embryonic loss occurs between days 8 and 16 after insemination [day 16 corresponding to the day of maternal recognition of pregnancy; reviewed by Diskin and Morris (1)]. There are many reasons, related to both the mother and the embryo, why implantation fails but there is increasing interest in the role of the maternal immune system. Disturbed or aberrant immune function has been identified as primary mechanisms behind infertility. In contrast to the extensive body of literature that exists for human and mouse, studies detailing the immunological interaction between the embryo and the maternal endometrium in cattle have primarily focused on the role of the maternal recognition factor, the type I antiviral cytokine, interferon tau (2, 3) in corpus luteum (CL) maintenance, and progesterone priming of the endometrium. Nevertheless, by reviewing the existing studies and extrapolating from sheep, pig, mouse, and human data, we can draw a reasonably comprehensive picture of immune cell involvement from follicle development, ovulation, gamete transfer, maternal recognition of pregnancy, implantation, and placentation. These events are reviewed below and key findings are discussed.

#### **OVARIAN FUNCTION**

The presence and temporal regulation of neutrophils, eosinophils, macrophages (MΦ), granulocytes, and T-lymphocytes in ovarian tissues has been characterized extensively during the menstrual cycle in women; a smaller body of data exists for several farm animal species, including cows, sheep, pigs, buffaloes, and horses (4).

#### **FOLLICLE DIFFERENTIATION**

Taken together, pre-ovulatory follicle differentiation and luteinization appear to be characterized by three phases of immune cell infiltration, which are illustrated in **Figure 1**: histological analysis of bovine dominant follicles shows that mast cell infiltration of the theca layer constitutes the first phase (5) (**Figure 1A**), luteinizing hormone (LH) triggered degranulation of the mast cells stimulates the second phase through the direct and indirect actions of TNF-alpha (TNFA), a constituent of the granules (**Figure 1B**). The second phase has been characterized in sheep and pigs as an influx of eosinophilic and neutrophilic granulocytes and Tlymphocytes (6, 7). The last phase of leukocyte migration consists of phagocytic monocytes (Mo); MΦ's increase in the sow and ewe follicles at the time of ovulation (6, 7) (**Figure 1C**), possibly in response to peak estradiol concentrations (8). The temporal changes in the influx of leukocytes appear to occur in response to various chemoattractant cues produced by the developing follicle (9), indeed leukocyte chemoattractant activity has been demonstrated in bovine, ovine, and human follicular fluid of ovulatory follicles (10–12). Immunohistochemical characterization of the immune cell repertoire of the bovine ovary has largely focused on the formation and regression of the CL, which will be discussed later. In contrast, there are many reports detailing the transcriptomic profile of ovarian follicle development in cattle: (13–19). In particular, the deep sequencing analysis of bovine follicular theca and granulosa tissue during pre-ovulatory follicle development, revealed a profound effect of ovarian follicle stage on the expression of many genes within immune-related pathways in these tissues: during follicle differentiation, bovine thecal tissue was characterized by the expression of immune factors associated with vascularization, angiogenesis, and cellular proliferation (15, 20), processes which are carried out by MΦ's in the theca layer during this time (21, 22). The bovine transcriptomic data

also concurred with the histological findings described for sheep and pigs, as factors with known inflammatory/chemotactic properties such as *AKT2*, *ARHGEF1*, *GNAI2*, *IL-1*, *IL-6*, and *IL-8b* (23–25) were upregulated and pathways associated with MΦ and neutrophil function were overpopulated in differentiating thecal tissue (15).

#### **FOLLICLE LUTEINIZATION AND OVULATION**

Findings from studies using rodent models indicate that the initiation of the ovulatory process occurs primarily in granulosa cells (26). Following the pre-ovulatory LH surge, morphological, endocrinological, and biochemical changes occur in the theca and granulosa cells, which redirect pre-ovulatory follicle development from differentiation to luteinization and thus the early stages of CL development (26). In particular, the post-LH deterioration of the basement membrane (BM) between the theca and granulosa-cell layers (GCs) (27), facilitates the movement of leukocytes into the granulosa tissue at luteinization, reflected by the peri-ovulatory granulosa-cell expression of factors involved in acute inflammation and immunosurveillance (15,26,28). It has been hypothesized that the dramatic increase in the expression of these signals in the follicle compartment activates the ovarian innate immune system (29) and that the damaged granulosa cells actively secrete alarmins or passively release them after death (30). Alarmins include acute

phase proteins (APP), S100 proteins, advanced glycation endproducts (AGE), high mobility group box-1 protein (HMGB1), defensins, and interleukin (IL)-1α, which are all present in follicle cells and the follicular fluid of pre-ovulatory follicles (26, 31–34) and can engage toll-like receptors (TLRs). In the case of ovarian granulosa cells, oxidized low density lipoprotein (oxLDL), which engages with TLR4 (34), has been proposed as a key alarmin in the pre-ovulatory cascade (29). This hypothesis is further supported by the identification of granulosa-cell exclusive expression of TLR signaling and NF-κB signaling pathways during luteinization in the bovine transcriptomic data (15). Furthermore, comparing the gene expression profiles of follicular tissue from heifers to that of lactating cows, it would appear that the recruitment of leukocytes to the differentiating follicle is delayed in cows. This is possibly a result of the demands of parturition/lactation in dairy cows, resulting in a reduced positive feedback loop, whereby lower steroid levels and chemoattractant signals recruit fewer leukocytes into the follicle, leading to lower steroid and chemoattractant levels (15).

#### **CORPUS LUTEUM FORMATION**

The CL is a transient organ established by cells of the follicle following ovulation; it is composed of a heterogeneous mixture of cell types that consist of not only steroidogenic luteal cells but also non-steroidogenic cells including vascular endothelial cells, fibroblasts, and immune cells such as lymphocytes and MΦ's (35). Studies in human, rat and sheep indicate that the immune cells of the developing CL are recruited during ovulation (36– 39) (**Figure 1D**), they were determined to have originated from the spleen (38), see Ref. (40), for review. Histological data from cattle indicate that they are primarily granulocytes, neutrophils, and eosinophils (35, 41, 42). However, as CL development and vascularization progresses, MΦ's and endothelial cells infiltrate (43), providing a source of TNF and TNFR, the presence of which have been demonstrated in the bovine CL (44). TNF is a potent stimulator of luteal prostaglandins (PG) including PGF2a, PGE2, and PG12 (45), TNF and TNF-induced PGE2 have been proposed as key regulators of CL vascularization (46), recent work in the mare supports this hypothesis (47). Exposure to seminal plasma has been shown to enhance CL development and ovarian steroidogenesis: gilts treated with seminal plasma had heavier CLs, higher plasma progesterone (P4) levels, which peaked earlier, without a concurrent increase in ovulation rate, suggesting that the number and output of steroidogenic luteal cells is greater in animals exposed to seminal components (48). Immunohistochemical analysis revealed a greater abundance of predominantly major histocompatibility complex (MHC) class II positive MΦ's and/or DCs in the stromal tissues and thecal cells of pre- and periovulatory follicles, implying greater leukocyte recruitment at the time of ovulation in seminal plasma treated animals (49).

Immune function is central to CL regression, which must occur in the absence of pregnancy in order for new follicular development to take place (40). The regressing CL is characterized by an increase in MΦ and Mo populations, which eventually constitute the major proliferating cell type of the late regressing CL (40). The number of T-lymphocytes appears to increase just prior to the onset of luteolysis (35, 50), analysis of the bovine CL T-lymphocyte population revealed that 25% of T-lymphocytes present in a functional CL were T helper cells (CD4+), 45% were cytotoxic T-cells (CD8+), and 30% were gamma delta (γδ+) T-cells and that this profile did not alter during luteolysis (51). However, decreased P4 levels and interruption of growth factor signaling in the CL appear to promote both MΦ and T-cell activation, leading to increased TNF and INF production, respectively (52–55). TNF and INF are likely to be key regulators of apoptosis and ovarian tissue remodeling (56), their receptors are expressed in bovine steroidogenic cells and luteal cells (57). It is probable that Fas expression is induced in luteal cells by leukocyte-derived cytokines and that Fas L expressed on T-lymphocytes transduces apoptotic signals to the luteal cells [see Ref. (46), for review]. This is likely to be a conserved action as both Fas and Fas L are expressed in theca cells in multiple species (58).

#### **GAMETE TRANSFER**

#### **INFLAMMATORY RESPONSE TO INSEMINATION**

The site of semen deposition is very much a species-specific location (59). In cattle, and also in humans, sperm enters the cervix canal rapidly after semen deposition. The stimulation of vaginal insemination ensures the migration of neutrophils in to the cervical and uterine tissues (60, 61) and has been proposed as the initial point to optimize pregnancy success (62). The early immune response to insemination appears to contribute to both the ovulatory process and sperm cell selection; as reports from several species, including cattle, suggest that neutrophilic granulocytes target dead or capacitated sperm, thus removing non-motile or damaged spermatozoa (63–65),rather than motile,fertile sperm (62). In both humans and mice, it has been clearly demonstrated that the post-mating inflammatory response is mainly caused by the seminal plasma, with sperm having a negligible part (66). The cytokine, transforming growth factor-β (TGFβ), is the principal inflammatory trigger found in seminal plasma; it is primarily present within the male seminal plasma fluid in latent form, which is activated in the female reproductive tract by plasmin and other enzymes after insemination (62, 67). Although TGFβ itself can be chemotactic for a variety of immune cell types (68), in the murine uterus it was reported to act indirectly, by inducing cytokine and chemokine expression (69).

#### **SPERM TRANSPORT**

The delivery of seminal fluid to the female reproductive tract at coitus represents the first exposure of the female immune system to paternal alloantigens (62), raising the possibility that the female activates an immune response to male antigens in seminal fluid that may ultimately confer immunological tolerance to paternal antigens (70). This theory is supported by data from mice, which show that chemoattractants, secreted by eosinophils and neutrophils, attract both Mos and DC's and shape the inflammatory status of MΦ's (71, 72). The response is not restricted to vaginal exposure; intrauterine horn insemination was shown to induce recruitment of MHC class II positive cells in gilts (73). Seminal plasma contains estrogen and testosterone, PG, and various signaling molecules, including IL-8, TGFβ, and IFNG, as well as bacterial lipopolysaccharide (LPS) (62). When murine uterine and cervical cells come into contact with the constituents of semen, they are stimulated to synthesize and release granulocytemacrophage colony-stimulating factor (GM-CSF), IL-6, and further chemokines (66, 74), which stimulate MΦ's, DC, and granulocyte infiltration of the uterine and cervical tissues (75). The induction of IL-6 is required for TGFβ to induce the generation of IL-17 producing, pro-inflammatory TH-17 cells, which in turn favor the induction of neutrophil-chemotactic IL-8 (76). In conjunction with IL-8, TGFβ induces the secretion of pro-inflammatory cytokines such as IL-1B, IL-6, and leukemia inhibitory factor (LIF) (77). Although, the expression of TGFβ was shown to increase in the bovine endometrium during the implantation period (78), the relatively high pregnancy rates achieved in cattle following artificial insemination or embryo transfer undermines the importance of maternal exposure to seminal plasma in cattle. The findings of studies designed to address this point indicate that neither exposure to seminal plasma nor TGFβ is critical for to the establishment of pregnancy in cattle (79, 80).

#### **IMMUNE TOLERANCE POST-FERTILIZATION**

Exposure to paternal antigens occurs in two waves in the reproductive process: initially during transmission of seminal fluid at coitus, and secondly when placental trophoblast cells come in contact with maternal tissues during embryo implantation (81). In sheep and cattle, morula-stage embryos enter the uterus around day 4–5 and blastocysts are formed by day 6 and 7, respectively, hatching occurs within 48–72 h. The hatched blastocyst subsequently elongates reaching 10 cm or more in length by day 14 and 25 cm or more in length by day 17 and the conceptus trophectoderm and endometrial luminal epithelium (LE) become closely apposed,see Ref. (82),for review. Implantation is a superficial,protracted affair in these species, commencing after attachment and adhesion of the trophectoderm to caruncular and intercaruncular areas on day 16 in sheep and day 19 in cattle. Again, in contrast to the volume of data that has been acquired in human and mouse studies, the number of investigations carried out in farm animal species on the involvement of the maternal immune system in the establishment of pregnancy is very limited, particularly, for early pregnancy. For several decades, human pregnancy was described as a Th1/Th2 dichotomy with an imbalance toward a Th2 type immune response (83, 84). However, this paradigm is considered a simplistic explanation of the molecular events occurring during pregnancy, as it does not account for reported endometrial expression of Th1-type cytokines during implantation (85, 86). In ruminants, studies investigating maternal immunomodulation by pregnancy have focused on the actions of the type 1 interferon, IFNT,which is secreted by the elongating conceptus and is the main signaling factor in maternal detection/recognition of pregnancy (87, 88). Initial studies demonstrated that endometrial luminal epithelial cell estrogen receptor and oxytocin receptor expression was down regulated in response to IFNT (89, 90). Critically for the continuation of pregnancy in cattle, this binding eventually results in the attenuation of endometrial PGF2a secretion, allowing CL production of P4 to be maintained (90). In addition to its anti-luteolytic properties, IFNT appears to be the key regulator of the maternal immune response in ruminants (91, 92), acting on the endometrium to induce or enhance the expression of genes hypothesized to regulate uterine receptivity to implantation and conceptus development (78, 93–95). The expression of IFNT is limited to the embryonic trophectoderm during the peri-implantation period (96). Additionally, there is significant evidence that the bovine conceptus does not endeavor to conceal itself immunologically, as MHC-I transcripts have been detected in early cleavage stage bovine embryos (97) and in first and second trimester and term trophoblast tissues (98). Furthermore, *MHC* class I mRNA expression by bovine embryos is both transcriptand embryo stage-specific (97) and can be regulated by a number of cytokines including IFNG, IL-4, and LIF (99, 100).

#### **MATERNAL RECOGNITION AND RESPONSE TO PREGNANCY MONOCYTES, MACROPHAGES, AND DENDRITIC CELLS**

Macrophage recruitment to the pregnant endometrium occurs in a wide range of mammalian species, including the mouse (101), human (102, 103), cynomolgus and vervet monkeys (104), sheep (105), and cattle (78, 106, 107). While their role has not been completely elucidated, functions include clearing of apoptotic cells, regulation of apoptosis (108), and regulation of placental lactogen concentrations at the fetal–maternal interface (109). Given the potential antigenicity of the conceptus due to paternal antigen and classical MHC protein expression (97), MΦ's may also feature in curtailing the activation of anti-conceptus immune responses (106). In cattle, the maternal immune response to the developing embryo is characterized by the expansion of Mo, MΦ's

(CD14+-cells), and DC (CD172a–CD11c+) populations in the endometrial stroma as early as day 13 of pregnancy (78). Interestingly, there was a parallel decrease in CD11b+-cells; CD11b is associated with Mo movement through the endothelium, which would imply that the Mo had acquired a stationary phenotype (78).

Dendritic cells have been shown to play an important role in decidua formation and the induction of immune tolerance in human and murine pregnancy (110, 111). Employing individual and combined CD172a and CD11c labeling of the bovine endometrium, it was determined that there was a high prevalence of immature cells within the endometrial DC population during early pregnancy (78). Immature DC's have been associated with the initiation and maintenance of peripheral tolerance (112) and their presence in large numbers in the uterine decidua has been associated with the establishment of healthy pregnancies in women (113). It is most likely that in cattle, IFNT induces this initial maternal response to the presence of the elongating embryo, either by attracting monocytes into endometrium or by modulating their differentiation into MΦ's or DC. Indeed, gene expression analysis of the same endometrial tissue revealed dramatic up-regulation of mRNA expression of IFN stimulated genes *IL12B*, *MCP1*, *MCP2*, *PTX3*, *RSAD2*, *ISG15*, and *TNFA* (78). Furthermore, MCP1 and MCP2 are members of the cellular chemoattractant chemokine β subfamily, which have highly potent MΦ recruitment and activation properties (114), thus increased *MCP1* and *MCP2* expression may be associated with the recruitment of Mo/MΦ from the systemic system into the endometrium. The up-regulation of the evolutionary conserved PTX3 is very interesting; gene deletion studies in mice have shown that it is essential for female fertility, participating in the assembly of the cumulus oophorus extra-cellular matrix (115). Moreover, PTX3 is involved in innate immunity, proposed roles include selected pathogen recognition, opsonization leading to enhanced phagocytosis, regulation of the inflammatory response, complement-mediated clearance of apoptotic cells, and control of autoimmunity (116–118).

#### **T-LYMPHOCYTES**

Despite the evidence from studies in humans and mice linking successful pregnancy with an imbalance toward a Th2 immune response type, data from cattle indicate that CD4+, CD8+, γδTCR+, and FoxP3 T-lymphocyte populations are not regulated temporally during estrus or early pregnancy in cattle (119). However, mRNA expression analysis on the same tissue revealed that the Th1 immune factors *IFNA*, *LIF*, *IL1B*, *IL8*, and *IL12A* were down regulated during the luteal phase of the estrus cycle, whereas the Th2 factors LIF and IL10 were upregulated, suggesting that the phenotypes/inflammatory status of Th cells are tightly modulated during the estrous cycle in anticipation of pregnancy. Additionally, LIF and IL-10 have been shown to regulate MΦ activation (120, 121). Similarly, endometrial *TGF*β*2* expression is down regulated during the ovine and bovine implantation period and is subsequently increased during placentation (78, 122), which may reflect TGFβ2 involvement in Mo recruitment and regulation of MΦ inflammatory status (123). Furthermore, the study in cattle showed TGFβ localization to the fetal–maternal interface of the bovine placentome, which may indicate TGFβ2 involvement in restricting trophoblast invasion

during the implantation phase, while enhanced expression during placentation and *in vitro* cell culture studies, suggest that TGFβ2 may play a mitogenic role during placentation, promoting caruncular growth, and coordinating epithelial cell development leading to placentome formation (123, 124).

Surprisingly, and in contrast to the situation in human and mouse models, where NK cells can constitute up to 70% of the endometrial lymphocyte population during the preimplantation phase of pregnancy (112), when uterine NK (uNK) cells are believed to play a pivotal role in local vascular remodeling and regulation of trophoblast invasion [for review see Ref.(125, 126)]; NK cells do not appear to play such a critical role during early pregnancy in cattle. Indeed, the only published data suggests the bovine endometrium population of CD335<sup>+</sup> NK cell population is not expanded as an immediate response to maternal recognition of pregnancy (119). The findings of an *in vitro* study which demonstrated anti-proliferative effects of recombinant IFNT exposure on immune and uterine cells, particularly leukocytes, infers that the IFNT secretion by the embryo may actively restrict NK cell expansion in early pregnancy (127), which is in keeping with the non-invasive nature of implantation in cattle [see review by Bazer et al. (128)]. However, further studies are required to determine if the NK cell population expands when IFNT secretion wanes and to what degree, if any, they are involved in placentation.

#### **CONCLUDING REMARKS**

Intensive cattle production systems have been associated with postpartum immunosuppression and subsequent reduced fertility; it is vital that basic research in the area of bovine reproductive immunology is expanded to generate new knowledge by which these issues can be overcome. However, although the number of studies investigating the contribution of the maternal immune system to reproductive function in cattle is a fraction of that carried out in human and mouse species, it is possible to conclude that maternal macrophage and dendritic cells play pivotal roles in key steps of the establishment of pregnancy, particularly, development of the CL and maternal immune response to the embryo.

#### **REFERENCES**


interval; rectal temperature; and uterine response to oxytocin in cyclic ewes. *Biol Reprod* (1997) **57**(3):621–9. doi:10.1095/biolreprod57.3.621


**Conflict of Interest Statement:** The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

*Received: 16 October 2014; paper pending published: 15 November 2014; accepted: 07 January 2015; published online: 28 January 2015.*

*Citation: Fair T (2015) The contribution of the maternal immune system to the establishment of pregnancy in cattle. Front. Immunol. 6:7. doi: 10.3389/fimmu.2015.00007 This article was submitted to Molecular Innate Immunity, a section of the journal Frontiers in Immunology.*

*Copyright © 2015 Fair. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

REVIEW ARTICLE published: 07 October 2014 doi: 10.3389/fimmu.2014.00493

## Bovine mastitis: frontiers in immunogenetics

### **Kathleen Thompson-Crispi 1,2, Heba Atalla1,2,3,4, Filippo Miglior 2,5 and Bonnie A. Mallard1,2\***

<sup>3</sup> Department of Biomedical Science, University of Guelph, Guelph, ON, Canada

<sup>4</sup> Department of Animal and Poultry Science, University of Guelph, Guelph, ON, Canada

<sup>5</sup> Canadian Dairy Network, Guelph, ON, Canada

#### **Edited by:**

Kieran G. Meade, Teagasc, Ireland

#### **Reviewed by:**

Miki Nakao, Kyushu University, Japan Hans-Martin Seyfert, Leibniz Institute for Farm Animal Biology, Germany

#### **\*Correspondence:**

Bonnie A. Mallard, Department of Pathobiology, University of Guelph, 50 Stone Road, Guelph, ON N1G 2W1, Canada e-mail: bmallard@uoguelph.ca

Mastitis is one of the most prevalent and costly diseases in the dairy industry with losses attributable to reduced milk production, discarded milk, early culling, veterinary services, and labor costs. Typically, mastitis is an inflammation of the mammary gland most often, but not limited to, bacterial infection, and is characterized by the movement of leukocytes and serum proteins from the blood to the site of infection. It contributes to compromised milk quality and the potential spread of antimicrobial resistance if antibiotic treatment is not astutely applied. Despite the implementation of management practises and genetic selection approaches, bovine mastitis control continues to be inadequate. However, some novel genetic strategies have recently been demonstrated to reduce mastitis incidence by taking advantage of a cow's natural ability to make appropriate immune responses against invading pathogens. Specifically, dairy cattle with enhanced and balanced immune responses have a lower occurrence of disease, including mastitis, and they can be identified and selected for using the high immune response (HIR) technology. Enhanced immune responsiveness is also associated with improved response to vaccination, increased milk, and colostrum quality. Since immunity is an important fitness trait, beneficial associations with longevity and reproduction are also often noted. This review highlights the genetic regulation of the bovine immune system and its vital contributions to disease resistance. Genetic selection approaches currently used in the dairy industry to reduce the incidence of disease are reviewed, including the HIR technology, genomics to improve disease resistance or immune response, as well as the Immunity+™ sire line. Improving the overall immune responsiveness of cattle is expected to provide superior disease resistance, increasing animal welfare and food quality while maintaining favorable production levels to feed a growing population.

**Keywords: disease resistance, genetic selection, genomics, immune response, mastitis**

#### **INTRODUCTION**

Mastitis, generally defined as the inflammation of the mammary gland, is a costly and complex disease associated with variable origin, severity, and outcome depending on the environment, pathogen, and host (1, 2). Mastitis is caused when pathogenic bacteria enter the sterile environment of the mammary gland, often as a result of disruption of physical barriers such as the teat, requiring prompt and appropriate host defenses to prevent colonization and subsequent disease pathology (3, 4). Mastitiscausing pathogens are commonly categorized as environmental or contagious, although this distinction has recently been disputed (5). Nonetheless, in general environmental pathogens have been grouped to include coliforms like *Klebsiella* or *Escherichia coli* (*E*. *coli*) and streptococci and are a major cause of clinical mastitis. On the other hand, those categorized as contagious pathogens can readily be spread from the infected quarters to other quarters of the same cow, or other cows and include *Staphylococcus aureus* (*S*. *aureus*) and *Streptococcus agalactiae* (6–8). Cow factors including age, stage of lactation, and somatic cell score (SCC) history

are known to influence the occurrence of mastitis infection (9, 10). The diverse pathogens that can cause mastitis induce different immune responses in the mammary gland, and therefore, the host requires highly specific pathogen-dependent responses for protection (11, 12).

Mastitis infections are described as either subclinical or clinical. Subclinical mastitis is the presence of infection without local inflammation resulting in an absence of visual signs (1). It may involve transient cases of inflammation and abnormal milk, and if this persists for longer than 2 months is termed chronic. Clinical mastitis, on the other hand, is an inflammatory response causing visibly abnormal milk. In the case of mild or moderate clinical mastitis, changes in the udder may include swelling, heat, pain, and redness. It is termed severe if the response includes systemic involvement such as fever, anorexia, and shock (13, 14). The diversity as well as the variation in prevalence and abundance of mastitis-causing organisms as well as the variation in host responses make mastitis a complex disease that continues to be a burden for the dairy industry.

<sup>1</sup> Department of Pathobiology, University of Guelph, Guelph, ON, Canada

<sup>2</sup> Center for Genetic Improvement of Livestock, University of Guelph, Guelph, ON, Canada

The bovine mammary gland is equipped with a non-immune anatomical barrier, and a plethora of immune-mediated defense mechanisms that include innate and adaptive immune responses. Innate immunity is relatively non-specific with rapid kinetics while the adaptive immunity offers a highly specific response with relatively delayed kinetics (15). Innate host-defenses depend on germline-encoded receptors that recognize conserved structures expressed by a wide range of microbes, and early induced cellular and soluble defenses. These natural defenses respond quickly to microbes during early stages of infection and are tightly integrated with the adaptive immune system. The innate host defenses of the mammary gland have been reviewed extensively elsewhere (16–18). The adaptive immune system uses a diverse repertoire of antigen specific receptors expressed by clonally expanded B and T-lymphocytes to regulate or eliminate the signal elicited by recognition events. Additionally, the induced adaptive immune response has the capacity to establish antigen specific memory for a rapid and augmented response upon subsequent exposure to the same antigen (19). For example, these various components of the immune system work in collaboration both locally and systemically in an attempt to control specific mastitis pathogens invading the mammary gland, but the details of the response is contingent upon the stage of infection and nature of the pathogen, as well as its interaction with the genetics of the host.

The interaction between mastitis pathogens and the host immune system is intricate, since both have the ability to co-evolve to recognize, respond, and adapt to the other. As such, microbial pathogens have developed various strategies to alter and evade host defenses in order to survive. Importantly, the host immune system is also adaptive and has a large arsenal to control or eliminate microbial threat. Even so, it is widely accepted that susceptibility of individuals within a given species differs to the same microbial pathogen. This variability in host–pathogen interaction is controlled by the inherent genetic make-up of the host, including innate and adaptive immune responses, particularly the acquired immunological memory, as well as the nature of the microbial pathogen (20).

Mastitis causing-bacterial pathogens are often well adapted to the bovine host resulting in clinical signs and, occasionally, subclinical infection before they lead to chronicity and persistence in the mammary gland. Persistent intramammary infections are frequently associated with recurrent clinical episodes and longterm increases in milk somatic cells counts. Persistent strains often express sets of genes that relate to their adaptation to the intramammary milieu and allow for intracellular survival and subsequent modulation of host-defense mechanisms (6, 21). *S. aureus* and *E*. *coli* are well-studied mastitis pathogens in the context of host–pathogen interaction and the elucidation of their genes, along with host immune response genes, is launching new studies in functional genomics (20). Understanding sequence data and locating functional SNPs in both the host and pathogen is expected to reveal relationships between immune function and the relevant genes that have the potential to advance resistance to specific pathogens.

Treatment of mastitis is given on the premise that treatment costs will be outweighed by production gains resulting from elimination of infection. Most farms have established mastitis management programs and include strategies such as routine whole herd antibiotic therapy, culling of chronically affected cows, post milking teat disinfection, as well as ensuring routine maintenance of milking machines (7, 14). Due to high treatment costs, lost income due to discarded milk, public health, and animal welfare concerns, it would be advantageous for dairy cattle to resist or mount effective immune responses to clear the wide variety of mastitis-causing pathogens. In the case of mastitis, the ability to control or tolerate the infection without actually clearing the pathogen, a phenomena known as resilience or tolerance (22), is not sufficient given that dairy products are consumed by human beings and are expected to be free of all potentially harmful pathogens. Antimicrobial treatment has the potential to increase the risk of antibiotic resistant strains of bacteria emerging in the environment (23), although it has been suggested that scientific evidence does not support emerging resistance in pathogens isolated from dairy cows (24). Nonetheless, other non-antibiotic treatment strategies are clearly warranted. Additionally, decreasing the incidence of mastitis would contribute to increased animal welfare as severe signs are associated with pain and discomfort for the cow (25).

Mastitis is a problem that plagues dairy cattle worldwide; however, this review will focus on the mastitis situation in the most economically developed countries. We highlight the genetic regulation of the bovine immune system and its vital contributions to disease resistance, in particular mastitis. Current genetic selection approaches used in the dairy industry to reduce the incidence of disease are reviewed, including the HIR technology; the Immunity+™ sire line, as well as genomics to improve disease resistance or immune response. While the complex interactions of the host and pathogen are fully acknowledged, they are only briefly discussed here.

#### **GENETIC REGULATION OF THE IMMUNE SYSTEM**

Robust, appropriate and timely host defense mechanisms are critical for prompt bacterial clearance and prevention of mastitis and mammary epithelial damage (14). Bacteria have a large repertoire of virulence factors that are produced at varying concentrations depending on the stage of infection (26), and these virulence factors in part determine differences in the magnitude and duration of host immune responses. Further, given the diversity of mastitiscausing pathogens, it is essential for the host to have a broad range of host-defense mechanisms as part of its immunological arsenal. Both innate and adaptive host defenses are required to protect the host from infection. Innate defenses against mastitis pathogens are rapid and include neutrophil recruitment to the mammary gland to facilitate bacterial clearance through phagocytosis, production of reactive oxygen species, antibacterial peptides, such as lactoferrin and β-lactoglobulin, and defensins, resulting in increases in the somatic cell count (18, 27). Mammary epithelial cells are known to play a role in early responses through the production of cytokines like IL-8 and other factors with antimicrobial activities (28, 29). If the bacteria survive these innate host defenses, adaptive immune responses mediated by T and B cells are required to clear the infection (30). The ideal immune response being one that appropriately recognizes epitopes on the invading pathogen to initiate swift and accurate clearance mechanisms while maintaining minimal pathological consequences. In some situations, such as experiments using *in vitro* or *in vivo* lipopolysaccaride challenge to measure bovine inflammatory responses, particularly IL-8, have noted that cows with lower IL-8 responses had quicker recovery in terms of somatic cell counts and milk production than those with high IL-8 production (31). This may relate to a more moderate inflammatory response generated in these low IL-8 responders. However, it is important to note that this does not mean that cows classified as low responders for other immune response mechanisms, particularly adaptive immune responses are advantageous. In fact, dairy cows classified as high responders (robust and balanced responses) for adaptive immune response traits have been demonstrated to have reduced disease incidence (32). The other thing worth noting in these experiments was the observation that the differences between high and low IL-8 responses seemed to be controlled by epigenetic effects (33). Epigenetic influences on bovine type 1 (Interferon-γ) and type 2 (IL-4) cytokine production have also been reported in cows classified as high or low responders based on their antibody and cell-mediated immune responses (34). Researchers are only beginning to dissect both the genetic and epigenetic mechanisms that control immunity.

Initiation and regulation of adaptive immune responses are critical to the resolution of infection. Cells of the innate immune system recognize conserved pathogen associated molecular patterns from the bacteria by binding pattern recognition receptors on antigen-presenting cells (APC) such as macrophages and dendritic cells (35). Such pattern recognition receptors include toll-like receptors (TLR) that are located on cell and endosomal membranes (27, 36). The association of a TLR with a pathogen associated molecular pattern initiates a downstream signaling cascade leading to the activation of transcription factors, such as NF-κβ, which enter the nucleus, bind target promoters, and may induce the production of cytokines and other endogenous mediators. The 10 mammalian TLRs are known to elicit unique responses through intracellular signaling pathways, which initiate inflammatory and antimicrobial processes to eliminate the pathogen (36, 37). For example, the recognition of lipopolysaccharide (LPS) from *E. coli* by TLR4, facilitated by additional proteins including CD14, LPS binding protein, and myeloid differentiation protein, is associated with production of TNF-α, IL-1β, IL-6, and IL-8. The lipoteichoic acid of Gram positive bacteria like *S. aureus* recognized by TLR2 is associated with only transient increases in TNF-α and IL-1β as well as IgG2 (27). It is well recognized that *E. coli* induces a stronger increase in the pro-inflammatory cytokines TNF-α and IL-1β compared to *S. aureus* (12, 27, 38), contributing to the severe clinical signs typically associated with *E. coli* mastitis as compared to *S. aureus* where the majority of cases go unnoticed. This draws attention to the fact that although the innate immune responses provide a first line of defense against invading microbial pathogens, including those that cause mastitis, and contours ensuing adaptive immune responses; innate responses have the potential to generate harmful pathology by driving inappropriate or soaring inflammatory cascades (31). These need to

be carefully considered and closely monitored when considering immunological interventions.

The major histocompatability complex (MHC) plays an essential role in the induction and regulation of immune responses (39). The bovine MHC, bovine lymphocyte antigen (BoLA), has been associated with resistance or susceptibility to mastitis (40–43), somatic cell count (42, 44, 45), and immune response (40, 41, 46). Genetic variation, such as single nucleotide polymorphisms (SNP) in other candidate genes associated with resistance or susceptibility to mastitis have been identified, including TLR4 (47, 48), TLR2, and caspase-recruitment domain 15 (49); IL-10 (50), osteopontin (51), IL-8 and its receptor CXCR1 (52–54), CCL2 and its receptor (55), as well as a variety of other genes (56). Other molecules important in host defense against mastitis-causing pathogens such as β-defensins have been identified and their complex genetic regulation is beginning to be understood (57). The feasibility of breeding for resistance based on one SNP or a combination of SNP depends on the degree of variation each SNP explains in resistance to mastitis. Since mastitis is a complex genetic trait a combination of many genes will ultimately be responsible for resistance to mastitis; however, certain major genes may contribute more benefit than others and it is important that these genes be elucidated.

Recent studies are beginning to uncover information about the epigenetic influences on bovine immune response genes (58). Some studies are now indicating that epigenetic changes are involved in the regulation of type I and II immune responses of mammals (59, 60), including cytokine profiles of dairy cows during the peripartum period when the risk of mastitis is the greatest (34). Epigenetic modifications have also been demonstrated to play a role in bovine innate immune responses to LPS stimulation (33, 61). Further, microRNA have been found to be differentially expressed upon challenge with mastitis-causing pathogens, suggesting a role for microRNA in regulating host responses to mastitis (62, 63). Indeed, many studies have demonstrated the bovine immune response to be under genetic and epigenetic control, and making use of this information in breeding strategies is anticipated to help improve udder health.

The important question is how to use this information regarding genetic associations with mastitis and the immune system to actually improve disease resistance. This is not necessarily a straight forward question given the plethora of genes, including their additive, dominant, epistatic, and epigenetic interactions. It is sometimes possible to make genetic gains in livestock health to a particular disease by selecting for or against a specific gene. Some examples of this include selection against Mareks Disease of poultry based on MHC haplotypes (64), bovine dermatopholosis using information on BoLA (65), brachyspina in cattle (66) among others (67). It is generally straightforward to make genetic gains for diseases caused by single recessive disorders, whereas information on single genes or clusters of genes may be less informative when trying to enhance resistance to complex traits, such as mastitis resistance, which is caused by a diverse set of pathogens controlled by a large variety of genes and gene interactions (68).

It is also worth noting that the immune system, which is the body's main host defense system, is regulated by thousands of genes (69). This points to the critical importance and complex nature of disease resistance as an overall fitness trait (70, 71). In fact, recent information from a human systems biology data base on immunity known as the immunogenetic-related information source – IRIS provides evidence for 1,535 immune response genes as of April 2013<sup>1</sup> . This list of genes was curated by IRIS with the following strict definition of a bona fide immune response gene, "a complete gene that produces a functional transcript and demonstrates at least one of the following defense characteristics: (i) known or putative function in innate or adaptive immunity, (ii) participates in the development or maturation of immune system components, (iii) induced by immunomodulators, (iv) encodes a protein expressed primarily in immune tissues, (v) participates in an immune pathway that results in the expression of defense molecules, (vi) produces a protein that interacts directly with pathogens or their products"<sup>2</sup> . When a broader definition of immune response genes are given that seeks to retrieve all genes that have some immune system or related functions, such as that provided by the Immunology Database and Analysis Portal (ImmPort), the list of genes is in the range of 6000<sup>3</sup> . Although these databases are based on human genes the newest version of the innate immunity database, InnateDB, does incorporate a list of bovine genes, including pathway and molecular interactions<sup>4</sup> . As pointed out by Karin Breuer and colleagues, as the experimental data from cattle research validates genetic interactions and immunological pathways this will allow for a deepened understanding of important bovine diseases, such as mastitis and tuberculosis (69). At the moment, these immunological databases rely largely on orthological-based approach to predict pathways. As of September 2012, the InnateDB contained more than 70,000 bovine interactions based on orthology and pathway analysis could assign to more than 7000 bovine genes (69). However, since the bovine immune system does contain some unique genetic features, such as a novel bovine type 1 interferon family known as IFNX, it will not always suffice to rely on orthogues from other species. Nonetheless, it is interesting to speculate about similar genetic pathways. For example, work in human beings has shown that following exposure to bacterial endotoxin a set of 3,714 unique genes were differentially expressed. These changes in genes of interest were confirmed in follow-up microarray experiments (72). Similar transcriptional changes might be predicted in cattle exposed to endotoxin from *E. coli* following intramammary exposure (73), as the complex plethora of genes involved in response to mastitis, such as that caused by *E. coli* is well known (74–76). The goal of this type of systems biology research is to provide a portrait of the entire "interactome between the innate and adaptive immune system, as well as its interconnection with other body systems in the hopes to enhance disease prevention and treatment strategies.

#### **GENETIC SELECTION FOR DISEASE RESISTANCE**

Current genetic selection approaches to improve mastitis resistance include both direct and indirect methods. With the exception of Nordic countries that have been selecting for disease resistance for over 35 years (77), most countries breed for mastitis resistance indirectly through SCC (78). More recently, France (79) and Canada (80) have launched routine national genetic and genomic evaluations for clinical mastitis. Problems associated with breeding directly for mastitis resistance include low heritability, the need for accurate health recording, and perhaps most importantly, the potential to skew the immune system causing individuals to be susceptible to other harmful pathogens. This skewing is thought to occur since antibody and cell-mediated immune responses are independent or slightly negatively correlated traits (81–84). This means that improvement for one of these traits does not translate into improvement of the other adaptive immune response trait. This concept will be discussed in more detail.

The heritability of mastitis resistance is low, with estimates ranging from about 0.02–0.10 (85, 86). SCS is genetically correlated (0.7) with mastitis and has a higher heritability of about 0.17, which is why it is used as an alternative trait to breed for resistance to mastitis (87–89). Divergent selection experiments based on SCS in sheep and cattle have been performed with the goal of creating lines of animals with an ability to resist intramammary infection (90, 91). Although these studies have shown a decrease in mastitis in the low SCS line, caution must be used in this approach to improve udder health. SCS tends to monitor subclinical cases (92) and although decreasing bulk tank counts has been associated with a decline in subclinical mastitis; clinical mastitis continues to be a problem (93). Further, since the cells that constitute the SCS are cells of the immune system, too low a SCS has been associated with an increased risk of clinical mastitis (94). In Canada, the approach will be to equally weight clinical mastitis and SCS in the LPI starting in August 2014. Other immune response traits known to associate with resistance to various diseases, including mastitis, may be added subsequently, although sires with improved immune responses are already available through the Canadian breeding company, the Semex Alliance since December 2012 (32).

In order to select directly for mastitis resistance, accurate disease records are essential. Many countries record disease on a voluntary basis, as is the current situation in the United States (86, 95) and Canada (85, 96). The use of voluntary producer records has brought into question the reliability of the estimates for disease resistance. By applying minimum lactation incidence rates to producer-recorded data to include only herds with regular recording, it has been found that although the heritability of disease resistance tends to be low (0.01–0.20) significant genetic variation exists to select for disease resistance (85, 95–97). Some research has demonstrated the use of genomics to improve the reliability of genetic estimates for disease resistance traits (86).

Selection against clinical mastitis has the potential to leave cattle susceptible to infection with other mastitis pathogens, since bacteria require unique immune responses for host protection (2), and mastitis pathogens have been demonstrated to change over time and geographically (7). Further, mastitis-causing pathogens tend to be extracellular in nature, requiring robust antibody responses (98). Since antibody- and cell-mediated immune responses tend to be negatively genetically correlated (83, 84) selection for mastitis resistance may potentially leave individuals with diminished capability to respond to intracellular pathogens generally controlled

<sup>1</sup>http://www.innatedb.com/curatedGenes

<sup>2</sup>http://www.innatedb.com/redirect.do?go=resourcesGeneLists

<sup>3</sup>http://www.immport.org

<sup>4</sup>http://www.innatedb.com

by the cell-mediated immune response. Cell-mediated responses have been demonstrated to be critical in controlling *Mycobacterium avium* spp *tuberculosis*, the causative pathogen associated with Johne's disease in cattle (99). Maintaining balanced immune responsiveness is an essential consideration in any breeding program to improve animal health. The other contributing factor is that different BoLA alleles have been shown to associate with antibody versus cell-mediated immune responses, as well as mastitis resistance (41). However, these are not the same alleles that associate with resistance to other viral or parasitic pathogens (100, 101). Therefore, caution must be exercised when selecting for resistance to one specific disease, particularly when it can be caused by multiple pathogens, as is the case with mastitis. Nonetheless, mastitis is such a costly disease that it is likely to be included in selection indices in conjunction with other health traits, such as SCS, until alternative approaches based on optimizing host defense mechanisms are more widely available. For example, in Canada an index for mastitis resistance was developed that includes both clinical mastitis and SCS traits and will be added to the Lifetime Profitability Index (LPI) in August 2014 (102, 103).

A combination of approaches is likely necessary to decrease mastitis occurrence,such as breeding for broad-based disease resistance based on immune response traits. Breeding for enhanced immune responsiveness is a solution to provide cows with an overall superior ability to respond to a variety of pathogen types requiring unique responses to provide broad-based disease resistance. Individuals with greater and optimally balanced antibody and cell-mediated immune responses breeding values are referred to as high immune responders (HIR) (**Figure 1**) and the method for identifying such individuals is referred to as the HIR technology (32, 104).

The HIR technology has been used to identify the ability of cows, calves, and bulls to mount antibody and cell-mediated adaptive immune responses (106, 107). These adaptive immune response traits are heritable, on average 0.25–0.35 (83, 84), considerably higher than estimates for specific clinical or subclinical disease resistance (**Table 1**). The heritability of immune response is similar to what has been found for milk production traits, indicating it would be possible to make significant genetic gain depending on how heavily health is weighted within the selection index. Cows with superior adaptive immune responses have been demonstrated to have substantially lower occurrence of diseases, including mastitis, metritis, displaced abomasums, retained fetal membranes (108) and are less likely to be seropositive for *Mycobacterium avium* spp *paratuberculosis* (109). It would, therefore, be feasible and desirable to breed dairy cows for enhanced immune responses to decrease the occurrence of diseases like mastitis (100, 110). Previously, this approach was shown to improve disease resistance of pigs (105). It should also be noted that producing robust adaptive immune responses requires appropriate priming via particular innate host defense pathways, such as TLR signaling (37). Priming the immune system with LPS in the udder has been shown to reduce bacterial load in experimentally induced mastitis via the TLR signaling (111, 112).

High immune responding cows have also been found to have an increased response to commercial *E. coli J5* mastitis vaccination (117), as well as improved colostrum quality as measured by

host immune response phenotype is ultimately determined by the interaction of the immune response genotype with the environment. The expression of the immune response genotype is also regulated by epigenetic effects. The innate immune response is relatively fast acting and non-specific, but is critical to signal appropriate adaptive cell-mediated and antibody-mediated immune responses. Dairy cattle with enhanced and balanced cell and antibody-mediated immune responses are known as high immune responders.

#### **Table 1 | Heritability estimates of immune response, mastitis resistance, and milk production and in Holstein dairy cattle**.


specific antibody (117), total immunoglobulin, lactoferrin, and βlactoglobulin (118). Differences in leukocyte populations between high and low immune responders have also been described, such that cows with superior antibody responses have a higher proportion of B cells in peripheral blood in response to immunization, whereas cows with high cell-mediated responses have a higher baseline proportion of gamma delta (γδ) T cells (119). These differences in the diverse phenotypes identified using the HIR technology suggest potential mechanisms that contribute to decreased disease occurrence among high immune responding individuals.

Multiple studies over many years have found beneficial associations between antibody responses and a lower occurrence of mastitis. A study that evaluated antibody-mediated immune responses to a specified test antigen found cows with superior antibody responses had lower occurrence of mastitis in two out of three herds tested (117). Subsequently, cows with greater antibody responses in a commercial herd in Florida were found to be 1.6–2.5 times less likely to get clinical mastitis compared to other cows in the herd (108, 120). Most recently, a nation-wide study in Canada evaluating the incidence rate of clinical mastitis over a 2-year study period found cows with superior antibody responses to have an incidence rate of 17.1 cases of clinical mastitis/100 cow years compared to average and low responding cows with 27.9 and 30.7 cases, respectively. The low responding cows were also found to have more severe mastitis compared to cows with better immune responses (98). Antibody-mediated immune responses have also been beneficially genetically correlated with some reproductive traits as well as longevity, suggesting that cows with better immune responsiveness and therefore, less disease remain in the herd longer (83).

Conversely, cows with greater cell-mediated immune responses have been found to be less likely to be seropositive for *Mycobacteria avium paratuberculosis* (109). Cell-mediated immune responses are also critical to provide protection against *S. aureus* small colony variants that can cause mastitis and have the ability to survive within host cells (6, 21). Antibody and cell-mediated immune responses have been found to be negatively genetically correlated (83, 84). Consequently, in order to ensure protection to a broad range of pathogens it is essential to identify and select individuals with the capacity to generate both effective antibody and cell-mediated immune responses (32).

The Semex Alliance utilizes the HIR Technology to identify dairy sires with superior immune responsiveness, termed *Immunity*+™*.* Daughters of *Immunity*+™ sires have been found to have lower disease occurrence and higher profitability compared to daughters of sires with either an unknown or an average or low immune response type. Specifically, daughters of *Immunity*+™ sires in a large herd in the US had a 44% reduction in mastitis, 25% less calf pneumonia, and an 8.5% reduction in all diseases in first lactation heifers (32). These results highlight the benefit and potential to improve disease resistance, in particular mastitis resistance, by improving overall immune responsiveness.

Genomic selection has allowed for the opportunity to include new phenotypes in breeding objectives, particularly those that may be relatively expensive to measure (121). Genomic selection refers to breeding decisions based on genomic estimated breeding values (GEBV), which are calculated using the joint effects of SNP markers across the entire genome (122–124). Using a large reference population with accurate phenotype information, the SNP or haplotype effects for a given trait are estimated. In subsequent generations, only information on the SNP or haplotypes are required to calculate the GEBV (123). Genomic selection has provided many substantial benefits to the dairy industry. Perhaps the most highlighted benefit is in the significant increase in the rate of genetic gain by decreasing generation interval, increasing, and selection intensity the accuracy of estimates (122).

The sequencing of the bovine genome and release of SNP arrays used for genomic selection has led to increases in the genome-wide association studies (GWAS). Many GWAS have been performed, which has lead to the identification of quantitative trait loci or SNP profiles associated with resistance or susceptibility to mastitis (125), or SCC as an indicator of mastitis (126–128). Using the approach, many genes involved in immune response have been found, including cytokines IL-4 and IL-13 as well as IL-17 (129). Recently, a series of GWAS have been performed for general immune responsiveness in dairy cattle and results have been validated in dairy sires (46). Results of this work have identified many genes associated with immune responses including the bovine MHC, the complement systems as well as cytokines including IL-17 and TNF in the genetic regulation of bovine immune system. Results of these GWAS on mastitis resistance and immune response suggest that it is possible to calculate GEBV for mastitis or immune response traits increasing the accuracy of estimates for genetic selection. The next critical steps are to create large reference populations with genotypes and accurate phenotypes for disease and immune response traits in order to improve dairy cattle health.

#### **CONCLUSION**

The ideal solutions to improve resistance to mastitis are likely to be those that focus on a large number of genes, by using information from GWAS, or selection based on breeding values of immune responses, which take into account complex genetic interactions between the innate and adaptive host defense mechanisms without the necessity of knowing all about each individual gene. Using selection indices also offers the advantage of being able to easily adjust the weights given to the various traits within the index as the selection proceeds. These two approaches may be best suited to help alleviate mastitis, at least until we gain more knowledge about genetic and epigenetic regulation of host defense mechanisms.

#### **ACKNOWLEDGMENTS**

The authors gratefully acknowledge the agencies that have provided research grants to Bonnie A. Mallard including the Semex Alliance, the Natural Sciences and Engineering Research Council (NSERC), the DairyGen Council of the Canadian Dairy Network (CDN), and the Ontario Ministry of Agriculture, Foods and Rural Affairs (OMAFRA).

#### **REFERENCES**


stages of a bacterial infection to *Staphylococcus aureus*. *Vet Res* (2014) **45**:16. doi:10.1186/1297-9716-45-16


polymorphisms and milk somatic cell score in Canadian Holsteins, and functional relevance of SNP c.3020A>T. *Dev Biol (Basel)* (2008) **132**:247–53. doi:10.1159/000317167


*10th International Colloquium on Paratuberculosis*, Vol. 1. Minneapolis (2009). p. 127.


129. Lewandowska-Sabat AM, Gunther J, Seyfert HM, Olsaker I. Combining quantitative trait loci and heterogeneous microarray data analyses reveals putative candidate pathways affecting mastitis in cattle. *Anim Genet* (2012) **43**:793–9. doi:10.1111/j.1365-2052.2012.02342.x

**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.

*Received: 29 July 2014; accepted: 23 September 2014; published online: 07 October 2014.*

*Citation: Thompson-Crispi K, Atalla H, Miglior F and Mallard BA (2014) Bovine mastitis: frontiers in immunogenetics. Front. Immunol. 5:493. doi: 10.3389/fimmu.2014.00493*

*This article was submitted to Molecular Innate Immunity, a section of the journal Frontiers in Immunology.*

*Copyright © 2014 Thompson-Crispi, Atalla, Miglior and Mallard. This is an openaccess article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.*

### ADVANTAGES OF PUBLISHING IN FRONTIERS

FAST PUBLICATION Average 90 days from submission to publication

COLLABORATIVE PEER-REVIEW

Designed to be rigorous – yet also collaborative, fair and constructive

RESEARCH NETWORK Our network increases readership for your article

#### OPEN ACCESS

Articles are free to read, for greatest visibility

#### TRANSPARENT

Editors and reviewers acknowledged by name on published articles

GLOBAL SPREAD Six million monthly page views worldwide

#### COPYRIGHT TO AUTHORS

No limit to article distribution and re-use

IMPACT METRICS Advanced metrics track your article's impact

SUPPORT By our Swiss-based editorial team

EPFL Innovation Park · Building I · 1015 Lausanne · Switzerland T +41 21 510 17 00 · info@frontiersin.org · frontiersin.org