# Dynamic transcriptomic changes of goat abomasal mucosa in response to Haemonchus contortus infection

Hadeer M. Aboshady, Nathalie Mandonnet, Yoann Félicité, Julien Hira, Aurore Fourcot, Claude Barbier, Anna M. Johansson, Elisabeth Jonas, Jean-Christophe Bambou

<figure><img src="image_1.png" /></figure>

## To cite this version:

Hadeer M. Aboshady, Nathalie Mandonnet, Yoann Félicité, Julien Hira, Aurore Fourcot, et al. Dynamic transcriptomic changes of goat abomasal mucosa in response to Haemonchus contortus infection. Veterinary Research, 2020, 51 (1), pp.1-12. 10.1186/s13567-020-00768-y. http://hal-02514818

HAL Id: hal-02514818

https://hal.science/hal-02514818

Submitted on 23 Mar 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

# Dynamic transcriptomic changes of goat abomasal mucosa in response to *Haemonchus contortus* infection

Hadeer M. Aboshady<sup>1,2,3,4</sup>, Nathalie Mandonnet<sup>3</sup>, Yoann Félicité<sup>3</sup>, Julien Hira<sup>3</sup>, Aurore Fourcot<sup>3</sup>, Claude Barbier<sup>5</sup>, Anna M. Johansson<sup>2</sup>, Elisabeth Jonas<sup>2</sup> and Jean-Christophe Bambou<sup>3*</sup>

## Abstract

Gastrointestinal nematode (GIN) infections are one of the major constraints for grazing sheep and goat production worldwide. Genetic selection for resistant animals is a promising control strategy. Whole-transcriptome analysis via RNA-sequencing (RNA-seq) provides knowledge of the mechanisms responsible for complex traits such as resistance to GIN infections. In this study, we used RNA-seq to monitor the dynamics of the response of the abomasal mucosa of Creole goat kids infected with *Haemonchus contortus* by comparing resistant and susceptible genotypes. A total of 8 cannulated kids, 4 susceptible and 4 resistant to GIN, were infected twice with 10 000 L3 *H. contortus*. During the second infection, abomasal mucosal biopsies were collected at 0, 8, 15 and 35 days post-infection (dpi) from all kids for RNA-seq analysis. The resistant animals showed early activation of biological processes related to the immune response. The top 20 canonical pathways of differentially expressed genes for different comparison showed activation of the immune response through many relevant pathways including the Th1 response. Interestingly, our results showed a simultaneous time series activation of Th2 related genes in resistant compared to susceptible kids.

## Introduction

Gastrointestinal nematodes (GIN) are an important constraint on grazing ruminants worldwide. These parasites can cause mortality especially in small ruminants but their main effect is reduced productivity [1, 2]. Anthelmintic treatments are the mainstay of current treatment but are threatened by the evolution of drug resistance in parasite populations [3]. Besides, the environmental side-effect of anthelmintic residues is no longer desirable for sustainable production and the increased demand for chemical-free animal products. Therefore, there is a need for additional control strategies. The introduction of resistance to GIN traits in small ruminants breeding

schemes, would be a promising sustainable method to control GIN infection [1, 4, 5].

Resistance against most of the common diseases are complex traits involving many genes, the detection of causative variations is therefore a complex task. Currently selection against GIN relies on indirect measures such as fecal egg count (FEC), packed cell volume (PCV) and blood eosinophilia [6–9]. A major disadvantage of these methods is that animals must be infected either naturally or experimentally for these measures. An alternative is the identification and the selection of genes that are responsible for resistance to GIN infection. Several studies investigated the molecular and cellular processes associated with GIN resistance in different tissues such as duodenal [10–12] and abomasal mucosa [13, 14] and draining lymph nodes [15–19] mainly in sheep. However, only a few studies have investigated the biological processes associated with GIN resistance in goats.

<footnote>*Correspondence: jean-christophe.bambou@inrae.fr</footnote>

<footnote n="3"><sup>3</sup> URZ Recherches Zootechniques, INRAE, 97170 Petit-Bourg, Guadeloupe, France
Full list of author information is available at the end of the article</footnote>
It has been shown that whole-transcriptome analysis via RNA-seq is a key tool to investigate the molecular mechanisms responsible for complex quantitative traits such as resistance to GIN infection [20]. A detailed understanding of the genes and biological mechanisms involved in resistance and protective immunity would provide new phenotypic and genetic markers for effective breeding schemes [21].

Previously, we investigated the transcriptome variation in response to GIN infection in goats at 42 days post-infection (dpi) [22]. The results indicated that the maintenance of the integrity of the mucosa was probably the priority for the host at this late infection stage (42 dpi). The present study aimed to identify the changes over time in the molecular pathways and immunity development in response of Creole goats to GIN infection by analyzing the transcriptome of abomasal mucosa of resistant and susceptible kids at different time points post-infection.

# Materials and methods

## Ethics statement

All animal care handling techniques and procedures as well as the procedures for experimental infection, tissue sampling and slaughtering were approved by the French Ethic Committee n°069 (Comité d'Ethique en Matière d'Expérimentation Animale des Antilles et de la Guyane, CEMEAAG) authorized by the French Ministry of Higher Education, Research and Innovation. The experiment was performed at the INRA Experimental Facilities PTEA (Plateforme Tropicale d'Expérimentation sur l'Animal), in Guadeloupe (French West Indies) (16° 20' latitude North, 61° 30' longitude West), according to the certificate number A 971-18-02 of authorization to experiment on living animals issued by the French Ministry of Agriculture.

## Animals and experimental design

The experimental design is described in Figure 1. Nine month-old Creole kids were chosen from the experimental flock of PTEA (Plateforme Tropicale d'Expérimentation sur l'Animal) in which the estimated

