Possible mechanisms of host resistance to Haemonchus contortus infection in sheep breeds native to the Canary Islands

Haemonchus contortus appears to be the most economically important helminth parasite for small ruminant production in many regions of the world. The two sheep breeds native to the Canary Islands display distinctly different resistant phenotypes under both natural and experimental infections. Canaria Hair Breed (CHB) tends to have significantly lower worm burden and delayed and reduced egg production than the susceptible Canaria Sheep (CS). To understand molecular mechanisms underlying host resistance, we compared the abomasal mucosal transcriptome of the two breeds in response to Haemonchus infection using RNAseq technology. The transcript abundance of 711 and 50 genes were significantly impacted by infection in CHB and CS, respectively (false discovery rate <0.05) while 27 of these genes were significantly affected in both breeds. Likewise, 477 and 16 Gene Ontology (GO) terms were significantly enriched in CHB and CS, respectively (P < 1.0 × 10−4). A broad range of mechanisms have evolved in resistant CHB to provide protection against the parasite. Our findings suggest that readily inducible acute inflammatory responses, complement activation, accelerated cell proliferation and subsequent tissue repair, and immunity directed against parasite fecundity all contributed to the development of host resistance to parasitic infection in the resistant breed.


Results
Parasitology. The total worms recovered from the infected groups of CHB and CS were 1,109.75 (± 1,547.73, SD) and 3,280.50 (± 2,398.03), respectively. The difference is statistically significant (P < 0.05, Fig. 1). Neither Haemonchus worms nor fecal eggs were recovered from the uninfected group of either breed, as expected. EPG values detected from infected CS sheep were 262.50 ± 287.54 (mean ± SD) while no fecal eggs were detectable in the infected group of CHB sheep at 20 days post infection (dpi). No parasite eggs in either group prior to the experimental challenge were observed.
Haemonchus infection induced distinctly different transcriptome patterns in the abomasal mucosa of CHB and CS breeds. In this study, approximately 79.91% of raw reads (± 7.08%; SD) were uniquely mapped to the ovine genome. Compared to their respective uninfected controls, the numbers of genes significantly impacted by infection in CHB and CS breeds at a stringent cutoff value (false discovery rate or FDR < 0.05), were 711 and 49, respectively (Fig. 2). The abundance of 27 genes was significantly changed by infection in both breeds (Table 1). Among them, 25 genes, such as arachidonate 15-lipoxygenase (ALOX15), collagen, type VI, α 5 (COL6A5), and serglycin (SRGN), were significantly upregulated while the expression of transthyretin (TTC) was repressed by infection. Intriguingly, the transcript abundance of cadherin 26 (CDH26) was significantly induced by infection in both breeds (adjusted P value or FDR < 1.63 × 10 −10 ); and is strongly correlated with worm counts only in CS (Fig. 3). However, infection had a bidirectional impact on the transcript abundance of a uncharacterized gene containing a unknown microRNA (ENSOARG00000023771), which was significantly upregulated in CHB but downregulated in CS. The genes significantly impacted by infection only in CS included mast cell proteinase-3, γ -glutamyltransferase 5 (GGT5), CD163 as well as those involved in smooth muscle contraction, such as tropomyosin (TPM2), myosin, light chain 9, regulatory (MYL9), and calponin 1, basic, smooth muscle (CNN1). Among the genes significantly impacted by infection in CHB sheep, several cytokine receptors and chemokines were strongly upregulated. Notable, the transcript of IL17 receptor beta (IL17RB) was 14.4 fold higher in infected animals than uninfected controls in CHB. IL2 receptor beta (IL2B) was also upregulated. Similarly, chemokine CXC ligand 12 (CXCL12) and chemokine (CXC motif) receptor 6 (CXCR6) were upregulated by infection in CHB. Among the well-known Th2 cytokines, the expression of IL6, IL10 and IL13 was upregulated by infection in both breeds. Moreover, while the extent of upregulation of IL6 by infection remained similar in both breeds (~6.8 fold), overexpression of both IL10 and IL13 mRNA molecules was more profound in the resistant breed (CHB) than in CS. On the other hand, the IL5 mRNA was upregulated by infection in CS but barely detectable in CHB at the sequencing depth in this study. The IL4 expression followed the similar trend: it was upregulated approximately 9 fold by infection in CS but was barely detectable in CHB. However, the IL9 mRNA level remained unchanged by infection in both breeds.
Of note, approximately 15% of the genes significantly impacted by infection are cell-cycle related. The expression of these cell cycle related genes was predominantly enhanced by Haemonchus infection in CHB. As Table 3 shows, at least 92 genes were significantly upregulated by infection, such as cyclin A2 (CCNA2), cyclin B3 (CCNB3), various centromere proteins (CENPL, CENPN, CENPT, and CENPW) and kinesin family (KIF) members, and at least 5 minichromosome maintenance complex (MCM) components (MCM3, MCM4, MCM5, MCM6, and MCM10). Nevertheless, the infection was also able to repress cell cycle related genes, such as cyclin G1 (CCNG1), regulator of cell cycle (RGCC), and synaptonemal complex protein 3 (SYCP3). Moreover, at least five transcription factors, such as the oncogene MYB, SMAD family members 6 and 9 (SMAD6 and SMAD9), and histone decetylase 5 (HDAC5), were significantly affected by infection in CHB.
Intriguingly, four genes known to regulate abomasal acid secretion and gastric function 16 were downregulated by Haemonchus infection in CHB, including ATPase, H+ /K+ exchanging, alpha polypeptide (ATP4A), progastricsin (pepsinogen C, PGC), appetite-regulating hormone precursor (GHRL), and forkhead box A2 (FOXA1). However, the transcript abundance of these four genes remained unchanged by infection in CS.  The RNAseq results of selected genes were validated by real-time RT-PCR (Fig. 4). For example, the expression of CFI, CXCR6, LGALS15, and MMP1 was significantly upregulated while TFF2 mRNA level was significantly repressed by infection only in the resistant breed (CHB), in a good agreement with the RNAseq analysis. A strong correlation in log 2 transformed fold values between the two platforms, qPCR and RNAseq, was evident (a correlation coefficient R = 0.946; Fig. 5).
Gene Ontology (GO) implicated in host resistance. Among 477 and 16 GO terms significantly enriched in CHB and CS at a P value cutoff 1.0 × 10 −4 , respectively, five were significantly enriched in both breeds (Table 4). Select GO terms that may be implicated in the development of host resistance to Haemonchus infection are listed in Table 5. Several GO related to complement activation (both classical and alternative pathways) and its regulation were significantly enriched only in CHB. Numerous cell cycle related GO were significantly enriched as well (Fig. 6). GO related to secretory granule and gastric acid secretion were also enriched, suggesting that the ability to regulate secretory and gastric function of the host may be involved in the development of host resistance. Furthermore, the regulation of inflammation at the site of infection (mucosa), including arachidonic acid metabolism, cyclooxygenase pathway, and positive regulation of MAPK cascade, as well as leukocyte migration were also implicated in host resistance. On the other hand, four of the 11 GO unique to CS were related to muscle contraction.

