Variation in the Ovine Abomasal Lymph Node Transcriptome between Breeds Known to Differ in Resistance to the Gastrointestinal Nematode

Texel lambs are known to be more resistant to gastrointestinal nematode (GIN) infection than Suffolk lambs, with a greater ability to limit infection. The objectives of this study were to: 1) profile the whole transcriptome of abomasal lymph node tissue of GIN-free Texel and Suffolk lambs; 2) identify differentially expressed genes and characterize the immune-related biological pathways and networks associated with these genes. Abomasal lymph nodes were collected from Texel (n = 6) and Suffolk (n = 4) lambs aged 19 weeks that had been GIN-free since 6 weeks of age. Whole transcriptome profiling was performed using RNA-seq on the Illumina platform. At the time of conducting this study, a well annotated Ovine genome was not available and hence the sequence reads were aligned with the Bovine (UMD3.1) genome. Identification of differentially expressed genes was followed by pathway and network analysis. The Suffolk breed accounted for significantly more of the differentially expressed genes, (276 more highly expressed in Suffolk v 162 in Texel; P < 0.001). The four most significant differentially expressed pathways were all related to immunity and were classified as: Role of Pattern Recognition Receptors in Recognition of Bacteria and Viruses, Activation of IRF by Cytosolic Pattern Recognition Receptors, Role of RIG-I-like Receptors in Antiviral Innate Immunity, and Interferon Signaling. Of significance is the fact that all of these four pathways were more highly expressed in the Suffolk. These data suggest that in a GIN-free environment, Suffolk lambs have a more active immune profile relative to the Texel: this immune profile may contribute to the poorer efficiency of response to a GIN challenge in the Suffolk breed compared to the Texel breed.


Introduction
Gastrointestinal nematode (GIN) infection of ruminants is a major economic and health concern, causing substantial loss to livestock producers worldwide. Increasing consumer concerns about drug residues in animal products and the emergence of anthelmintic-resistant nematode species have stimulated efforts to develop alternative strategies to anthelmintic therapy [1], including genetic selection for natural resistance to nematodes. The latter approach is justified by evidence of genetic variation in resistance to GIN infection both within and between sheep breeds [2][3][4][5][6].
The identification of genetically resistant animals through DNA analysis has been focused, to date, on the identification and characterization of candidate genes such as the major histocompatiblity complex (MHC) DRB [5,7], MHC class I and II genes [8,9] and interferon gamma [10,11], all of which have alleles that have been associated with resistance to GIN. A weakness of the candidate gene approach is that resistance/susceptibility to GIN is a complex trait with genetic control being polygenic. Our hypothesis is that biochemical pathways and networks are central to differences in resistance and, thus, identification of these will provide a more secure foundation for the elucidation of gene profiles associated with resistance to GIN. Others have shown that a number of genes/biological processes are differentially expressed in duodenal tissue from two lines of lambs selected for divergence in resistance to GIN [12,13].
Previous studies from our group have shown that the Texel breed is more resistant to GIN infection than the Suffolk breed [14]. Comparing the gastrointestinal lymph tissue transcriptome between Texel and Suffolk in the absence of a GIN challenge should allow the identification of gene pathways and networks that are differentially expressed between the two breeds and thus likely to be involved in the response to a GIN challenge. An understanding of the molecular basis for such breed variation may reveal markers associated with resistance to pathogens. Hence, the objectives of the present study were: 1) to profile the whole transcriptome of abomasal lymph node tissue of GIN-free Texel and Suffolk lambs; and 2) to identify the differentially expressed genes and characterize the immune-related biological pathways and networks associated with the difference in resistance between the Texel and Suffolk breeds. As a well annotated Ovine genome was not released at the time of conducting this study we aligned the transcriptome with that of the UMD3.1 Bos taurus genome.

Ethics statement
All procedures described in this experiment were conducted under experimental license (B100/ 2584) from the Irish Department of Health in accordance with the Cruelty to Animals Act 1876 and the European Communities (Amendments of the Cruelty to Animals Act 1976) Regulations, 2002 and and approved by the Teagasc Ethics Committee (5909).