<figure><img src="image_5.png" /><figcaption>**Figure 1** Experimental design. The Creole kids reared at pasture were chosen from the experimental flock at 9 month-old. The animals were drenched and housed indoors under worm-free conditions in a single pen. Kids were orally infected with a single dose of 10 000 *Haemonchus contortus* third-stage larvae (L3) (challenge 1 and 2, lasted for 5 weeks each).</figcaption></figure>
breeding value (EBV) was calculated regularly since 1995 for each animal for FEC at 11 months of age following natural mixed infection on pasture taking into account the FEC of its ascendants and pedigree. Before the experiment the kids were reared at pasture with a limited level of GIN contamination (FEC < 500). The FEC of the 8 kids (n=4 resistant and n=4 susceptible), chosen on the basis of their extreme EBV in their cohort, were not statistically different. The EBV was estimated by taking into account the FEC of their ascendants and their pedigree. The averages EBV of the 2 groups were distant by 1.04 genetic standard deviation. The animals were drenched with moxidectine (Cydectine®, Fort Dodge Veterinaria S.A., Tours, France, 300 µg/kg) and housed indoors under worm-free conditions in a single pen, 1 month before the start of the experiment. Kids were orally infected with a single dose of 10 000 *Haemonchus contortus* third-stage larvae (L3) in two consecutive challenges. Each challenge lasted for 5 weeks with 8 weeks interval between the end of challenge 1 and the start of challenge 2. At the end of the challenge 1, the kids were drenched with moxidectin (Cydectine®, Fort Dodge Veterinaria S.A., Tours, France, 300 µg/kg) and 4 weeks later a fistula was surgically implanted in the abomasum of each animal to allow abomasal mucosa sampling at 0, 8, 15 and 35 dpi. After another period of 4 weeks, the animals were orally infected with a single dose of 10 000 *H. contortus* L3 (challenge 2). For FEC measurements during the experimental infection, approximately 10 g of faeces were collected in plastic tubes directly from the rectum of each animal, and transported from the experimental facility to the laboratory in refrigerated vials. The samples were individually analysed using a modified McMaster method for rapid determination and FEC was expressed as the number of eggs/g faeces [9].

## Surgical procedure

The custom designed abomasal cannula consisted of a flexible plastic tube with a length of 7 cm and a diameter of 2 cm with a rounded base of 4 cm in diameter. This flexible plastic was chosen to limit the possibility of mechanical abrasion of the mucosal surface of the abomasum. The animals were fasted 16 h before cannula insertion surgery. The animals were premedicated with ketamine (2 mg/kg IV, Le Vet Pharma, Wilgenweg, Netherlands), xylazine (0.2 mg/kg IM, Le Vet Pharma, Wilgenweg, Netherlands) and oxytetracycline (20 mg/kg IM, Eurovet Animal Health, Handelsweg, Netherlands). The animals were positioned in left lateral recumbency. Skin over the surgical site was shaved and prepped with povidone iodine (Vétédine, Laboratoire Vetoquinol S.A., Lure, France). A ventral midline incision was made to locate and externalise the abomasum. A 3 cm purse-string

suture (Silk 2-0) was placed midway between the lesser and greater curvature and a stab incision was made in the center to insert the cannula. Then, the purse-string suture was tightened and tied off. To maintain the abomasum in an anatomically correct position, another stab incision was made in the abdominal wall at 10 cm from the laparotomy incision on the right paramedian area to enable the cannula to be passed freely through. An external flange was placed over the external part of cannula and fixed with adhesive fabric plaster strip. A sterile compress was inserted into the cannula as stopper. After the surgical procedure, all the animals were housed individually with free access to fresh water and hay.

## Biopsy sampling procedure

Biopsy specimens were taken from the abomasal mucosa using a flexible endoscope (FG-24 V, Pentax, France). The biopsies samples of 2 × 2 × 2 mm taken with the endoscopic forceps with window (model KW1815S) were quickly snap frozen into liquid nitrogen and stored at -80 °C until RNA extraction. The animals were restrained in a harness made with a surgical drape allowing animal legs to protrude and which exposed the cannula. No sedation was used since no signs of discomfort or pain were observed during or after the procedure. The sterile compress inserted into the cannula was removed and the abomasal contents collected. The endoscope was introduced into the abomasal lumen and 3 biopsies per animal and per time points were taken from the abomasal folds of the fundic mucosa. At each time point the whole fundic mucosa was observed and no sign of mucosal injury due to the previous sampling was observed.

## RNA extraction and sequencing

Total RNA was extracted using the NucleoSpin® RNA isolation kit (Macherey-Nagel, Hoerdt, France) following the manufacturer's instructions, except that DNase digestion was performed with twice the indicated amount of enzyme. The total RNA concentration was measured with NanoDrop 2000 (ThermoScientific TM, France). The RNA integrity was verified using an Agilent Bioanalyzer 2100 (Agilent Technologies, France) with a RNA Integrity Number of >7.5. The extracted total RNA was stored at -80 °C until sequencing.

High-quality RNA from all samples was processed for the preparation of cDNA libraries using an Illumina TruSeq RNA sample prep kit for mRNA analysis following the Illumina's protocols. After quality control and quantification, cDNA libraries were pooled in groups of 6 and sequenced on 5 lanes on the HiSeqTM 2000 (Illumina® NEB, USA) to obtain approximatively 30 million reads (100 bp paired-end) for each sample with insert sizes ranging from 200 to 400 base pairs.
## Bioinformatics and data analysis

The quality control check on raw reads in FASTQ format were processed using FASTQC and the Q20, Q30 and GC contents of the clean data were calculated. The Salmon software (version 0.9.1) was used for transcript quantification [23]. NCBI RefSeq reference transcript of the *Capra hircus* genome (assembly ARS1) was used to build the index within Salmon. The reads from each sample were mapped to the same index and quantified. Unix commands were used to obtain corresponding gene and transcript identifiers from the NCBI RefSeq annotation of the *Capra hircus* (ARS1). Using these identifiers, the tximport (version 1.8.0) package was used to import data into the R software (v3.5.1) and summarize the TPM estimates obtained from the Salmon tool of all samples at the gene level [24]. This process produced a global count file on which the statistical analyses were performed. A threshold of greater than or equal to 5 counts across samples was applied in order to remove genes showing low expression.