Discussion
Parasite resistance refers to the ability of the host to avert infection, resulting in reduced worm burden 17 . Numerous factors affect this trait. Among them, host genetics play a predominant role in controlling the development of resistance while host sex, age, and prior exposure are also important 18 . Differences in parasite resistance and susceptibility existing in various sheep breeds have been long recognized 3 . Moreover, inter-and intra-host variations in resistance are evident in certain sheep populations 18 . Identifying genetics components controlling inter-, and intra-breed differences in parasite resistance has both pragmatic and theoretical implications. Towards   this end, numerous efforts have been made over the decades to unravel genes and/or genetic variants responsible for resistance, partially driven by strong desires to breed farm animals with strong resistance traits. Traditional QTL analysis and Genome-wide Association Studies (GWAS) have led to reports of dozens of QTL or markers on almost every ovine chromosome that are associated with various resistance phenotypes, such as fecal egg counts, packed cell volume, and parasite-specific antibody titers 9,11,13 . Nevertheless, the development of parasite resistance relies upon the precise control of expression of the host genome. Understanding these regulatory elements will be crucial towards unraveling their functional relevance. As a result, while much progress has been made to identify genes associated with nematode resistance in sheep during the past few years 19,20 , an in-depth comparison and characterization of transcriptome responses of various breeds and populations, especially those local indigenous breeds harboring varying degrees of parasite resistance and susceptibility, is urgently needed. The two indigenous breeds of sheep native to the Canary Islands, CHB and CS, display unique and distinct differences in parasite resistance and susceptibility. When co-grazing together on the same pasture under natural infections, differences in fecal trichostrongylid egg counts between CHB and CS are consistently observed 14 . Under experimental infections with H. contortus, CHB has a significantly lower, by approximately 50%, worm burden than CS, a undeniable trait of parasite resistance 14,15 , which is confirmed in this study. Moreover, worms recovered in CHB tend to have significantly shorter body length than those in CS. A significantly lower EPG value is consistently observed in the feces of CHB sheep than those of CS animals during experimental infection. For example, at 27dpi, the mean EPG in CS is 5 fold higher than in CHB 14 . CHB sheep not only shed significantly fewer parasite eggs but also tend to have a delayed egg production, indicating an anti-fecundity effect of the immune response in this breed. The results from this study show that at 20 dpi, no parasite eggs were recovered in the feces of infected CHB animals while EPG in the feces of infected CS sheep reached 262.50 (± 287.54, SD). This observation is in agreement with the previous findings 14 . Haemonchus contortus infection generally elicits a potent Th2 immune response in small ruminants. A strong upregulation of several well-known Th2 cytokines by infection in CHB were observed in this study. Previous studies in the Canary Island breeds suggest that divergence in immune response mechanisms exist between CHB and CS. Among various immune cells, abomasal eosinophil numbers are 2 fold higher in CHB than in CS, suggesting that CHB sheep may have developed abilities for enhanced recruitment of eosinophils to the site of infection (abomasal mucosa). Furthermore, CHB sheep have evolved mechanisms attacking the adult stage of the Haemonchus parasite, especially its reproduction, as evidenced by the fact that fecundity is negatively correlated with eosinophils and γ δ T cells in the abomasal mucosa 15 . However, the precise molecular mechanisms of the parasite resistance in CHB breed remain largely unclear.
In this study, we identified a total of 477 and 16 Gene Ontology (GO) terms that are significantly enriched in the transcriptome of resistant and susceptible sheep breeds in responses to Haemonchus infection, respectively.