Animals
The lambs described here were part of a larger experiment concerning the response of Suffolk and Texel lambs to experimental infection with Teladorsagia circumcincta [15]. Lambs (10 in total; 6 Texel and 4 Suffolk) were born (28th and 29th March) indoors, and moved to pasture until 6 weeks of age (mid-May), at which point they were weaned and moved and kept indoors until mid-September. Given that Nematodirus battus was the predominant species observed in FEC in lambs after housing (all lambs, data not shown), it was the predominant species on pasture. Animals were maintained indoors on a concentrate-based diet, with free access to water, from mid-May to mid-September. The lambs were faecal sampled per rectum at weaning but sufficient material was obtained from only 7 individuals (4 Texel and 3 Suffolk). Nematodirus eggs were detected in all cases (17,50,150, 1700 e.p.g. for Texel and 22, 300, 1000 e.p.g. for Suffolk) while eggs from other Trichostrongyle spp (excluding S. papillosus) were detected in 4 cases (1, 17, 100 e.p.g. for Suffolk and 50 e.p.g. for one Texel). All lambs were then treated with Ivermectin (Oramec, Merial Animal Health Limited) according to the manufacturer's instructions and quarantined in a slatted pen for 48 h prior to being penned on straw. Faecal samples were collected on 3 consecutive days 5 weeks post-housing/anthelmintic treatment to determine GIN infection status. No eggs were detected in the faecal material of any of the 10 lambs at this stage. The lambs were faecal sampled again one week prior to slaughter to confirm their trichostrongyle infection-free status for at least 13 weeks prior to slaughter.

Tissue samples
At the end of the experiment, all ten animals were slaughter by electrical stunning followed by exsanguination. Immediately after slaughter, animals were dissected and approximately 1 g of abomasal lymph node was stored from each animal at room temperature for 24 h in 10 mL RNAlater (Ambion, Austin, TX, USA). RNAlater was then discarded and samples were stored at -80°C until RNA was isolated.

RNA extraction
A representative subsample (40 to 50 mg) of the abomasal lymph node tissue was homogenized in 1 mL of TRIreagent (MRC, Cincinnati, OH, USA), using the Qiagen TissueLyzer, and total RNA was extracted following the manufacturer's instructions. DNase I (Sigma-Aldrich, St Louis, MO, USA) treatment of the RNA was then performed and further purification was carried out using the GenElute mammalian total RNA miniprep kit (Sigma-Aldrich, St Louis, MO, USA). The quantity and quality of total RNA were assessed using a Nano Drop spectrophotometer (Thermo Fisher Scientific, USA) and an Agilent 2100 Bioanalyzer (Agilent Technologies, Inc. CA, USA), respectively. All RNA samples used in this study had an RNA integrity value 8.0.

RNA-seq library preparation
An RNA-seq library was prepared from~2 μg of total RNA using an Illumina TruSeq RNA Sample Prep Kit (Cat. No FC-122-1001; Illumina Inc., San Diego, CA, USA) according to the manufacturer's protocol. A multiplexing-capable kit was used for sample multiplexing. The final concentration of the libraries was assessed using a Qubit fluorometer (Invitrogen, Paisley, UK) and quality was assessed using a high sensitivity DNA kit (Agilent Technologies, Inc. CA, USA). The high sensitivity DNA chip and the Agilent 2100 Bioanalyzer were used for library sizing and final validation. One RNA-seq library was prepared from abomasal lymph node tissue from each animal.

Sequencing
Equi-molar quantities (10 nM) of the bar-coded cDNA libraries were multiplexed and run in 6 lanes across three flow cells. The sequencing was carried out using Illumina's Genome Analyzer II (Illumina Inc., San Diego, CA, USA) according to the manufacturer's instructions. The sequencing products consisted of single-end reads of 36 nucleotides plus the 6-nucleotide index marker (42 bases in total). FastQC [16] software was used to assess the quality of the sequence data; sequences with a Phred mapping quality score 30 were used for further analysis. This quality score threshold corresponds to base call accuracy of 99.9% [17]. As the quality of the sequence data was very good we did not trim the raw sequences during the processing of data, and thus the length of raw and processed reads were same.