Partial least squares discriminant analysis (PLSDA) had been conducted using the mixomics package within R [25]. In this analysis, x was the matrix of gene expression values (count table) and the classes of y were given as resistant and susceptible. Each row of the x matrix represented the gene expression values for a sample, and each column corresponded to a gene.

Differentially expressed genes (DEG) of read counts were identified using the Bioconductor package DESeq2 within R [26]. Ten comparisons were performed; three comparing day 0 with day 8, 15 or 35 post-infection in the susceptible group, another three comparing the same days in the resistant group and four comparing samples from resistant versus susceptible animals at day 0, 8, 15 and 35 post-infection. To account for multiple testing, genes were filtered using a Benjamini and Hochberg false discovery rate (FDR) of <0.001. Final DEG were determined on the basis of their fold change values to be $log_2 \ge 1.0$ for up-regulated genes and $\le -1.0$ for down-regulated genes. Gene ontology (GO) analysis for the biological processes was performed to identify the biological function classification of the genes, which describes properties of genes and their products. DEG are functionally grouped into the biological processes looking for significantly enriched functions compared to the human genomic background due to the lack of goat (*C. hircus*) GO data. GO enrichment analysis and GO annotations plotting were performed using the clusterProfiler R package [27]. All enriched GO terms that possessed a *p*-value <0.01 were displayed and the top 5 biological processes for each comparison were plotted. Analysis of canonical pathways and regulator effects were performed using Ingenuity pathway analysis (IPA)

software (Ingenuity Systems, Redwood City, CA, USA) for DEG in each comparison.

Faecal egg counts (FEC) were measured twice a week after infection from 21 to 36 dpi. The FEC variance was normalized using log transformation. PROC MIXED procedure (v. 9.4, SAS Inst. Inc., Cary, NC, USA, 2012) was used to test statistical differences. The differences were considered significant when *p* < 0.05. The results are presented after back transformed.

## Quantitative real-time PCR (qRT-PCR) validation

To validate the results of the RNAseq analysis, the gene expression for a total of 9 genes ($n=6$ for each comparison: resistant vs susceptible at 0, 15 and 35 dpi, and resistant and susceptible for 0 vs 8 dpi, 0 vs 15 dpi and 0 vs 35 dpi) was determined by qRT-PCR. The endogenous control for all reactions was goat *ACTB* (*actin beta*) gene whose expression remained stable among the samples. The cDNA was synthesize with a total of 2 µg of high quality total RNA (RIN > 7.5) by using M-MLV Reverse Transcriptase (Promega, Charbonières, France) according to the manufacturer's instructions. All qRT-PCR reactions were carried out in 48-well plates in a Prime Pro 48 Real-Time PCR System and analyzed with the ProStudy Software v5.2.10 (Technel, Staffordshire, UK). Taqman® predesigned gene expression assay (Table 1) and the universal PCR master mix were purchased from Applied Biosystems and the analysis were performed according to the manufacturer's instructions (ThermoFisher Scientific, Applied Biosystems, Courtaboeuf, Villebon-sur-Yvette, France). Samples were analyzed in duplicate in a total volume of 20 µL containing: 4 µL of cDNA, 10 µL of 2X TaqMan® Fast Advanced Master Mix, 1 µL of TaqMan® Gene Expression Assays 20X (ThermoFisher Scientific, Applied Biosystems, Courtaboeuf, France) and 5 µL of distilled RNAse DNAse-free water. Relative gene

<table><caption>**Table 1 List of target genes for qRT-PCR validation and assay IDs according to the manufacturer**</caption><thead><tr><td>Gene symbol</td><td>Gene description</td><td>Assay IDs</td></tr></thead><tbody><tr><td>ACTB</td><td>actin beta</td><td>Ch04810274_s1</td></tr><tr><td>CYP4F2</td><td>phyloquinone omega-hydroxylase</td><td>Ch04672252_m1</td></tr><tr><td>DUOXA2</td><td>dual oxidase maturation factor 2</td><td>Ch04786286_m1</td></tr><tr><td>CCL20</td><td>C-C motif chemokine ligand 20</td><td>Ch04791475_m1</td></tr><tr><td>IFI6</td><td>interferon alpha inducible protein 6</td><td>Ch04807049_g1</td></tr><tr><td>LST1</td><td>leukocyte specific transcript 1</td><td>Ch04741898_m1</td></tr><tr><td>NKX6-3</td><td>NK homeobox 3</td><td>Ch04677616_m1</td></tr><tr><td>OLFM4</td><td>olfactomedin 4</td><td>Ch04796577_m1</td></tr><tr><td>TFF3</td><td>trefoil factor 3</td><td>Ch04767901_m1</td></tr><tr><td>TLR4</td><td>toll-like receptor 4</td><td>Ch04654181_m1</td></tr></tbody></table>
expression values were determined using relative quantification ($2^{-\Delta\Delta Ct}$ method, [28]).

# Results

## Parasitological measures

A significant effect of the group (i.e. resistant vs susceptible), the dpi and their interaction (*p* < 0.001) was observed for FEC (Figure 2). At 21 dpi no difference was observed between groups. Thereafter the FEC was significantly lower in resistant compared to susceptible animals whatever the dpi.

## RNA sequencing and variance analysis