Table 4. Gene Ontology (GO) biological processes (BP) significantly enriched in both resistant (CHB) and susceptible (CS) breeds.
Scientific RepoRts | 6:26200 | DOI: 10.1038/srep26200 Among them, only five enriched GO were shared by both breeds. These GO, including leukotriene metabolic process, eicosanoid biosynthesis process, adaptive immune response, and unsaturated fatty acid biosynthesis, likely represents the basic mechanisms of host immune responses to helminth infection in sheep. Indeed, local inflammatory responses are known to be involved in the development of host resistance 21 . The enriched GO unique to the susceptible CS breed were predominantly muscle contraction-related. In cattle, our previous results suggest that smooth muscle hypercontractility induced by primary infection of the intestinal worm Cooperia oncophora represents an important aspect of host responses 22 , as in several other host-parasite systems 23 . In rodent models, helminth infection results in an increase in thickness of jejunal smooth muscle layers. Other studies also support the idea that enhanced muscle contractility appears to be associated with more rapid worm expulsion and stronger host immune responses 24 . In addition, granzyme-mediated apoptotic signaling pathway (GO:0008626) may play an important role in protecting the host from H. contortus infection in the susceptible CS breed. Complement activation as one of the earliest events in host immune responses to helminth infection plays an important role in the development of host resistance 25 . At least 11 complement related genes, such as CFI and C7, were significantly impacted by infection in the resistant CHB breed compared to uninfected controls while none of these genes were affected by infection in the susceptible breed. As a result, both classical and alternative   complement pathways appeared to be activated in the resistant breed. Furthermore, two GO molecular functions related to C5a (GO:0031714) and C5L2 anaphylatoxin chemotactic receptor binding (GO:0031715) were significantly enriched in the resistant breed. It is conceivable that these peptides play a critical role in subsequent recruitment of effector cells, such eosinophils and mast cells, to the site of infection.
Intriguingly, approximately 15% of the 711 genes whose transcript abundance were significantly altered by infection in the resistant breed were cell cycle related. The vast majority of these genes were significantly upregulated (Table 3). These genes included several cyclins, minichromosome maintenance complex components, and various kinesin family members (Table 3). In addition, a large class of genes significantly impacted in the resistant breed was ECM related (Table 2). ECM related genes are required during the classical stages of wound repair, inflammation, new tissue formation, and remodeling 26 . Previous studies identified essential roles of Th2 cytokines in limiting tissue damage during helminth infection in rodent models, especially the involvement of IL17 in the early stage of tissue repair via its role in neutrophil recruitment 27 . In this study, the transcript abundance of IL17RB was increased approximately14 fold by Haemonchus infection only in the resistant breed. Of note, upregulation of Th2 cytokines IL10 and IL13 by infection was more profound in CHB than in CS. Together, our findings suggest that the accelerated tissue repair ability, likely mediated by Th2 cytokines, has evolved in the resistant CHB breed. Recently, a significant SNP marker on ovine chromosome (OAR) 6 was reported to affect one of key resistant phenotypes in sheep, EPG 9 . This maker explains approximately 4% of the variance observed for EPG. It is suggested that there may exist up to 3 QTL within the 5 Mb region of this locus (73. , in addition to a fourth QTL at 55.9-62.6 Mb on OAV6. Several earlier reports also indicate the presence of the QTL lined to EPG in various sheep breeds 28,29 . Among 21 differentially expressed genes located on OAV6 identified in our study in CHB breeds, at least 5 genes are located within 15 Mb of this marker. The expression of four genes were significantly induced whereas one, albumin, was repressed by infection. Of note, mast cell stem cell factor (SCF) receptor KIT gene (chr6:70189728:70234612) is the closest to the SNP marker. Two major receptors, KIT and the high affinity receptor for IgE, are responsible for regulating various mast cell functions, including chemotaxis, proliferation, apoptosis, and cytokine releases 30 . The critical roles of mast cells in host immune responses to helminth infection have long been recognized 31 . Neutralization of KIT and its ligand, SCF, using monoclonal antibodies completely abrogates the mast cell hyperplasia generated by T. spiralis infection in mice, resulting in drastically delayed worm expulsion and a reduced mucosal eosinophilia 32 . This finding suggests that KIT plays an important in host-parasite interaction. In the past few years, increasing evidence suggests that the epidermal growth factor like molecule, amphiregulin (AREG), plays critical roles in regulating immunity and inflammation as well as in enhancing host resistance to helminth parasites 33,34 . In rodent models, T. suis infection increases AREG expression, in parallel with the expression of Th2 cytokines IL4 and IL13 33 . Furthermore, worm clearance is significantly delayed at 14 dpi in AREG deficient mice, which correlates with reduced proliferation of colonic epithelial cells. Recent studies show that AREG is critical for efficient regulatory T cell function 35 and may play an important role in orchestrating immunity, inflammation, and tissue repair 34 . In this study, AREG transcript abundance was significantly enhanced by infection only in the resistant breed, suggesting that this gene may play an important role in the development of host resistance. It would be intriguing to identify SNPs in both coding and promoter regions of the genes located within or closer to the QTL related to parasite resistance on OAV6, including AREG, and correlate the observed genetic variations with various resistant phenotypes. Moreover, dissecting mechanisms of transcriptional regulation of AREG and understanding how it promotes epithelial cell proliferation and regulates host immunity in the gastrointestinal tract warrant further investigation.
In conclusions, the two sheep breeds native to the Canaria Island displayed a distinct difference in several Haemonchus contortus resistant phenotypes under both natural and experimental infections. CHB tends to have significantly reduced worm burden, delayed egg production, and decreased fecal egg yield (counts) than the susceptible Canaria Sheep. A broad range of mechanisms have evolved in resistant CHB to provide protection against H. contortus. Readily inducible acute inflammation responses, complement activation, accelerated cell proliferation and subsequent tissue repair, and innate and acquired immunity directly against worm fecundity are likely to contribute to the development of host resistance to gastrointestinal nematode infection in the CHB breed.

Methods
Animals and parasitology. Male lambs of CHB (11 animals) and CS breeds (12 animals) were obtained from local farms in the Gran Canaria Island (Spain), weaned, and kept in pens at the Faculty of Veterinary Science, University of Las Palmas de Gran Canaria until they were approximately one year old. The animals were fed with a commercial pelleted sheep ration ad libitum and had free access to water throughout the experimental period. The animals were drenched upon arrival with levamisole (Cyber, Fort Dodge, Spain) at the recommended dose (1 ml/10 kg bodyweight) and remained free of parasites (as determined by fecal egg counts) until experimental parasite inoculation. Seven CHB and eight CS animals were inoculated intraruminally with 20,000 H. contortus infective L3 larvae. Four age-matched animals of each breed remained uninfected and served as controls. The experimental infection was allowed to progress for 20 dpi. The time point chosen for this study was based on the results from a previous report that the difference in resistance phenotypes, especially mean EPG values, is most profound between the two breeds 14 . Furthermore, no fecal parasite eggs become detectable at this time point in CHB. At 20 dpi, both infected and control animals were sacrificed. EPG values were monitored twice during the experiment, prior to the experimental inoculation and immediately prior to necropsy using the MacMaster technique. Adult worms as well as immature larvae from both contents and the tissue of the abomasum were isolated and counted. The fundic abomasum tissue was then sampled and snap frozen in liquid nitrogen prior to storage at − 80 °C until total RNA was extracted. The Haemonchus strain used in this trial was initially donated by Drs. Knox and Bartley (Moredun Research Institute, Edinburgh, Scotland) and passaged through successive inoculations in sheep at the premises of the Faculty of Veterinary Science, University of Las Palmas de Gran Canaria (Spain). During the experiment, all animal protocols were approved by the Animal Care and Use Committee of University of Las Palmas per the Institutional Animal Care and Use Committee (IACUC) guidelines. All experimental procedures were carried out in accordance with the approved protocols.
RNA extraction and sequencing using RNA-seq. Total RNA from fundic abomasal samples of both CHS and CS sheep breeds was extracted using Trizol (Invitrogen, Carlsbad, CA, USA) followed by DNase digestion and Qiagen RNeasy column purification (Qiagen, Valencia, CA, USA), as previously described 22,36 . The RNA integrity was verified using an Agilent BioAnalyzer 2100 (Agilent, Palo Alto, CA, USA). High-quality RNA (RNA integrity number or RIN > 7.5) was processed using an Illumina TruSeq RNA sample prep kit following the manufacturer's instructions (Illumina, San Diego, CA, USA). Pooled RNAseq libraries were sequenced at 2 × 101 bp/sequence read using an Illumina HiSeq 2000 sequencer, as described previously 37 . Approximately 56 million paired-end sequence reads per sample (mean ± SD = 55,945,621 ± 41,305,493.24; N = 23) were generated. The metadata and raw sequences files related to this project were deposited in the NCBI Sequence Read Archive (Accession #SRP059627).
Scientific RepoRts | 6:26200 | DOI: 10.1038/srep26200 Data analysis and bioinformatics. Raw sequence reads were first checked using FastQC (http://www.bioinformatics.bbsrc.ac.uk/projects/fastqc/). The effect of trimming of low-quality nucleotides on genome alignment was examined using STAR algorithm 38 . Raw sequence reads (FASTQ files) of 23 samples were mapped against the ovine reference genome Oar_v3.1 using STAR (v2.3.1t) with default parameters. The uniquely mapped read were used to count against the Ensembl annotation Oar_v3.1 for calculating the number of reads per gene. The counts of all samples were tabulated. This table was then inputted to DESeq 39 for normalization and identification of differentially expressed genes between infection and control groups of both CHS and CS using the standard workflow as described 38 . To correct for multiple hypothesis testing, the Benjamini-Hochberg procedure 40 was used with an FDR cutoff of 0.05. Gene Ontology (GO) analysis over differentially expressed genes was performed using Fisher's exact test.
Real-Time RT-PCR (qPCR) validation. In order to validate the results obtained in the RNAseq analysis, the expression of 7 genes (see Supplementary table for their primer sequences) was determined by qPCR as previously described 36 . Ovine ribosomal protein L19 gene (RPL19), whose expression remained stable among the experimental samples, was used as an endogenous reference gene for all reactions. cDNA was synthesized from high quality total RNA (RIN > 7.5) using Superscript II reverse transcriptase (Invitrogen, Carlsbad, Carlsbad, CA) according to the manufacturer's instructions. All qPCR reactions were carried out in 96-well plates in a 7500 Real-Time PCR System and analysed with a 7500 Software v2.0.6 (Applied Biosystems, NY). Samples were run in duplicate in a total volume of 25 μ l containing the following: 5 μ l of cDNA (100 ng), 1 μ l of primer mix (forward and reverse, 10 nM each), 0.1 μ l of ROX, 12.5 μ l of SYBR GreenER qPCR SuperMix (Invitrogen) and 6.4 μ l of dd water. The amplification reactions were subjected to a holding stage of 50 °C for 2 min, followed by an initial denaturation at 95 °C for 10 min. The reactions were then followed by 40 cycles of 95 °C for 30 sec, 60 °C for 30 sec and 72 °C for 32 sec. Melting curves were obtained from 60 °C to 95 °C. Relative gene expression values were determined using a standard curve method. Briefly, eight 10-fold serial dilutions of a pool of cDNA samples were used to generate standard curves for each gene to calculate relative gene expression levels. These results were then normalized to RPL19 gene expression levels for each sample.