Data analysis
As a well annotated Ovine genome was not released at the time of conducting this study we aligned the transcriptome with that of Bovine. The alignment software Bowtie 0.12.7 [18] was used to align the reads to the UMD3.1 Bos taurus genome assembly allowing up to two mismatches per read. Only uniquely aligned reads were used for the analyses. The binary alignment/map (BAM) files from the Bowtie mapping were sorted and filtered for duplications (possibly resulting from PCR bias) using SAMtools [19]. The raw counts, per gene were estimated by HTseq (v0.4.7) (http://www-huber.embl.de/users/anders/HTSeq/doc/overview. html#). The data were normalized in edgeR [20] using the trimmed mean of M values (TMM) method [21]. Only genes with 5 reads in total, across all samples, were included in the final analysis. The normalized raw counts per gene were used to identify genes that were differentially expressed (DE) between breeds using edgeR. Gene ontology (GO) of DE genes was performed using GO-Elite [22]. Classification of functional processes of biological importance and of canonical pathways and networks of DE genes were performed using Ingenuity Systems Pathway Analysis (IPA; Ingenuity Systems, Redwood City, CA, USA; http://www.ingenuity. com). Data has been made publically available: http://www.ncbi.nlm.nih.gov/geo/query/acc. cgi?acc=GSE43241.
As some of the most significant differentially expressed pathways suggested that the animals may have been infected with a virus, the RNA-seq data were further aligned with the reference genomes of a number of different viral species [orf virus (HM133903.1) virus, jaagsiekte sheep retrovirus (NC_001494.1), Schmallenberg (NC_018459.1), bluetongue virus (NC_006023.1)]. Moreover, the findings of alignment (presence of orf virus) were further confirmed by using conventional polymerase chain reaction. Genomic DNA was prepared from 100 mg of abomasal lymph node tissue samples. Amplification of orf virus was carried out according to Inoshima et al. [23]. The forward-(5'GCGAGTCCGAGAAGAATACG3') and reverse (5'TACGTGG GAAGCGCCTCGCT3') primers bind to a target region specific to the major envelope protein (B2L) gene of orf virus isolate SV178/12 and yielding a PCR product of 594 bp. The reaction conditions and PCR cycling conditions are as described by Inoshima et al. [23].

Validation by quantitative real time PCR (qPCR)
To validate the results of RNA-seq analysis, a panel of differentially expressed genes was selected from the RNA sequence database and qPCR assays were designed to confirm the direction and magnitude of the expression profiles. Ten genes were randomly selected from the panel of 27 DE genes identified and classified as part of the genetic network 'Infectious Disease, Antimicrobial Response, Inflammatory Response'. Gene specific primers were designed using Primer Express (v2.0) (PE Applied Biosystems) software and were synthesized by Sigma Aldrich, USA. The details of the primers are in Table 1.
Complementary DNA (cDNA) was synthesized from 1 μg of purified total RNA using cDNA Synthesis Kit (Fermentas, Vilnius, Lithuania) and oligo dT primers following the manufacturer's instructions. The primer efficiency was determined using a serial dilution (1:4 dilution series over 7 points) of a cDNA pool, prepared by pooling an equal quantity of cDNA from all of the experimental samples; the efficiency of all primers was shown to be between 90% to 110%. The reaction was carried out in a total volume of 20 μL, which comprised 10 μL 2X Fast SYBR green PCR Master Mix (Life Technologies), 1 μL of forward and reverse primer mix (final concentration 300 nM each), 1 μL of diluted cDNA template (equivalent to 2.5 ng of RNA) and 8 μL of nuclease-free water (Sigma-Aldrich, St Louis, MO, USA), in an ABI Prism 7500 FAST sequence detection system (Applied Biosystems, Warrington, UK). The PCR cycling conditions were 95°C for 10 min for 1 cycle followed by 95°C for 3 sec and 60°C for 30 sec for 40 cycles. The specificity of the product was confirmed by melt curve analysis.
Following analysis of the relative quantities of the gene transcripts of 8 potential reference genes (GAPDH, B2M, PGK1, ATPsynth, RPL19, TBP, GUSB, ACTB) with the GeNorm application within the qBase PLUS software package [24] (Biogazelle, Zwijnaarde, Belgium), two robust reference genes (PGK1 and RPL19) were selected for normalization of the relative expression of the target genes. Calibrated normalized relative quantities of gene expression for each target gene were generated using qBase PLUS . RNAseq and qPCR data were compared by calculating the correlation coefficient for specific genes after accounting for breed effects using SAS v9.1.3 (SAS Institute, Cary, NC, USA).

Results and Discussion
An average of 13 (s.d. 1.9) million reads were generated per sample. Due to the lack of a wellannotated Ovine genome at the time of conducting this study, reads were aligned to the bovine sequence, which was the closest species with a well-annotated genome. Alignment of these reads to the Bos taurus genome yielded an average of 7 (s.d. 1.1) million uniquely aligned reads per sample, which were subsequently used in the analysis.

Differentially expressed genes and validation of RNA-seq data
Based on a false discovery rate < 0.05, a total of 437 DE genes were identified. The gene expression patterns, in terms of direction and magnitude, of all 10 genes chosen for validation were reproducible by qPCR analysis (Fig 1). Of the DE genes, 276 were more highly expressed in Suffolk while 161 were more highly expressed in Texel. The proportion of DE genes was significantly associated with breed (P < 0.001). Gene-ontology analysis of differentially expressed genes Functional annotation classification of the 276 genes more highly expressed in the Suffolk breed yielded a total of 59 significant (P< 0.05) functional annotation clusters; the 10 most significant clusters are presented in Table 2. The majority of the functional annotation clusters in the Suffolk relate to the immune system, including: ISG15-protein conjugation, Immune response, Antigen processing and presentation of peptide antigen via MHC class I, Response to other organism, Defence response, MHC class I protein complex, Regulation of interleukin-10 production, Chemokine activity and Type I interferon biosynthetic process. The most significant group, ISG15-protein conjugation, contains 4 genes, whereas the two largest clusters Immune response and Defence response comprise 31 and 21 genes, respectively. In the Texel breed, a total of 26 significant (P < 0.05) functional annotation clusters were identified and the 10 most significant clusters are presented in Table 3. The GO terms associated with highly expressed genes in the Texel breed describe normal cellular processes and include: Vitamin transport, Regulation of smooth muscle cell proliferation, Collagen fibril organization, Proteinaceous extracellular matrix, Extrinsic to plasma membrane, Extracellular space, Glucose metabolic process and Cell-cell adhesion.

Identification of pathways and gene interaction networks
Pathway analysis. A total of 23 significantly over-represented canonical pathways were identified as differentially expressed between GIN-free Texel and Suffolk sheep ( Table 4). The four most significant pathways all relate to processes of anti-viral /anti-bacterial innate immunity and were all more highly expressed in the Suffolk lambs relative to the Texel lambs. The four pathways and the DE genes assigned to each pathway are as follows:  • Role of RIG-I-like Receptors in Antiviral Innate Immunity: RIG-I, MDA-5, TBK1, IRF3, TRIM2, (Fig 4); • Interferon Signaling: IFIT1, IFIT3, IRF9, OAS1, MX1, STAT, (Fig 5).
Gene interaction networks. The overall gene interaction networks for the differentially expressed genes were constructed using IPA. A total of 25 significant (P < 0.05) gene interaction networks were identified. The details of the 10 most significant networks are presented in Table 5. The gene interaction network with the highest number of focus molecules (n = 27) was designated Infectious Disease, Antimicrobial Response, Inflammatory Response. This network included 24 genes that were more highly expressed in the Suffolk breed and 3 genes that were more highly expressed in the Texel breed (Fig 6).
Exploring the roles of these specific genes, pathways and the network in more detail, it is evident that in a GIN-free environment, the Suffolk lambs had greater expression of a number of cytosolic pattern recognition receptors, membrane bound receptors, transcription factors and IFN-induced proteins than Texel lambs. Interestingly, a number of DE genes are assigned to more than one of these pathways (e.g., RIG-I, MDA-5, STAT2, IRF-3, IRF-9, OAS1), indicating that there is overlap between the pathways. These molecules are normally associated with the T helper cell type 1 (Th1) response to viral and/or bacterial exposure [25][26][27][28]. A variety of pattern recognition receptors (PRRs) were more highly expressed in the Suffolk lambs relative to the Texel lambs. PRRs are a primitive component of the innate immune system. They recognise microbe-specific molecules (pathogen associated molecular patterns; PAMPS), including viral RNA/DNA, lipopolysaccaride, mannose, bacterial peptides, peptidoglycans, lipoproteins and fungal glucans [26][27][28][29]. They can be located intracellularly, extracellularly, or membrane bound. The gene expression of three cytosolic PRRs (RIG-I, MDA-5 and DAI) was higher in the Suffolk lambs than in Texel lambs. Both RIG-I and MDA-5 belong to the RIG-1-like receptor family, whose members act as sensors of viral RNA. As represented in Fig 3, upon binding to viral RNA, these receptors induce type-1 interferon (IFN) production via serine/threonineprotein kinase (TBK1), interferon regulatory factor-3 (IRF-3) and ISG-15 [26]. The PRR DAI binds to microbial dsDNA, which induces type-1 interferon production via TBK1 and IRF-3 [27]. Hence, activation of any of these three PRRs leads to the production of IFN, which in turn activates a cluster of IFN-induced proteins, including a number of those genes that were more highly expressed in the Suffolk lambs: IFIT1, IFIT3 and MX1 which are activated via the transcription factors STAT2 and IRF9 [28,30,31]. It has previously been reported that the IRF9 protein binds with STAT2 protein to produce a set of IFN-induced proteins IFIT3, IFIT1, OAS1, OAS2 and MX1 to activate non-specific anti-viral immunity [32,33]. In fact polymorphisms in the OAS1 gene are associated with host susceptibility to various diseases [34,35] and with resistance to viral infections [36]. The genes coding for two membrane-bound PRRs (DECTIN-2 and C3AR1) were also more highly expressed in the Suffolk. DECTIN-2 is a PRR that stimulates non-specific antifungal immunity [37,38] and upon activation mediates a Th2 immune response [39]. There was, however, no evidence that the downstream effector molecules of DECTIN-2 (ERK, NFkB or JNK) were differentially expressed between the breeds (Table 4). C3AR1 is a receptor of the complement cascade [40,41], and interestingly, C1QB and C1QC, which are the B and C chains of the first component of the classical complement pathway were more highly expressed in the Suffolk lambs.
The spectrum of genes in the Th1/IFN network that are being transcribed in the Suffolk breed, suggests that the Suffolk animals are either innately more prepared for a microbial challenge than the Texel animals, or were actually actively responding to a microbial challenge. In exploring the possibility of recent microbial challenge, two approaches were adopted. Firstly, the haematological parameters of both Suffolk and Texel animals were assessed, and secondly, the RNA-seq data was searched for evidence of viral RNA from viral species known to infect sheep, including Schmallenberg virus, bluetongue virus, orf virus and jaagsiekte sheep retrovirus (JSRV). There was no evidence for an active infection based on the haematological parameters evaluated (unpublished data). Interestingly, there was evidence for two types of viral RNA in all the animals with a similar average number of reads in both breeds (JSRV: Suffolk 62±34, Texel 112±32; orf virus: Suffolk 2040±121, Texel 2190±152). It has previously been reported that the exogenous infectious form of JSRV has an endogenous counterpart, which was integrated into the ovine genome prior to the evolution of sheep and goats [42]. The sheep genome has approximately 27 integrated copies of endogenous JSRV that are genomically closely related to JSRV [43]. It is possible that the Suffolk lambs are responding to the products from the integrated JSRV genome and reacting in a manner described for autoimmune diseases [44]. However, the high number of reads aligned with the orf virus genome is highly suggestive of recent or active orf infection. Orf virus specific polymerase chain reaction in the sheep genomic DNA samples from the lymph node showed the presence of orf virus in all the animals (Fig 7). Viral DNA may be integrated with the genome of the animals and carried over within the cell during the life of the animals All of the animals used in this study received a fully virulent live virus vaccine (Scabivax, Intervet Ireland) during the first week of life. A recent study detected orf viral DNA in reindeer lymph nodes 4 weeks post experimental inoculation [45]. Within the confines of the current experiment, it was not possible to determine if either of these viruses are responsible for the immune profile in the Suffolk lambs. However, it is worth considering in our future studies particularly the alignment of this sequence data with Ovine genome are likely to reveal greater details on the nature of adaptive and innate immune system of these two breeds.
In conclusion, the analysis of the RNA transcriptome of Suffolk and Texel lambs maintained under GIN-free conditions, suggests that Suffolk sheep have a more active antiviral/antibacterial immune profile than Texels. It is likely that this profile contributes to the variation between the two breeds in their capacity to mount a timely and effective immune response when exposed to a gastrointestinal nematode challenge.  BPI, RIG-I, EIF2AK2, ENPP1, EPCAM, HERC5, IFI6, IFI27,  Ifi47, MDA-5, IFIT1, ISG-54, IFIT3, IFIT5, IRF3, IRF9, IRG, ISG15,  MX1, OAS1, OAS2