Alignment of RNA sequencing to the reference *Capra hircus* genome (assembly ARS1) resulted in an average of 4.5±0.1 million reads per sample. These reads correspond to 23 258 genes of the goat genome. A total of 15 188 out of the 23 258 annotated genes (65%), showed at least 5 read counts per row and were used in the subsequently analysis. The multilevel PLSDA for gene expression of infected resistant and susceptible kids explained more than 20% of the variance in its two-dimension components (Figure 3). Component 1 represented 11% of the whole variability and component 2 represented also 11% of the variation.

## Differential gene expression

The numbers of DEG for each comparison are shown in Table 2. The numbers of DEG were low for the comparison between groups (R vs S) whatever the time point. For the comparison within infected resistant or infected susceptible, the numbers of DEG were lower for 0 versus 15 dpi (678 and 1748, respectively) compared with

0 versus 8 or 0 versus 35 dpi. Meanwhile the highest number for DEG was recorded for the comparison of 0 versus 35 dpi of infected susceptible (3316) and infected resistant (2263). The fold change was on average higher when comparing different time points within each group (from −11.15 to 24.17 and from −11.83 to 9.30 for R or S respectively) than between groups at different days. Human orthologues were mapped for 72–85% of the DEG (Table 2).

## Validation of expression by qRT-PCR

qRT-PCR for nine genes was performed to validate RNA sequencing results. For the comparison of resistant

<table><caption>**Table 2** Number of differentially expressed genes (n) for the different comparisons including log₂ fold change and the number of human orthologues (including proportion of genes with human orthologues)</caption><thead><tr><th>Comparison</th><th>n</th><th>Log<sub>2</sub> fold change</th><th>Human orthologues</th></tr></thead><tbody><tr><td>Inf. R 0 vs 8 dpi</td><td>1336</td><td>-11.15, 24.17</td><td>1017 (76.12%)</td></tr><tr><td>Inf. R 0 vs 15 dpi</td><td>678</td><td>-10.81, 4.57</td><td>549 (80.97%)</td></tr><tr><td>Inf. R 0 vs 35 dpi</td><td>2263</td><td>-10.58, 6.66</td><td>1881 (83.12%)</td></tr><tr><td>Inf. S 0 vs 8 dpi</td><td>2221</td><td>-10.60, 9.30</td><td>1744 (78.52%)</td></tr><tr><td>Inf. S 0 vs 15 dpi</td><td>1748</td><td>-11.84, 8.82</td><td>1439 (82.32%)</td></tr><tr><td>Inf. S 0 vs 35 dpi</td><td>3316</td><td>-11.83, 9.23</td><td>2811 (84.77%)</td></tr><tr><td>R vs S 0 dpi</td><td>456</td><td>-7.39, 6.00</td><td>337 (73.90%)</td></tr><tr><td>R vs S 8 dpi</td><td>679</td><td>-4.27, 27.7</td><td>490 (72.16%)</td></tr><tr><td>R vs S 15 dpi</td><td>318</td><td>-5.1, 7.82</td><td>247 (77.67%)</td></tr><tr><td>R vs S 35 dpi</td><td>758</td><td>-7.34, 8.48</td><td>579 (76.39%)</td></tr></tbody></table>

Inf. R infected resistant, Inf. S infected susceptible, R vs S resistant versus susceptible, dpi days post-infection.

– Log₁₀ p-value: – Log₁₀ (p-value).

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 3** Multilevel PLS-DA of the gene expression of infected resistant and susceptible animals. Infected Resistant: Inf. R, Infected susceptible: Inf. S.</figcaption></figure>
versus susceptible animals at 0, 8, 15 and 35 dpi, the genes selected randomly among the DEG were: *DUOXA2*, *IFI6*, *CYP4F2*, *OLFM4* and *TFF3*. For the comparison of 0 versus 8, 15 and 35 dpi within the resistant and the susceptible animals the genes were respectively: *IFI6*, *CYP4F2*, *OLFM4*, *TFF3*, *TLR4* and *NKX6-3*, *CCL20*, *OLFM4*, *LST1*, *TFF3*. The log₂ fold change levels of the selected genes measured by qRT-PCR were in good agreement with the values from the sequencing data (Figure 4). The gene expression patterns from qRT-PCR were highly correlated with the sequencing results: the correlation coefficients were respectively 0.91, 0.96 and 0.81 for the comparison of resistant versus susceptible animals at different time points and the comparison of 0 versus other time points within the resistant and the susceptible animals.

## Functional classification analysis

### Gene ontology (GO)

An enriched GO term analysis for biological processes was performed using the DEG from each comparison. The top 5 significant biological processes in each term are presented in Figure 5. Comparing 0 versus 35 dpi, four out of the top 5 biological processes were the same for the resistant and the susceptible kids; meanwhile leukocyte differentiation was in the top biological process only for the resistant kids. The comparison of infected resistant at 0 versus 8 dpi showed biological processes related to the immune response within the top 5 significant processes (e.g. T cell activation, leukocyte cell-cell adhesion and lymphocyte differentiation). Positive regulation of the innate immune response was in the top 5 biological processes when comparing susceptible with resistant at 35 dpi.

### Pathway enrichment analysis

The Ingenuity Pathway Analysis was used to compare results from different comparison over time. The top canonical pathways (Figure 6) and the top upstream regulators (Figure 7) were compared. When comparing day 0 versus 35 post-infection, the top 20 canonical pathways showed a high activation of the immune response through dendritic cell maturation, *IL-8* signaling, Leukocyte extravasation signaling, *NFAT* in regulation of the immune response, *P13K* signaling in B lymphocytes, Th1 pathway and B cell receptor signaling pathways. In resistant compared with susceptible kids the B cell receptor signaling pathway was activated at 8 dpi while dendritic cell maturation and Th1 pathways were activated at 35 dpi.

<figure><img src="image_8.png" /><figcaption>**Figure 4** Fold change of differentially expressed genes measured by RNA-Seq (white) and qRT-PCR analyses (black).</figcaption></figure>

RNA-Seq: white bars, qRT-PCR: black bars. The results are presented according to the comparisons: resistant versus susceptible at 0, 8, 15 and 35 days post-infection (dpi), 0 versus 8, 15 and 35 dpi for resistant and susceptible animals respectively.

The top 10 upstream regulators of the DEG for different comparisons showed that some genes like *TGF-β1*, *TNF-α*, *IFN-γ*, *IL1-β* and *IL-6* were in the group of the top significant upstream regulators in both infected resistant and susceptible kids specially when comparing 0 versus 35 dpi. These genes were still significantly differently expressed between resistant and susceptible kids at 35 dpi (Figure 7). The *TGF-β1* gene was the top significant upstream regulator that was differently expressed in resistant compared with susceptible kids in the abomasal mucosa.
<figure><img src="image_9.png" /></figure>

## Differential of CD4+ T cell

Genes related to the CD4+ T cell activation and the fold change comparing resistant versus susceptible kids at 0, 8, 15 and 35 dpi are presented in Figure 8. The CD4+ T cell differentiation pathway showed a significant difference and a positive fold change for the majority of genes controlling the Th1 pathway when comparing resistant versus susceptible kids at 35 dpi. The expression of genes controlling the Th2 pathway showed time series activation in resistant compared with susceptible kids at different dpi: *IL2RG* activated at 8 dpi, *IL4R* and *STAT6* at 15 dpi, *GATA3* and *CCR4* at 35 dpi. Meanwhile the expression of *IL4R* and *STAT6* at 35 dpi is higher in susceptible kids. The expression levels for genes controlling the Th17 pathway showed a positive fold change for *STAT3* and *RORC* in resistant kids at 15 dpi, then for *IL17F* at 35 dpi while for *STAT3* the expression was higher in susceptible kids at 35 dpi.

Comparing resistant versus susceptible at 0 dpi (before the experimental infection), the expression of *IL17F* was three times higher in resistant kids. No difference of *FOXP-3* expression was observed between resistant and susceptible whatever the dpi, while the expression of *TGF-β1* was significantly higher in resistant kids at 8 dpi and lower at 35 dpi.

## Discussion

This study aimed to investigate the kinetic changes in mucosal molecular pathways and immunity development of resistant and susceptible Creole kid goats in response to *H. contortus*. The classification of the animals as resistant or susceptible was explained at 22% by the gene expression profile. *H. contortus* infection induced a high number of DEG in the mucosa of both resistant and susceptible animals whatever the time points while the numbers of DEG were much lower when comparing resistant versus susceptible animals at the different time points of
<figure><img src="image_10.png" /></figure>

infection. This result indicates that most genes involve in the host response against *H. contortus* infection were similar in susceptible and resistant animals.

GO of enriched biological processes showed an earlier activation of immune biological processes in resistant kids. Indeed, the top biological processes at 8 dpi were T cell activation, leukocyte cell-cell adhesion and lymphocyte differentiation. One of the top significant pathways was B cell receptor signaling. In keeping with this results, McRae et al. reported an early immune response to *Tela-dorsagia circumcincta* in resistant sheep at 7 dpi [19]. The same top four biological processes were observed in resistant and susceptible animals when comparing 0 and 35 dpi. However, none of these processes appeared in these top biological processes when comparing susceptible with resistant animals at 35 dpi, suggesting that at 35 dpi the host priority at the abomasal mucosa interface would be similar for resistant and susceptible kids.

The Th1 pathway was one of the top pathways identified in most of the comparison performed in this study. Upstream regulators of the genes involved in the Th1

processes include *TNF-α* and *IFN-γ*, which were also identified as DEG. In accordance with this result, a transient increase of the expression of *TNF-α* and *IFN-γ* was observed earlier after *H. contortus* infection in sheep both in the abomasal mucosa and the draining lymph nodes [29–31]. However, a non-protective Th1 response associated with an increased expression of cytokines, as *TNF-α* and *IFN-γ*, was observed respectively in susceptible and primary infected sheep infected with *H. contortus* [32, 33]. Indeed, studies on murine models demonstrated for a long time that the protective response against GIN parasites is better associated with the Th2 polarization of the immune response [34], while host susceptibility is associated with a Th1 response [35, 36]. In ruminants, the Th1/Th2 dichotomy remains controversial despite studies showing a correlation between host resistance and a polarized Th2 immune response [37–39]. A simultaneous increased expression of Th1- and Th2-type cytokines was shown in cattle infected with *Ostertagia ostertagi* [40–42]. Similarly, looking at differential of CD4+ T cell, we found signals for Th1 and Th2 activation
<figure><img src="image_11.png" /><figcaption>**Figure 7** Top 10 upstream regulators of differentially expressed genes for infected resistant and infected susceptible animals. Infected Resistant: Inf. R, Infected susceptible: Inf. S. Comparison of day 0 with 8, 15 and 35 dpi and resistant versus susceptible animals (R vs S) at 8, 15 and 35 dpi. The color gradient moves from white [no significant difference, -Log₁₀ (p-value) = 0] to red [significant difference, -Log₁₀ (p-value) ≥ 2].</figcaption></figure>

at 35 dpi in resistant animals when comparing them with susceptible animals. Caucheteux et al. [43] reported that the expression of IL1-β gives rise to inflammatory Th2 cells that are specialized to induce allergic inflammatory responses, whereas Th2 primed in the absence of IL1-β are more important as regulatory cells, that is amplifiers of Th2 cells and antibody response by B cells. Our results showed IL1-β in the top upstream regulator genes controlling infection response.

Transforming growth factor beta (TGF-β) is a multi-functional cytokine known for its regulatory activity and the induction of peripheral tolerance [44]. We found that the gene expression profile of TGF-β1 was the top significant upstream regulator when comparing the dynamics of infection in resistant and susceptible animals. TGF-β1 was activated in susceptible and inhibited in resistant animals at 35 dpi. The same was previously reported in other studies in goats [22, 45] and also a study on sheep infected with *H. contortus* [15]. The underlying mechanisms could be a manipulation of the host immune response by *H. contortus*, notably through the induction of the secretion of *IL-10* and *TGF-β1* by goat monocytes to promote an anti-inflammatory environment favorable for worm survival [46]. This hypothesis needs to be investigated.

A gene expression profiling study of the abomasal mucosal and lymph nodes of resistant and susceptible

goats in response to *H. contortus* infection at 42 dpi has previously reported that the maintenance of the integrity of the mucosal barrier is one of the priorities of the host response at the late stage of infection [22]. The study presented here studied the dynamics of the gene expression in the goat abomasal mucosa in response to *H. contortus* infection using information from the whole transcriptome of resistant and susceptible kids. A time series activation of Th2 genes was identified for resistant animals compared with the susceptible ones. The later activation of some genes in susceptible animals indicated that the Th2 response was activated earlier in resistant kids compared to susceptible kids. Transcriptional profiling of the abomasal lymph node from Scottish Blackface lambs showed that resistant animals are generating an earlier immune response to *T. circumcincta* infection compared to susceptible animals [19]. This difference was through pathways relating to the inflammatory response, migration of T lymphocytes and synthesis of reactive oxygen species [19].

IL17 is the leading inflammatory cytokine in the Th17 cell populations [47]. Neither the IL17A nor the IL17F genes have been described in studies analyzing the resistance to GIN in sheep. Nonetheless, IL17 transcripts have been shown to be upregulated in the bovine abomasal mucosa after 24 days of single *O. ostertagi* challenge and 60 days of trickle experimental or natural infection
<figure><img src="image_12.png" /></figure>

[42]. However, the positions of these interleukin genes have been found to be relatively close to the DRB1 gene in sheep [48], which has been reported to be associated with GIN resistance in sheep [49, 50]. In the present study, *IL17F* was the gene showing the most significant expression difference at day 0 of infection, having an expression three times higher in resistant compared with susceptible kids. Future experiments should investigate the potential of this gene as a pertinent biomarker in a selection program.

The present study showed that *H. contortus* infection in goat induces a marked immune response at the mucosal level in resistant animals, which is characterized by the simultaneous upregulation of Th1 and Th2 genes. Our results suggested differences in the time

series activation for Th2 genes, indicating that the immune response is activated earlier in resistant kid goats compared to the susceptible ones. We also found that TGF-β1 has a major regulator role during GIN infection in goats.

### Acknowledgements

The authors want to give thanks to the Duclos team for care and handling of the animals: F. Pommier, F. Nimirf, F. Periacarpin, C. Deloumeaux and M. Jean-Bart. A special thanks to Dr Harry Archimede for his valuable contribution for the surgical procedure.

### Authors' contributions

Conceptualization: NM and JCB. Formal Analysis: HMA and JCB. Funding acquisition: NM, EJ, JCB. Investigation: YF, JH, AF, CB. Methodology: HMA, JCB. Supervision: AMJ, EJ, JCB. Writing (original draft): HMA. Writing (review and editing): AMJ, EJ, JCB. All authors read and approved the final manuscript.
## Funding

This study was funded by the Project MALIN (La Région Guadeloupe and Fonds Européens FEDER). H.M.A was supported by a doctoral fellowship from the project European Graduate School in Animal Breeding and Genetics.

## Availability of data and materials

All data generated during this study are available in the NCBI SRA repository. All other relevant data are included in this published article.

## Competing interests

The authors declare that they have no competing interests.

## Author details

¹ AgroParisTech, Paris, France. ² Department of Animal Breeding and Genetics, Swedish University of Agriculture Science, Uppsala, Sweden. ³ URZ Recherches Zootechniques, INRAE, 97170 Petit-Bourg, Guadeloupe, France. ⁴ Department of Animal Production, Faculty of Agriculture, Cairo University, Cairo, Egypt. ⁵ UEPTEA Plateforme Tropicale d'Expérimentation sur l'Animal, INRAE, 97170 Petit-Bourg, Guadeloupe, France.

Received: 10 October 2019 Accepted: 11 February 2020

Published online: 16 March 2020

## References

1. Bishop SC (2012) Possibilities to breed for resistance to nematode parasite infections in small ruminants in tropical production systems. Animal 6:741-747

2. Silva MVB, Sonstegard TS, Hanotte O, Mugambi JM, Garcia JF, Nagda S, Gibson JP, Iraqi FA, McClintock AE, Kemp SJ, Boettcher PJ, Malek M, Van Tassell CP, Baker RL (2011) Identification of quantitative trait loci affecting resistance to gastrointestinal parasites in a double backcross population of Red Maasai and Dorper sheep. Anim Genet 43:63-71

3. Jabbar A, Iqbal Z, Kerboeuf D, Muhammad G, Khan MN, Afaq M (2006) Anthelmintic resistance: the state of play revisited. Life Sci 79:2413-2431

4. Stear MJ, Boag B, Cattadori I, Murphy L (2009) Genetic variation in resistance to mixed, predominantly *Teladorsagia circumcincta* nematode infections of sheep: from heritabilities to gene identification. Parasite Immunol 31:274-282

5. Coutinho RMA, Benvenuti CL, de Andrade Junior ALF, Silva FC, Neves MRM, do Navarro AMC, da Vieira LS, Zaros LG (2015) Phenotypic markers to characterize F2 crossbred goats infected by gastrointestinal nematodes. Small Rumin Res 123:173-178

6. Stear MJ, Henderson NG, Kerr A, McKellar QA, Mitchell S, Seeley C, Bishop SC (2002) Eosinophilia as a marker of resistance to *Teladorsagia circumcincta* in Scottish Blackface lambs. Parasitology 124:553-560

7. Dawkins HJS, Windon RG, Eagleson GK (1989) Eosinophil responses in sheep selected for high and low responsiveness to *Trichostrongylus colubriformis*. Int J Parasitol 19:199-205

8. Davies G, Stear MJ, Bishop SC (2004) Genetic relationships between indicator traits and parasitic nematode infection in sheep. Meet EAAP, Bled, Slov 5-9 Sept

9. Bambou JC, Larcher T, Cei W, Dumoulin PJ, Mandonnet N (2013) Effect of experimental infection with *Haemonchus contortus* on parasitological and local cellular responses in resistant and susceptible young creole goats. Biomed Res Int 2013:902759

10. Diez-Tascón C, Keane OM, Wilson T, Zadissa A, Hyndman DL, Baird DB, McEwan JC, Crawford AM (2005) Microarray analysis of selection lines from outbred populations to identify genes involved with nematode parasite resistance in sheep. Physiol Genomics 21:59-69

11. Keane OM, Zadissa A, Wilson T, Hyndman DL, Greer GJ, Baird DB, McCulloch AF, Crawford AM, McEwan JC (2006) Gene expression profiling of naïve sheep genetically resistant and susceptible to gastrointestinal nematodes. BMC Genomics 7:42

12. Keane OM, Dodds KG, Crawford AM, McEwan JC (2007) Transcriptional profiling of *Ovis aries* identifies *Ovar-DQA1* allele frequency differences between nematode-resistant and susceptible selection lines. Physiol Genomics 30:253-261

13. Rowe A, Gondro C, Emery D, Sangster N (2008) Genomic analyses of *Haemonchus contortus* infection in sheep: abomasal fistulation and two

Haemonchus strains do not substantially confound host gene expression in microarrays. Vet Parasitol 154:71-81

14. Knight PA, Griffith SE, Pemberton AD, Pate JM, Guarneri L, Anderson K, Talbot RT, Smith S, Waddington D, Fell M, Archibald AL, Burgess ST, Smith DW, Miller HR, Morrison IW (2011) Novel gene expression responses in the ovine abomasal mucosa to infection with the gastric nematode *Teladorsagia circumcincta*. Vet Res 42:78

15. MacKinnon KM, Burton JL, Zajac AM, Notter DR (2009) Microarray analysis reveals difference in gene expression profiles of hair and wool sheep infected with *Haemonchus contortus*. Vet Immunol Immunopathol 130:210-220

16. Andronicos N, Hunt P, Windon R (2010) Expression of genes in gastrointestinal and lymphatic tissues during parasite infection in sheep genetically resistant or susceptible to *Trichostrongylus colubriformis* and *Haemonchus contortus*. Int J Parasitol 40:417-429

17. Gossner A, Wilkie H, Joshi A, Hopkins J (2013) Exploring the abomasal lymph node transcriptome for genes associated with resistance to the sheep nematode *Teladorsagia circumcincta*. Vet Res 44:68

18. Ahmed AM, Good B, Hanrahan JP, McGettigan P, Browne J, Keane OM, Bahar B, Mehta J, Markey B, Lohan A, Sweeney T (2015) Variation in the Ovine abomasal lymph node transcriptome between breeds known to differ in resistance to the gastrointestinal nematode. PLoS ONE 10:e0124823

19. McRae KM, Good B, Hanrahan JP, McCabe MS, Cormican P, Sweeney T, O'Connell MJ, Keane OM (2016) Transcriptional profiling of the ovine abomasal lymph node reveals a role for timing of the immune response in gastrointestinal nematode resistance. Vet Parasitol 224:96-108

20. Costa V, Aprile M, Esposito R, Ciccodicola A (2013) RNA-Seq and human complex diseases: recent accomplishments and future perspectives. Eur J Hum Genet 21:134-142

21. McRae KM, Stear MJ, Good B, Keane OM (2015) The host immune response to gastrointestinal nematode infection in sheep. Parasite Immunol 37:605-613

22. Aboshady HM, Mandonnet N, Stear MJ, Arquet R, Bederina M, Sarry J, Tosser-Klopp G, Klopp C, Johansson AM, Jonas E, Bambou JC (2019) Transcriptome variation in response to gastrointestinal nematode infection in goats. PLoS ONE 14:e0218719

23. Patro R, Duggal G, Love MI, Irizarry RA, Kingsford C (2017) Salmon: fast and bias-aware quantification of transcript expression using dual-phase inference. Nat Methods 14:417-419

24. Soneson C, Love MI, Robinson MD (2015) Differential analyses for RNA-seq: transcript-level estimates improve gene-level inferences. F1000Research 4:1521

25. Rohart F, Gautier B, Singh A, Cao K-AL (2017) mixOmics: an R package for 'omics feature selection and multiple data integration. PLoS Comput Biol 13:e1005752

26. Love MI, Huber W, Anders S (2014) Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol 15:550

27. Yu G, Wang L-G, Han Y, He Q-Y (2012) clusterProfiler: an R package for comparing biological themes among gene clusters. Omi A J Integr Biol 16:284-287

28. Winer J, Jung CKS, Shackel I, Williams PM (1999) Development and validation of real-time quantitative reverse transcriptase-polymerase chain reaction for monitoring gene expression in cardiac myocytes in vitro. Anal Biochem 270:41-49

29. Pernthaner A, Shaw RJ, McNeill MM, Morrison L, Hein WR (2005) Total and nematode-specific IgE responses in intestinal lymph of genetically resistant and susceptible sheep during infection with *Trichostrongylus colubriformis*. Vet Immunol Immunopathol 104:69-80

30. Pernthaner A, Cole SA, Morrison L, Green R, Shaw RJ, Hein WR (2006) Cytokine and antibody subclass responses in the intestinal lymph of sheep during repeated experimental infections with the nematode parasite *Trichostrongylus colubriformis*. Vet Immunol Immunopathol 114:135-148

31. Robinson N, Pleasance J, Piedrafita D, Meeusen EN (2011) The kinetics of local cytokine and galectin expression after challenge infection with the gastrointestinal nematode, *Haemonchus contortus*. Int J Parasitol 41:487-493

32. Zaros LG, Neves MRM, Benvenuti CL, Navarro AMC, Sider LH, Coutinho LL, Vieira LS (2014) Response of resistant and susceptible Brazilian Somalis
crossbreed sheep naturally infected by *Haemonchus contortus*. Parasitol
Res 113:1155–1161

33. Schallig HD (2000) Immunological responses of sheep to *Haemonchus contortus*. Parasitology 120(Suppl):S63–S72

34. Urban JF Jr, Madden KB, Svetica A, Cheever A, Trotta PP, Gause WC, Katona IM, Finkelman FD (1992) The importance of Th2 cytokines in protective immunity to nematodes. Immunol Rev 127:205–220

35. Else KJ, Finkelman FD (1998) Intestinal nematode parasites, cytokines and effector mechanisms. Int J Parasitol 28:1145-1158

36. Finkelman FD, Shea-Donohue T, Goldhill J, Sullivan CA, Morris SC, Mad-
den KB, Gause WC, Urban JF Jr (1997) Cytokine regulation of host defense
against parasitic gastrointestinal nematodes: lessons from studies with
rodent models. Annu Rev Immunol 15:505-533

37. Gill HS, Altmann K, Cross ML, Husband AJ (2000) Induction of T helper 1- and T helper 2-type immune responses during Haemonchus contortus infection in sheep. Immunology 99:458–463

38. Lacroux C, Nguyen THC, Andreoletti O, Prevot F, Grisez C, Bergeaud J-P, Gruner L, Brunel J-C, Francois D, Dorchies P, Jacquiet P (2006) *Haemonchus contortus* (Nematoda: trichostrongylidae) infection in lambs elicits an unequivocal Th2 immune response. Vet Res 37:607–622

39. Arsenopoulos K, Symeonidou I, Papadopoulos E (2017) Immune and other factors modulating host resistance against gastrointestinal nematode parasites in sheep Arsenopoulos. J Hell Vet Med Soc 68:131–144

40. Almeria S, Canals A, Zarlenga D, Gasbarre L (1997) Quantification of cytokine gene expression in lamina propria lymphocytes of cattle following infection with *Ostertagia ostertagi*. J Parasitol 83:1051–1055

41. Claerebout E, Vercauteren I, Geldhof P, Olbrechts A, Zarlenga DS, Godde-
eris BM, Vercruysse J (2005) Cytokine responses in immunized and non-
immunized calves after Ostertagia ostertagi infection. Parasite Immunol
27:325-331

42. Mihi B, Van Meulder F, Vancoppernolle S, Rinaldi M, Chiers K, Van den Broeck W, Goddeeris BM, Vercruysse J, Claerebout E, Geldhof P (2014) Analysis of the mucosal immune responses induced by single and trickle infections with the bovine abomasal nematode *Ostertagia ostertagi*. *Parasite Immunol* 36:150–156

43. Caucheteux SM, Hu-Li J, Guo L, Bhattacharyya N, Crank M, Collins MT, Paul WE (2016) IL-1β enhances inflammatory TH2 differentiation. J Allergy Clin Immunol 138:898–901.e4

44. Nakao A, Afrakhte M, Morén A, Nakayama T, Christian JL, Heuchel R, Itoh S, Kawabata M, Heldin N-E, Heldin C-H, ten Dijke P (1997) Identification of Smad7, a TGFbeta-inducible antagonist of TGF-beta signalling. Nature 389:631–635

45. Bhuian AA, Li J, Wu Z, Ni P, Adetula AA, Wang H, Zhang C, Tang X, Bhuyan AA, Zhao S, Du X (2017) Exploring the genetic resistance to gastrointestinal nematodes infection in goat using RNA-sequencing. Int J Mol Sci 18:E751

46. Wang Y, Wu L, Liu X, Wang S, Ehsan M, Yan R, Song X, Xu L, Li X (2017) Characterization of a secreted cystatin of the parasitic nematode *Haemonchus contortus* and its immune-modulatory effect on goat monocytes. *Parasit Vectors* 10:425

47. Korn T, Bettelli E, Oukka M, Kuchroo VK (2009) IL-17 and Th17 Cells. Annu Rev Immunol 27:485–517

48. Benavides MV, Sonstegard TS, Van Tassell C (2016) Genomic regions associated with sheep resistance to gastrointestinal nematodes. Trends Parasitol 32:470–480

49. Schwaiger FW, Gostomski D, Stear MJ, Duncan JL, McKellar QA, Epplen JT, Buitkamp J (1995) An ovine major histocompatibility complex DRB1 allele is associated with low faecal egg counts following natural, predominantly *Ostertagia circumcincta* infection. Int J Parasitol 25:815–822

50. Sayers G, Good B, Hanrahan JP, Ryan M, Angles JM, Sweeney T (2005) Major histocompatibility complex DRB1 gene: its role in nematode resistance in Suffolk and Textel sheep breeds. Parasitology 131:403–409

**Publisher's Note**

Springer Nature remains neutral with regard to jurisdictional claims in pub-
lished maps and institutional affiliations.

Ready to submit your research? Choose BMC and benefit from:

* fast, convenient online submission

* thorough peer review by experienced researchers in your field

* rapid publication on acceptance

* support for research data, including large and complex data types

* gold Open Access which fosters wider collaboration and increased citations

* maximum visibility for your research: over 100M website views per year

At BMC, research is always in progress.

Learn more [biomedcentral.com/submissions](https://biomedcentral.com/submissions)

BMC