Multi-omics elucidated parasite-host-microbiota interactions and resistance to Haemonchus contortus in sheep

doi:10.21203/rs.3.rs-3457910/v1

Download PDF

Research Article

Multi-omics elucidated parasite-host-microbiota interactions and resistance to Haemonchus contortus in sheep

https://doi.org/10.21203/rs.3.rs-3457910/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Background:

The integration of molecular data from hosts, parasites, and microbiota can enhance our understanding of the complex biological interactions underlying parasite resistance. Haemonchus contortus, the predominant sheep parasite species in the tropics, results in significant production and economic losses, which are further compounded by the diminishing efficiency of chemical control measures due to anthelmintic resistance. Knowledge of how the host responds to infection and how the parasite, in combination with microbiota effects, modulates host immunity can guide selection decisions to breed more resistant animals. This understanding can also refine management practices and inform the development of new therapeutics for long-term helminth control.

Results:

Egg per gram (EPG) counts were obtained in Morada Nova sheep subjected to two artificial infections with H. contortus, and used as a proxy to select animals with high resistance or susceptibility. The GAST, GNLY, IL13, MGRN1, FGF14, and RORC genes and transcripts were differentially expressed between groups based on RNA-seq of the abomasum. From 50K SNP genotyping, EPG heritability estimate was 0.12, and a genome-wide association study (GWAS) identified regions on chromosomes 2 and 11 harboring candidate genes for resistance, immune response, body weight, and adaptation. Trans-eQTLs between significant variants and differentially expressed transcripts were found. Amplicon sequence variants (ASVs) from PCR amplification and sequencing of bacterial and archaeal 16S rRNA genes in sheep feces and rumen generated functional co-expression modules correlated with resistance to H. contortus, showing enrichment in pathways of response to bacterium, immune and inflammatory responses, and hub features of the Christensenellaceae, Bacteroides, and Methanobrevibacter genera, Prevotellaceae family, and Verrucomicrobiota. In RNA-seq of H. contortus, some mitochondrial, collagen- and cuticle-related genes were expressed only in parasites retrieved from susceptible sheep.

Conclusions:

This study identified chromosome regions, genes, transcripts, and pathways involved in the elaborate interactions between the sheep host, its gastrointestinal microbiota and the H. contortus parasite. These findings can assist with the development of animal selection strategies for parasite resistance and interdisciplinary approaches to control H. contortus in sheep.

GWAS

RNA-seq

16S rDNA-sequencing

eQTL

parasite resistance

Helminths are responsible for large economic and production losses in small ruminant livestock systems [1]. Haemonchus contortus, the barber pole worm, is the most prevalent and pathogenic nematode species affecting small ruminants in tropical regions [2]. Adversely, chemical treatments, which are the main approach currently used for its control, are becoming less effective due to the occurrence of anthelmintic resistance [3].

There is a need to develop new alternative control strategies to reduce the parasite burden on production systems. One interesting animal-centric control strategy is to use sheep breeds that are already naturally more resistant to gastrointestinal nematodes or, alternatively, to select for more resistant individuals within breeds (reviewed by [4]). Some indigenous breeds have already adapted to survive with low maintenance requirements and in the presence of high parasitic challenges [5]. The Morada Nova is one of these breeds, a native Brazilian hair sheep breed that has high gastrointestinal nematode resistance/resilience compared to both other tropical-adapted and high-productive sheep breeds [6]. A better molecular characterization of Morada Nova can support the use of this indigenous breed in production systems as a long-term helminth control strategy [7] and increase our understanding of the mechanisms of immune responses and parasite resistance to be used for other sheep breeds.

Knowledge of haemonchosis pathology is essential to establish new therapeutics and control strategies. This can be assisted by omics technologies, which can be used to increase our understanding of helminths and host responses [8]. Multi-omics approaches are able to aggregate different layers of information to elucidate complex host-parasite interactions [9], particularly those that are regulated by several genes with minor effects and influenced by diverse environmental conditions. Genomics can reveal genetic mechanisms underlying parasite resistance in hosts, while transcriptomics can identify the mechanisms through which hosts respond to infections or those used by parasites to establish infection and modulate the host’s immune system. Integrating both host and parasite data quantified through genomics and transcriptomics layers, such as those used in expression quantitative trait loci (eQTL) studies [10], can provide more in-depth knowledge about interactions.

Another protagonist in the parasite-host molecular scenario is the host microbiota, which has grown in importance in recent years due to the microbiome effect on resistance to gastrointestinal nematodes [11], H. contortus burden [12], and Teladorsagia circumcincta female fertility [13] in sheep. The symbiosis between microorganisms and their hosts, ruminants in particular, is crucial to provide nutrients, carbohydrates, and proteins [14]. Gastrointestinal nematode infections, however, alter the typical host microbiota, resulting in microenvironmental dysbiosis [15, 16]. It has been reported that H. contortus infections change the bacterial genera in the abomasum and rumen [15], and that T. circumcincta increases the bacterial population involved in pro-inflammatory processes, which may contribute to the immunopathology of helminth disease [16]. Furthermore, the microbiota may influence parasite and host gene expression, as they can interact with host/parasite cells and affect gene regulation, turning on or off certain cellular processes.

The aim of this work was to increase our understanding of the parasite-host-microbiota interactions that lead to increased resistance to H. contortus in Morada Nova sheep using an integrated multi-omics approach with RNA-seq data from sheep abomasum and from H. contortus, 50K sheep genotypes and 16S rRNA gene-sequencing of the microbiota from sheep feces and rumen content.

Sample collection and phenotypes

The Morada Nova animals and phenotypes used in this study were previously described [17]. Briefly, 274 three-month-old lambs (136 males and 138 females) from the Embrapa Pecuária Sudeste flock, progeny of four sires, born in two breeding seasons (March to May, 2017 and 2018), and raised under the same environmental conditions, were treated with monepantel (2.5 mg/kg; Zolvix®) to remove natural infection with gastrointestinal nematodes. After two weeks, animals were experimentally infected with 4,000 H. contortus third-stage larvae (L₃). For counts of eggs per gram (EPG), 2 g of feces collected from the rectum were mixed with 28 mL sodium chloride-saturated solution and evaluated in a McMaster chamber [18], and EPG values were obtained on days 0, 21, 28, 35, and 42 after infection by multiplying the total number of counted Trichostrongylidae eggs by 50. On day 42, lambs were dewormed with monepantel and, after 15 days, subjected to a new parasitic infection with H. contortus, following the same sampling protocol. On day 42 of the second parasitic challenge, the EPG mean was calculated from the 10 collection dates, and animals were ranked for resistance/susceptibility to H. contortus. Dams and sires were similarly evaluated in two parasite challenges (September, 2018 to January, 2019).

The additive genetic relationship matrix (A matrix) was built with the AGHmatrix [19] and pedigreemm [20] R packages (version 4.2.1 [21]) for 366 individuals (sires, dams, and lambs) of the Morada Nova flock with EPG records, and heritability (h²) was estimated with the NAM package [22].

A total of 44 lambs from the top- and bottom-ranked animals (Table 1) were split into two groups with high and low EPG means, with a 21.5x difference in EPG means between the groups. Lambs were also selected to balance the sex distribution between groups and, as much as possible, the number of offspring per sire. Venous blood samples from the four sires and the 44 lambs were collected through jugular puncture, and DNA was extracted from the leukocytes. All 48 DNA samples were then genotyped on the Illumina 50K SNP array.

Table 1

Mean and standard deviation of egg per gram (EGP) counts on collection days during two parasite challenges (day.challenge) in male (M) and female (F) Morada Nova lambs ranked as resistant (R) or susceptible (S) to *Haemonchus contortus*.
ID	Sex	Sire	Dam	0.1	21.1	28.1	35.1	42.1	0.2	21.2	28.2	35.2	42.2	Mean	SD	Rank
819	F	3	210	0	0	0	0	0	0	0	0	0	100	10	32	R
783	F	1	166	0	0	0	50	0	0	0	0	0	100	15	34	R
886	F	3	401	0	0	0	50	50	0	50	0	0	0	15	24	R
803	F	3	370	0	0	0	0	0	0	0	0	50	100	15	34	R
795	F	4	230	0	50	0	0	100	0	0	50	.	.	25	38	R
836*	F	1	289	0	50	0	0	150	100	0	100	150	200	75	75	R
897*	F	2	489	0	100	100	50	0	0	200	350	50	0	85	113	R
893	F	2	389	0	100	0	1000	350	0	0	50	0	0	150	318	R
883*	F	4	532	50	0	200	100	100	50	100	250	450	600	190	194	R
949	F	4	439	0	0	0	0	0	0	50	0	100	2100	225	660	R
818*	M	3	210	0	0	0	100	550	0	50	850	700	0	225	337	R
805	F	1	283	0	100	50	900	400	0	0	50	700	500	270	333	R
732*	F	1	382	0	0	250	1850	1450	50	0	0	450	50	410	676	R
762	M	4	7	150	0	800	150	2850	0	50	250	550	650	545	858	R
891	M	1	249	0	100	350	2950	2600	500	400	100	300	650	795	1065	R
988	M	3	333	200	400	1050	1950	2450	350	50	150	1500	500	860	841	R
789	M	4	177	0	100	150	600	1300	50	350	2100	2700	2200	955	1034	R
917	M	4	254	200	100	900	3100	1800	550	1000	1950	800	500	1090	932	R
797	M	3	403	0	300	350	700	3000	2000	1700	850	600	2100	1160	979	R
918	M	4	508	0	1750	4700	2750	100	350	1550	800	900	1400	1430	1422	R
894	M	2	389	150	400	3250	5600	4200	400	500	1600	700	850	1765	1908	R
906	F	2	503	100	700	5350	9400	11300	600	3200	2250	1500	1050	3545	3923	S
944	F	2	542	150	2950	7850	17300	8850	400	700	250	200	500	3915	5732	S
985	F	3	368	0	550	7950	21750	26450	1600	0	0	100	1200	5960	9916	S
921	M	2	122	50	8000	12500	12500	7400	1400	6350	1200	5800	6100	6130	4347	S
898	M	2	431	50	14800	17100	12000	16650	650	1200	1000	800	400	6465	7589	S
740	F	3	349	50	1450	9100	6200	4000	0	3350	15700	13200	11950	6500	5692	S
1038	F	3	252	50	7500	12100	18300	19350	400	500	700	3550	3800	6625	7461	S
907	M	3	415	400	10450	13000	12100	3850	650	4100	2550	900	18400	6640	6340	S
756	F	3	527	0	15000	14550	100	100	0	2650	14000	15600	13000	7500	7376	S
852*	M	1	391	100	3500	19000	20000	22000	0	3000	1100	8700	2750	8015	8871	S
739	M	3	349	50	3700	26250	10550	6200	150	3350	9000	9700	14200	8315	7811	S
1053*	F	1	299	0	5400	6550	18750	33550	350	500	4000	4950	9500	8355	10447	S
718	M	1	252	0	3150	12700	14950	10600	50	4700	12450	6500	20700	8580	6821	S
923*	F	1	223	0	8800	11700	15850	10200	600	5150	13600	8000	12300	8620	5309	S
742	M	1	401	0	50	12100	10550	25100	0	1450	13050	8350	20000	9065	8870	S
727	M	1	184	50	100	7350	7350	9200	0	4100	10500	25000	37200	10085	12072	S
1013*	M	3	509	300	6000	21200	20500	35000	9150	1500	3300	9800	7400	11415	10931	S
858	F	4	387	50	10450	25000	15850	38450	0	5450	8700	15000	9000	12795	11729	S
863	F	4	282	.	7300	16000	36500	56100	0	3000	3550	12200	12500	16350	18404	S
781*	M	4	339	100	150	20350	21900	8900	0	11000	36500	45950	31100	17595	16361	S
1049	F	4	486	0	1400	12400	34450	44200	5300	14450	40000	14450	.	18517	16793	S
1052	M	4	296	300	16000	28100	46150	49300	9450	12000	16100	30900	25800	23410	15772	S
1048	M	4	382	50	16250	22600	39500	42500	7350	22500	27500	81950	.	28911	24112	S

*Animals selected for RNA-seq.

For transcriptome (RNA-seq) and metabarcoding (16S rRNA) sequencing, out of the 44 lambs 10 animals – 5 resistant and 5 susceptible (showing a 54.8x difference in EPG means between the groups) were selected (Table 1). Aiming to better characterize the host response in the early phase of infection, dewormed lambs allocated to cemented stalls were subjected to a third challenge with 4,000 H. contortus L₃ and euthanized 7 days later (in the prepatent period). Abomasum fundic region samples were collected, immediately frozen in liquid nitrogen, and stored at -80°C for RNA extraction. Furthermore, 10 g of stool from the rectal ampulla were collected two weeks before slaughter, and 50 mL of rumen content were collected after slaughter of these 10 sheep and frozen in liquid nitrogen for microbiota investigation.

RNA-seq

Total RNA was extracted from the fundic region of the abomasum with QIAzol Lysis Reagent (Qiagen, Germany) using TissueRuptor (Qiagen, Germany) and purified in a silica column with RNeasy Mini kit (Qiagen, Germany). RNA concentration (ng/µL) and purity (A260/A280 ratio) were estimated in Qubit with the HS RNA kit (Thermo Fisher Scientific, USA) and in NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, USA), respectively. RNA samples were submitted to RNA-seq after quality confirmation by RNA Integrity Number ≥ 7.5 in Agilent 2100 Bioanalyzer (Agilent Technologies, USA). RNA-seq libraries were prepared with TruSeq Stranded mRNA Illumina kit (Illumina, USA) and sequenced on Illumina HiSeq2500. For each sample, approximately 40 million reads in a 100 bp paired-end (PE) protocol were sequenced.

After sequencing, fastq files were visualized with FastQC (version 0.11.7 [23]) to assess RNA-seq quality. Trimmomatic (version 0.39 [24]) was used to remove adapters, the first 10 initial bases and 2–3 final bases, reads smaller than 50 bases, and reads with phred score < 15 in a 4 bp window.

Initially, three sheep reference genomes were compared: Oar_rambouillet_v1.0 (GCA_002742125.1), Oar_v3.1 Texel (GCA_000298735.1), and ARS-UI_Ramb_v2.0 (GCF_016772045.1) (data not shown). As ARS-UI_Ramb_v2.0 (version 104) resulted in higher alignment rates and number of transcripts compared to the other ones, this genome assembly was selected as the reference for subsequent analysis. HISAT2 (version 2.1.0 [25]) was used to index the sheep reference genome and align the trimmed paired reads. The resultant SAM files were converted into BAM files with SAMtools (version 1.11 [26]), and statistics were generated with MultiQC (version 1.7 [27]). Aligned and sorted BAM files were assembled with StringTie (version 2.1.3 [28]). Gene and transcript read counts were extracted using the Python (version 3.6.6 [29]) prepDE.py3 script from the StringTie pipeline.

For differential gene expression analysis between the susceptible and resistant animals, an initial filtering step in R was used to remove lowly expressed genes with a count sum below 50 and the overexpressed genes with a mean count above 100,000. Then, the gene count matrix was analyzed in DESeq2 [30], with a Wald test adjusted p-value (q-value) < 0.1, and in edgeR [31], with a false discovery rate (FDR) < 0.05. Differential transcript expression was analyzed in ballgown [32], with q-value < 0.1, and in DESeq2, with q-value < 0.05. In all analyses both |fold change (FC) | > 1.5 and adjusted p-value were considered as thresholds. Regarding p-values, in the absence of results at < 0.05, an adjusted p-value < 0.1 was used. Differential expression for the overexpressed genes was investigated separately using t test, DESeq2, and edgeR. The biological roles of the differentially expressed and overexpressed genes were investigated with GeneCards (version 5.15 [33]) and through a literature review.

Variant calling and annotation of RNA-seq data

For variant calling from the abomasum RNA-seq data [34], trimmed PE reads were mapped to the indexed sheep reference genome using the 2-pass mode of STAR (version 2.7.3a [35]). The STAR aligner was used because of its higher tolerance for accepting mismatches and soft-clipping in comparison to HISAT2 [36]. SAMtools was used to index the BAM files, and Picard (version 2.25.0 [37]) was used to add read groups, mark duplicates, and create the sequence dictionary. GATK (version 4.2.0.0 [38]) with SplitNCigar, BaseRecalibrator including known sheep variants [39], ApplyBQSR, and HaplotypeCaller was used for calling the variants. GATK GenomicsDBImport was used to combine the files from different animals, and GenotypeGVCFs performed the joint genotyping. GATK SelectVariants was used to select the SNPs, which were then filtered with VariantFiltration with --cluster-window-size 35, --cluster-size 3, QD < 2.0, FS > 30.0, SOR > 3.0, MQ < 40.0, MQRankSum < -12.5, and ReadPosRankSum < -8.0. BCFtools (version 1.9.64 [26]) was used to retrieve the hard filtered variants and remove poor quality calls per variant (GQ < 20, depth < 3 reads, no calls, MAF < 0.05, missing genotypes > 10%, and multiallelic sites). A total of 51,434 SNPs were retained.

Afterwards, a database for the ARS-UI_Ramb_v2.0 sheep reference genome was created, and the variants were annotated with SnpEff (version 5.1 [40]). SnpSift [41] was used to detect high impact variants, and the caseControl function was applied to select variants that were homozygous between resistant and susceptible sheep (named here as RNA-seq genotypes).

GWAS

DNA from the 4 sires and 44 extreme phenotype lambs, including the 10 lambs used in the RNA-seq study, was extracted by saline precipitation and stored at -20°C. DNA integrity in 1% agarose gel electrophoresis, concentration (ng/µL) and purity (260/280 absorbance ratio between 1.8 and 2.0) in NanoDrop were assessed. DNA samples were genotyped with the Illumina OvineSNP50v3 chip. Genotyping data were subjected to quality control with snpQC [42], and SNPs with GCScore < 0.5, MAF < 0.05, HWE < 10e-16, call rate < 0.90, and samples with call rate < 0.75 were assigned as missing values. No sample with more than 10% missing values was detected. A total of 45,070 SNPs and all 48 samples were included in the following steps. The genotype frequency of SNPs with missing data was imputed by the mean frequency from the other samples. The paternity was verified from the matrix of opposing homozygous built in R [43].

These data were previously used for a case-control GWAS – genome-wide association study [44], but it was based on the Oar_v3.1 Texel reference genome. In the present study, the data from the 44 lambs were reanalyzed with the same reference genome used in the RNA-seq analysis. A custom script was used to lift the 50K genotypes to the ARS-UI_Ramb_v2.0 assembly by mapping the probe sequence to the reference and then re-aligning the nucleotides to the reference strand. EPG mean normalized by cube root transformation was used as phenotype in a snpBLUP GWAS to identify significant SNPs at FDR < 0.05, and h² for EPG was estimated from the genomic data with the NAM package in R. Only the genotypes from the arrays were used for the GWAS, as no significant results (data not shown) were obtained using the RNA-seq genotypes from the 10 lambs.

Genes located 1 Mbp upstream and downstream of each significant SNP from the GWAS analysis were retrieved from the NCBI Genome Data Viewer [45], and gene function was investigated with GeneCards and literature review to identify candidate genes of resistance to H. contortus.

eQTL analysis

The R MatrixEQTL package [46] was used for eQTL analysis using all SNP data (45,070 from the array and 51,434 from the RNA-seq) and the 10 most extreme animals with RNA-seq and genotype data. Gene and transcript expression were the RNA-seq read counts normalized by the DESeq2 variance stabilizing transformation. Sex and sire were included as covariates. Cis-eQTLs were defined as SNPs located up to 1 Mb from the regulated gene, while trans-eQTLs considered distances > 1 Mb. From all significant cis- and trans-eQTLs (FDR < 0.05) retrieved, only those comprising significant SNPs from GWAS or RNA-seq genotypes and differentially expressed data were reported here.

Analysis of unmapped RNA-seq reads

The RNA-seq study was originally designed to assess gene expression in sheep hosts, but the presence of genetic material from the H. contortus parasite in the RNA-seq data was also investigated using the reads unmapped to the sheep reference genome. Trimmed PE reads were aligned to the H. contortus reference genome (haemonchus_contortus.PRJEB506.WBPS17.genomic.fa) and assembled with the annotation file (haemonchus_contortus.PRJEB506.WBPS17.annotations.gff3) retrieved from the WormBase ParaSite database (version 17 [47]).

Due to the small amount of data and resultant low read counts, no additional filtering was performed, and differential expression analysis was implemented in DESeq2, edgeR, ballgown, and t test using the same parameters previously described for sheep.

Microbiota from sheep feces and rumen content

DNA from the stool and rumen content of the 10 lambs selected for RNA-seq was extracted with the Quick-DNA™ Fecal/Soil Microbe Miniprep Kit (Zymo, USA) following the standard protocol. PCR amplification of bacterial and archaeal 16S rRNA genes was performed with a primer set previously described [48]. Libraries were sequenced on an Illumina MiSeq (2 × 250 bp) using the Illumina V3 sequencing kit (Illumina, USA). Raw reads were quality checked, filtered for quality (> Q25), and trimmed at positions 220 (forward) and 175 (reverse) based on aggregated quality plots generated by QIIME 2 [49]. The remaining data were submitted to DADA2 [50] to resolve the sequences into amplicon sequence variants (ASVs) and exclude chimeric sequences. Bacterial and archaeal sequences were classified using the SILVA database (version 132 [51]). The conditional quantile regression (ConQuR) approach [52] was used to remove batch effects caused by multiple sequencing runs and sex. Only ASVs present in at least five lambs were considered. A total of 526 and 607 bacterial ASVs and 22 and 28 archaeal ASVs from feces and rumen, respectively, were considered for downstream analysis.

Multi-omics integration

The final step of the analysis was the integration of fecal and ruminal microbiota abundance identified in the 16S rRNA data with the expression identified in RNA-seq. A co-expression network analysis between ASVs, host genes/transcripts, and H. contortus gene/transcript expression was performed using the R CEMiTool package [53]. CEMiTool uses an unsupervised filtering method based on the inverse gamma distribution for each gene selection used in the analyses and employs a soft-thresholding power β to determine a similarity criterion between pairs of features [53]. The CEMiTool parameters were apply_vst = TRUE, cor_method = "spearman", network_type = "signed", and tom_type = "signed" to construct co-expression networks. The minimum module size was set to 50 or 100, for gene or transcript expression levels, respectively. Genes within significantly associated modules were then subjected to Metascape [54] for biological pathway and gene ontology enrichment analysis, with the input species set to Homo sapiens.

From the A matrix of the pedigree data, the inbreeding coefficient was estimated at 0.0088, and EPG h² = 0.05.

RNA-seq

The sheep abomasum RNA-seq reads ranged from 31.1 M to 41.3 M per sample, with 97.9–98.5% of the reads mapped to the sheep ARS-UI_Ramb_v2.0 reference genome (64.5 to 70% of proper pairs). For the differential gene expression analysis, 9,099 genes with zero counts and 5,891 genes with count sums below 50 were removed from the total of 33,800 genes. Fourteen overexpressed genes (Table 2) were removed and analyzed separately, but no differentially expressed gene (DEG) was detected between groups. Then, 18,796 genes (55.6%) remained for subsequent analysis.

Table 2

Overexpressed genes in RNA-seq from abomasum of Morada Nova sheep.
Gene ID	Gene name
LOC101105864	Pepsin A
LOC443320	Lysozyme C-1-like
LOC443322	Lysozyme 3a precursor
COX1	Cytochrome c oxidase subunit I
PIGR	Polymeric immunoglobulin receptor
COX3	Cytochrome c oxidase subunit III
LGALS15	Lectin, galactoside-binding, soluble, 15
LOC114108841	Immunoglobulin lambda variable 2–14
LOC101122151	Intelectin 2
ATP6	ATP synthase F0 subunit 6
ND4	NADH dehydrogenase subunit 4
LOC121818276	NADH-ubiquinone oxidoreductase chain 1-like
LOC101107475	Ig alpha-1 chain C region-like
COX2	Cytochrome c oxidase subunit II

Eleven DEGs (Table 3) were detected between H. contortus-susceptible and H. contortus-resistant animals using edgeR (Table S1 and Figure S1). Four genes (CCDC85B, LOC121819799, GAST, and LOC101121371) were upregulated in H. contortus-resistant sheep, whereas seven genes (LOC114116426, LOC101116991, LOC101104728, LOC101107420, GNLY, GRP, and IL13) were upregulated in susceptible sheep. Using DESeq2 (Table S2), only two DEGs were detected as upregulated in susceptible sheep (LOC101116991 and LOC114116426; both were also DEGs in the edgeR results) (Table 3).

Table 3

Differentially expressed genes in the abomasum between *H. contortus* susceptible and resistant Morada Nova sheep using edgeR and DESeq2.
Gene ID	Log FC	FDR/ q-value	Gene name
edgeR
CCDC85B	-8.95	0.0011	Coiled-coil domain containing 85B
LOC121819799	-6.67	0.0060	ncRNA
GAST	-13.43	0.0158	Gastrin
LOC101121371	-3.15	0.0257	60S ribosomal protein L37a
*LOC114116426	7.33	0.0004	Protein C1orf43 homolog (pseudogene)
*LOC101116991	2.17	0.0060	Ribonuclease K6-like
LOC101104728	8.36	0.0060	Antimicrobial peptide NK-lysin
LOC101107420	3.98	0.0074	Phospholipase A2, membrane associated
GNLY	3.45	0.0141	Granulysin
GRP	9.09	0.0232	Gastrin releasing peptide
IL13	2.10	0.0474	Interleukin 13
DESeq2
*LOC101116991	2.17	0.060	Ribonuclease K6-like
*LOC114116426	6.81	0.055	Protein C1orf43 homolog (pseudogene)

*Detected in both edgeR and DESeq2 analysis. Negative values of log fold change (FC) indicate upregulation in resistant sheep.

A subset of ten genes (expressed in at least two samples and count sum > 50 reads) were found to be exclusively expressed in either the susceptible or resistant animals. Of these, six genes (CCDC85B, LOC121819799, GAST, ABCC11, LOC101106468, and LOC121818066) were expressed only in resistant sheep, and four genes (LOC101104728, GNGT1, LOC114116426, and GRP) were expressed only in susceptible sheep (Table 4). Six of these (CCDC85B, LOC121819799, GAST, LOC101104728, LOC114116426, and GRP) were previously identified as DEGs (Table 3).

Table 4

Mean counts per million (CPM) of genes exclusively expressed in the abomasum of Morada Nova sheep resistant or susceptible to *H. contortus*.
Gene ID	CPM in resistant	CPM in susceptible	Gene name
*CCDC85B	4.06	0.00	Coiled-coil domain containing 85B
*LOC121819799	0.82	0.00	ncRNA
*GAST	87.29	0.00	Gastrin
ABCC11	0.53	0.00	ATP binding cassette subfamily C member 11
LOC101106468	0.74	0.00	Transcription initiation factor TFIID subunit 9-like
LOC121818066	1.16	0.00	ncRNA
*LOC101104728	0.00	2.67	Antimicrobial peptide NK-lysin
GNGT1	0.00	0.75	G protein subunit gamma transducin 1
*LOC114116426	0.00	1.31	Protein C1orf43 homolog
*GRP	0.00	4.41	Gastrin releasing peptide

*Detected in edgeR and/or DESeq2 analysis.

From the total of 76,900, 478 differentially expressed transcripts (DETs) from 435 genes (including the LOC101116991, LOC114116426, LOC101107420, and IL13 DEGs) were detected by DESeq2 (Table S3), whereas using ballgown (Table S4 and Figure S2), only three of them (transcripts of TECPR1, LOC114111361, and LOC121818463 genes) were detected.

Variant calling and annotation of RNA-seq data

A total of 884,918 variants were detected in the sheep abomasum RNA-seq data. After filtering, 51,434 autosomal SNPs remained, and they were associated with 264,604 effects, of which 0.036% was of high impact and 24.79% were missense effects. Case-control selection of RNA-seq genotypes between resistant and susceptible sheep detected 16 variants (Table S5) in 3 genes: TRAPPC6B (Chr18), MGRN1 (Chr24), and SPCS3 (Chr26).

The genes from differential expression and RNA-seq genotypes are related to response or resistance to gastrointestinal nematodes and parasites (e.g., GAST, GNLY, IL13, MGRN1, FGF14, and RORC), growth, body conformation, and weight (e.g., LRRC8B and IGF1), microbiota (e.g., RAPGEF2), and adaptation (e.g., ABCC11, GNGT1, DTX3, and ZNF789).

GWAS

After quality control filtering of the genotype data, 45,070 SNPs from the 44 lambs were used in the snpBLUP GWAS. The observed SNP average heterozygosity was 0.38, and the diagonal mean value of the genomic relationship matrix was 0.95, confirming no inbreeding in the population. The opposing homozygotes matrix identified two discordant sire-offspring pairs that had to be reassigned. In the final data of 44 lambs, 11 were offspring of sire 1 (with a 28.1x difference in EPG means between resistant and susceptible half-siblings), 7 of sire 2 (7.5x difference in EPG means), 13 of sire 3 (19.9x difference in EPG means), and 13 of sire 4 (30.8x difference in EPG means). EPG h² = 0.12 was estimated using the GRM.

The snpBLUP GWAS (Fig. 1) detected three significantly associated SNPs (FDR < 0.05): rs419988472 (intron variant in the DHRS9 gene) at position Chr2:140323047, OAR11_27473453.1 (rs161041632 - exon variant in the MED11 gene) and s08310.1 (rs427555933 - intergenic variant) at positions Chr11:26501006 and Chr11:26595120, respectively, which overlapped in Fig. 1.

Figure 1. snpBLUP GWAS for EPG in Morada Nova sheep. Dashed line: FDR < 0.05.

Within the ± 1 Mbp window of the significant SNP, there were 171 annotated genes (Table S6). The candidate genes included genes related to parasite resistance (e.g., ALOX12, ALOX15, CXCL16, DHRS9, DNAH2, GP1BA, KDM6B, LRP2, NLRP1, RABEP1, and TNFSF13), immune response (e.g., ALOX12B, CD68, and CHD3), microbiota (e.g., ABCB11 and NLGN2), muscle development and growth (e.g., ENO3, TRNAC-GCA, and WSCD1), and adaptation (e.g., DVL2, SLC2A4, and TRNAG-CCC).

eQTL analysis

For gene expression, 48,682 significant (FDR < 0.05) eQTLs consisting of 426 cis-eQTLs (Table S7) and 48,256 trans-eQTLs (Table S8) were detected, while transcript expression resulted in 9,999,256 eQTLs: 2,540 cis-eQTLs (Table S9) and 9,986,716 trans-eQTLs (Table S10). From the significant eQTLs, there was no DEG regulated by significant variants, but trans-eQTLs were detected for DETs (Table 5), including for transcripts of functional candidate DEGs, such as DTX3, FGF14, LRRC8B, RAPGEF2, RORC, and ZNF789.

Table 5

Significant eQTLs between SNP variants and differentially expressed transcripts (DET) upregulated in *H. contortus*-resistant or -susceptible Morada Nova sheep.
SNP (gene - chromosome)	DET (gene - chromosome)	Upregulated
OAR11_27473453.1 and s08310.1 (MED11 and intergenic - Chr11)	XM_042241996.1 (TMEM187 - ChrX)	Susceptible
NC_056071.1:46893094_C/A (TRAPPC6B - Chr18)	XM_042256817.1 (RORC - Chr1), XM_042243310.1 (LOC121818624 - Chr2), XM_027967527.2 (DTX3 - Chr3), XM_042249264.1 (AGAP3 - Chr4), XM_027970457.2 (THG1L - Chr5), XM_042252750.1 (PTGR2 - Chr7), XM_042255599.1 (CHRNE - Chr11), XM_042230308.1 (CDC123 - Chr13), XR_006056188.1 (LOC101122718 - Chr14), XM_042234867.1 (C17H22orf15 - Chr17), XM_042237011.1 (BRD2 - Chr20), XR_003586743.2 (KDSR - Chr23), XM_042240704.1 and XM_042240703.1 (TECPR1 - Chr24)	Resistant
NC_056071.1:46893094_C/A (TRAPPC6B - Chr18)	XM_015092088.3 (LRRC8B - Chr1), XM_042251691.1 (HGFAC - Chr6), XM_015097958.3 (FGF14 - Chr10), XM_027977036.2 (FAM110A - Chr13), XM_042234479.1 and XM_042234478.1 (EP400 - Chr17), XM_004021043.5 (ZNF789 - Chr24), XR_003587432.2 (LOC114111361 - ChrX), XR_006058459.1 (LOC121818463 - unknown)	Susceptible
NC_056077.1:4352395_C/A and NC_056077.1:4352403_C/T (MGRN1 - Chr24)	XM_042234385.1 (RAPGEF2 - Chr17)	Resistant
NC_056079.1:6981437_A/G (SPCS3 - Chr26)	XM_015097958.3 (FGF14 - Chr10)	Susceptible

Unmapped read analysis (RNA-seq from H. contortus)

Between 0.4% and 1.0% of the total abomasum RNA-seq data reads mapped to H. contortus, resulting in 0.1 to 0.4 M of reads/sample. Most of the 19,776 annotated genes had a zero count, and only 40 genes (0.2%) were retained for further analysis.

The most highly expressed gene, which was detected in all samples, was HCON_00026760, a non-annotated pseudogene. No DEG using DESeq2 (Table S11), edgeR (Table S12), and t test, and no DET using DESeq2 (Table S13) and ballgown (Table S14) were detected between the H. contortus reads retrieved from resistant and susceptible Morada Nova sheep.

Of the 40 genes, 26 were expressed exclusively in H. contortus recovered from susceptible Morada Nova sheep (Table S15), but only four genes were expressed in at least two samples: HCON_00174250 (Epsin-1), HCON_00667215 (Cytochrome c oxidase subunit 1), HCON_00667245 (NADH:ubiquinone reductase (H(+)-translocating)), and HCON_00667255 (NADH dehydrogenase subunit 6). Despite the expression in only one sample each, three genes related to collagen and cuticle composition were expressed exclusively in H. contortus recovered from susceptible sheep (Table S15). For H. contortus genes collected from resistant sheep, nine genes were exclusively expressed (Table S15), but only in one sample each.

Gene co-expression network analysis

The co-expression network analysis detected four functional modules based on the sheep RNA-seq gene data and ASVs (Fig. 2). Among them, three modules (Fig. 2 and Table 6) exhibited a significant association with resistance to H. contortus, while one module was not associated with the trait. Notably, the most active module, referred to as M1, encompassed a total of 302 correlated features (being 77 ASVs) and displayed a positive normalized enrichment score (NES) of 4.87 in the resistant group. Conversely, M2, consisting of 284 features (being 117 ASVs), exhibited a negative NES of -3.75 in the resistant group, and M3, containing 240 features (being 36 ASVs), also exhibited a negative NES of -2.32 in the resistant group.

Figure 2. Co-expression network analysis of sheep RNA-seq genes and amplicon sequence variants (ASVs). The figure displays the normalized enrichment score (NES) of modules (red represents higher activity and blue represents lower activity), with circle size and color intensity proportional to the NES values.

Table 6

Significant functional co-expression modules of sheep RNA-seq genes and amplicon sequence variants (ASVs).
Module	Number of features	Hubs
M1	Genes = 225; archaea stool = 3; archaea rumen = 5; bacteria stool = 27; bacteria rumen = 42	LOC114116187 (CALHM6), CD8B, LOC114109611 (UBD), CD8A, LOC101122689 (TRGC1)
M2	Genes = 167; archaea rumen = 9; bacteria stool = 54; bacteria rumen = 54	ASV_Bac_1273_R, ASV_Bac_422_S, CCDC150, ASV_Bac_396_R, HKDC1
M3	Genes = 204; archaea stool = 1; archaea rumen = 4; bacteria stool = 16; bacteria rumen = 15	SLC26A7, HRH2, TNNI3, KCNJ16, PPP1R1A

Hub genes/genera representative of each of the three modules associated with H. contortus resistance were identified (Table 6 and Tables S16 and S17). In the case of M1, the hub features were genes associated with the immune system, such as CALHM6, CD8B, CD8A, and TRGC1. On the other hand, M2 had three ASVs as hub features, i.e., ASV_Bac_1273_R, classified as belonging to the Christensenellaceae R-7 group genus, ASV_Bac_396_R, classified as the Kiritimatiellae class (both found in the rumen), and ASV_Bac_422_S, classified as the Bacteroides genus, from stool samples.

The annotated genes were subjected to pathway analysis to capture the biological information through enriched terms (Fig. 3 and Table S18). For genes from M1, the top pathways and processes included response to bacterium and positive regulation of immune and inflammatory response. For genes from M2, the enriched terms were related to regulating insulin-like growth factor (IGF) transport and uptake by insulin-like growth factor binding proteins (IGFBPs), post-translational protein phosphorylation, and degradation of the extracellular matrix. Additionally, for genes in M3, gastric acid secretion was highlighted.

Figure 3. Metascape enrichment analysis of statistically enriched ontology terms of hub genes/genera from three functional co-expression modules of sheep RNA-seq genes and amplicon sequence variants (ASVs). Heatmap of enriched terms across input gene lists, colored by p-values.

An enrichment of genera classified as Prevotella, Treponema, Christensenellaceae R-7 group, and Methanobrevibacter was found in the M1 and M2 modules. Four ASVs were found in both the rumen and stool for M1, ASV_98 and ASV_23, classified as the p-251-o5 and [Eubacterium] coprostanoligenes group bacterial families, and two archaea, i.e., ASV_13 (Methanobrevibacter genus) and ASV_6 (Candidatus Methanomethylophilus). For M2, five archaeal ASVs, i.e., ASV_3, ASV_4, ASV_8, ASV_31, and ASV_57, all classified as Methanobrevibacter, were enriched.

Transcript co-expression network analysis

The co-expression network analysis with the RNA-seq transcript data identified 26 functional co-expression modules (Figure S3 and Table S19) and hub features (Table S20). The top modules and hub features (Table 7) were M2 (294 features with 55 ASVs, NES = -3.85), M25 (123 features with 7 ASVs, NES = -3.29), M6 (263 features with 31 ASVs, NES = 3.66), and M10 (231 features with 8 ASVs, NES = 3.4).

Table 7

Top significantly functional co-expression modules of sheep RNA-seq transcripts and amplicon sequence variants (ASVs).
Module	Number of features	Hubs
M2	Isoforms = 239; archaea stool = 2; archaea rumen = 3; bacteria stool = 24; bacteria rumen = 26	ASV_Bac_396_R, ASV_Bac_250_R, XM_042238037.1 (SYTL2), ASV_Bac_1273_R, ASV_Bac_422_S
M6	Isoforms = 232; archaea stool = 3; archaea rumen = 4; bacteria stool = 11; bacteria rumen = 13	ASV_Bac_1254_R, ASV_Bac_433_R, ASV_Bac_115_R, ASV_Bac_269_R, ASV_Bac_218_S
M10	Isoforms = 223; archaea rumen = 1; bacteria rumen = 7	XM_027958292.2 (SFMBT1), XM_042256129.1 (MAPT), XM_042230866.1 (ZNF570), ASV_Bac_239_R, XM_027964183.2 (SEMA4D)
M25	Isoforms = 116; bacteria stool = 4; bacteria rumen = 3	XM_042250691.1 (PWWP2A), XM_027978299.2 (COQ8B), XM_012180311.3 (PTPN13), XM_042247173.1 (RIC8B), XM_042237043.1 (SSR1)

Regarding M2, several ASVs were identified as hubs, including the same ASVs found in the gene-level M2 module. ASV_Bac_250_R, classified as Verrucomicrobiota, was also identified as a hub for M2 at the transcript level. For M6, a total of five ASVs were identified as hubs, including ASV_Bac_1254_R (Prevotella genus), ASV_Bac_433_R (Ruminococcus gauvreauii group genus), ASV_Bac_115_R (Prevotellaceae NK3B31 group genus), ASV_Bac_269_R (Prevotella genus), and ASV_Bac_218_S (Lachnospiraceae family).

Regarding pathway enriched terms (Fig. 4 and Table S21), the genes associated with each transcript showed enrichment in MAPK, Wnt, and EGF/EGFR signaling pathways, morphogenesis, and cellular processes, including regulation of cell morphogenesis, chromatin organization, and actin cytoskeleton organization. In addition, transcripts representing the same genes (i.e., ABO, ASAP2, ADGRF5, ARHGAP12, CAPN15, and FN1) were found in different modules.

Figure 4. A) Metascape enrichment analysis of statistically enriched ontology terms. Heatmap of enriched terms across input gene lists, colored by p-values. B) Overlap between gene lists - identical genes are linked with purple curves, genes that hit multiple lists are colored dark orange, and genes unique to a list are shown light orange.

None of the H. contortus genes or transcripts correlated significantly with the identified modules.

The origin of the indigenous Brazilian Morada Nova sheep breed is still under investigation. Although it was assumed to have descended from the Bordaleira Portuguese breed [55], a study using genomic data showed no close relationship with European or African breeds and suggested that it originated from a single migration event of Sidaoun Algerian hair sheep [56]. The animals used in this study were raised in the state of São Paulo, which holds 10.42% of the Morada Nova population registered in Brazil [57]. There was no evidence of inbreeding, calculated from either the pedigree or the genomic data, in these animals. High inbreeding levels have, however, been reported in animals from the northwest region of Brazil [56]. Production conditions under extensive or semi-intensive system and animal selection practices vary substantially between Brazilian states and consequently, so does the genetic background of the breed, e.g., as reported for frequencies of a SNP associated with litter size between Morada Nova from Ceara (a state from the Brazil's northwest region) and São Paulo [58]. Nonetheless, Morada Nova animals have enough genetic variation to be used in breeding programs [56].

Reported EPG heritability estimates vary considerably between 0.01 and 0.65, as it is a complicated trait that is easily affected by breed, age, parasite, climate, natural versus artificial parasite infection, and other experimental conditions (reviewed by [4]). Here, the estimated h² of EPG was 0.12 with the GRM and 0.05 from the pedigree, noting the very limited number of animals used to estimate h². Further studies with larger numbers of individuals are warranted to address the reliability of selecting Morada Nova individuals resistant to H. contortus based solely on EPG. On the other hand, even if not used for animal selection, the genomic estimated breeding value (GEBV) for EPG is considered a good criterion to identify animals to target for selective anthelmintic treatment [59], which can contribute to the control of gastrointestinal nematodes in flocks.

Regarding the sheep RNA-seq data, the overexpressed genes we identified were already previously reported in abomasum transcriptome studies, are involved in relevant gastrointestinal biological processes and responses to infection, and showed no differential expression between resistant and susceptible animals. The PIGR, COX1, and LGALS15 genes were detected in sheep abomasum in a study comparing several tissues [60], and PIGR was the most abundant transcript found in cattle abomasum [61]. Pepsin is the principal acid protease of the stomach [62]; lysozymes are involved in bacterial cell wall cleavage [63]; COX1 produces cytoprotective prostaglandins for the stomach [64]; LGALS15 (also known as ovgal11), associated with immune/inflammatory responses and protection against infection [65], and intelectin-2, associated with the gastrointestinal mucus [66], were detected in the abomasum of sheep infected with H. contortus; immunoglobulin lambda expression was detected in gastric cancer [67]; and ND4 mutations are considered biomarkers of preneoplastic lesions of the gastrointestinal tract [68].

Among the differentially expressed genes, upregulation of the IL13 gene in H. contortus-susceptible Morada Nova is in accordance with a previous study using infected young sheep of a native breed [69]. However, IL13 was also detected as upregulated in adult resistant sheep [70]. In addition to animal age, the infection length should be considered, as genes related to the inflammatory response, T lymphocyte attraction, and leukocyte binding are upregulated in the abomasal lymph nodes of resistant animals at the beginning of T. circumcincta infection but downregulated at later stages, evidencing a delayed Th2 response in susceptible animals [71]. The roles of IL13 in resistance to H. contortus are through smooth muscle hypercontraction and increased mucus production, resulting in detachment of parasites from the gut wall [70].

Other differentially expressed genes and transcripts detected here were previously associated with the response to H. contortus, such as GAST [72] and FGF14 [73] in sheep and RORC in goats [74], and to Trichostrongylus colubriformis, such as GNLY [75]. Some of the genes were associated with weight, body conformation, and growth traits, such as LRRC8B [76] and IGF1 [77]; adaptation, such as ABCC11 for desert adaptability [78], GNGT1 for hypoxia [79], ZNF789 for disease resistance and climate adaptation [80], and DTX3 for tropical adaptation in cattle [81]; and gut microbiota composition in mice, such as RAPGEF2 [82]. MGRN1, identified in the RNA-seq genotypes, was a candidate gene identified in a GWAS for gastrointestinal nematode resistance in sheep [83].

Candidate genes in these regions were associated with resistance to gastrointestinal nematodes, such as NLGN2 in goats [84], CXCL16 and CD68 in cattle [85], specific response and resistance to H. contortus, such as ALOX15 [69] and LRP2 [86], to T. circumcincta, such as RABEP1 [87], to the cestode Echinococcus granulosus, such as TNFSF13 [88], and to the protozoan Toxoplasma gondii in mice, such as NLRP1 [89]. In addition, candidate genes were related to macrophage regulation during helminth infection, such as KDM6B [90], response against O. ostertagi vaccination in cattle, such as DHRS9 [91], and lipoxygenases (ALOXE3, ALOX12, ALOX12B, and ALOX15 genes) associated with the inflammatory response in mice [92]. There were also genes related to body size and feed efficiency in sheep, such as ALOX12 and TP53 [93], and in cattle, such as TP53 [94], WSCD1 [95], and TRNAC-GCA [96], to adaptability, such as SLC2A4 for hypoxia [97] and CHD3 for disease resistance and climate adaptation in sheep [80], and DVL2 [98] for thermal adaptation in cattle. From the candidate genes involved in microbiota regulation, ABCB11 (also known as BSEP) has a role in the transport of bile salts, whose direct microbial activity shapes the gut microbiota [99], and NLGN2 gene expression in horses is affected by anthelmintic treatment, probably through an effect on the microbiota [100].

The lack of cis-eQTL for the significant variants suggests the absence of regulation among genes closely located on the same chromosome; however, the small sample size of this study has limited power to identify smaller effects – a larger sample size generally results in gains for cis-eQTL mapping [101]. Trans-eQTLs were detected between significant 50K or RNA-seq genotypes and differentially expressed transcripts, validating the relationship between genomic and transcriptomic data. Additionally, the identification of genomic regions controlling the differential expression of genes and transcripts may contribute to the comprehension of molecular responses associated with resistance to H. contortus in sheep and be used for animal selection.

The Epsin-1 gene was expressed in H. contortus recovered from susceptible sheep. Epsin-1 is an ENTH-domain protein involved in endocytosis and lysosomal protein trafficking, whose silencing in Heterodera avenae, a cereal cyst nematode, resulted in a 71% reduction in females and eggs [102]. Cytochrome c oxidase, NADH:ubiquinone reductase, and NADH dehydrogenase were also expressed. These are mitochondrial genes, which are considered potential targets for anthelmintic treatments due to the unique energy transducing and anaerobic systems developed by nematode parasites in their adaptation to the low oxygen concentration of the mammalian host gastrointestinal tract [103, 104]. Additional genes were cuticle- and collagen-related; the nematode cuticle interacts with the host immune system, and its major protein is collagen [105]. Genes associated with collagen and cuticle development were upregulated in the transition from the L3 to L4 stage in H. contortus [106]. Potential roles of epsin-1, mitochondrial, collagen and cuticle genes in the establishment and maintenance of infection and parasite fitness in sheep hosts should be further investigated and considered as potential targets for the development of new therapeutics against H. contortus.

The early establishment of the gut microbiome is crucial for immune system development and for the maintenance of a healthy gut, including barrier function and mucosal immunity (reviewed by [107]). In fact, the presence of some microorganisms, such as Prevotella, might stimulate the production of specific inflammatory cytokines [108], whereas others, such as Clostridium, have anti-inflammatory effects (reviewed by [109]). Moreover, the nematodes also modulate the microbiome to provide an adequate environment for their survival [110].

The M1 module of the gene-ASV network with higher activity (i.e., which contains many gene-AVSs that present a strong and positive correlation with each other and with the phenotype) in resistant sheep was related to positive regulation of the immune system, response to bacteria, and inflammatory response, which are traits intrinsically related to parasite resistance. Nematodes evolved immunomodulatory mechanisms to suppress host immune responses and promote infection [111]. Knowledge about the mechanisms and target molecules involved in the inflammatory response may provide effective means to control nematode parasites [112]. Four of the five hub genes played a role in immune responses. The CD8B and CD8A genes encode subunits of the CD8 protein found on the surface of cytotoxic T cells. A DNA vaccine conferring partial protection against H. contortus infection in goats resulted in increased CD8 + T lymphocyte production and a reduction in EPG and worm burden in the abomasum [113]. TRGC1 encodes the gamma-constant region of the T-cell receptor surface that recognizes and binds to specific antigens, playing a multifaceted role in tissue homeostasis, autoimmunity, pro- and anti-tumor activity, and innate and adaptive immune responses [114, 115]. However, no relationship between TRGC1 and nematode infection has been described so far.

SLC26A7 and KCNJ16 were hub genes in the M3 module, which had lower activity in resistant animals. These genes were related to homeostasis and pH balance [116] and, in addition to HRH2, to impaired gastric acid secretion [117, 118], an enriched pathway from M3. Some studies suggest that gastric acid protects against nematode infections, as shown by reduced gastric acid secretion and predisposition to infection with a variety of nematodes, including Ostertagia and T. colubriformis in sheep [119, 120]. H. contortus infection causes significant increase in abomasal pH at the early and late stages of infection in goats [121]. Our findings indicate a correlation between gastric acid secretion and susceptibility to H. contortus. It is plausible to infer that in comparison to resistant animals, susceptible animals are less able to respond to H. contortus infection through gastric acid secretion, which, in turn, results in higher infection rates. The HKDC1 hub gene from M2, which had lower activity in resistant animals, is a protein related to glucose use and homeostasis [122].

Among the ASVs, the p-251-o5 family, identified in M1, was described in the fecal microbiota of pigs [123] and in the rumen of ewes and cows [124, 125]. In addition, the [Eubacterium] coprostanoligenes group, an anaerobic bacteria potentially impacting host lipid metabolism [126], had lower abundances in Tibetan sheep infected with gastrointestinal nematodes [127], corroborating our results and suggesting that it could be a promising target for H. contortus resistance. Christensenellaceae R-7, a hub feature in M2, was associated with amino acid and lipid metabolism and with human metabolic health in different disease contexts, including obesity and inflammatory bowel disease [128]. The relative abundance of Christensenellaceae was reported to decrease in goats infected with H. contortus [121]. The Bacteroides genus, an ASV from stool in the M2 module, is considered beneficial to hosts by preventing potential pathogens from colonizing the gut and by processing complex molecules to simpler ones [129].

The hub archaeal ASVs classified as Methanobrevibacter genera in the M2 module play a role in the gut microbial ecosystem [130], contributing to the efficient digestion of polysaccharides by consuming the end products of bacterial fermentation [131] and being considered potential therapeutic targets for managing gastrointestinal disorders [132]. Anthelmintic treatment of adult ewes significantly impacted the archaeal community, resulting in increased relative abundances of different Methanobacteria [133]. An increase in Methanobrevibacter occurred during chronic Trichuris trichiura infection in humans [134]. The hypothesis is that H. contortus infection could increase mucus secretion [135], which would provide energy for adapted microorganisms, including the mucus colonizer Methanobrevibacter [134]. Thus, the potential role of Methanobacteria is of interest for further research to control H. contortus infection.

The M2 and M25 modules of the transcript networks were less active in resistant animals. The ASV classified as Verrucomicrobiota, a hub feature in M2, was associated with mucin-degradation, glucose homeostasis, and immunity regulation [136]. Previous studies described a symbiotic association between Verrucomicrobia and soil nematodes [137], and the abundance of the Verrucomicrobiota phylum was decreased in the abomasum of lambs infected with H. contortus [138], suggesting potential mechanisms for therapeutic interventions.

Enriched ASV hubs in M6, with higher activity in the resistant group, were from the Prevotellaceae family (ASV_1254, ASV_115, and ASV_269). Previous studies showed an increase in Prevotella abundance in the ruminal, fecal, and abomasal microbiota of sheep and goats infected with H. contortus [15, 139]. Here, a decrease in Prevotella abundance associated with H. contortus susceptibility was observed, which is consistent with more recent findings after H. contortus infection [121, 138] and after other nematode infections in humans [140] and mice [141]. In addition, as the Prevotella genus are known to specialize in complex polysaccharide degradation, such as starch and cellulose, and contribute to the metabolism of dietary fiber [142], its reduction may affect feed digestion, resulting in low weight gain in infected animals [141] and, as a consequence, lower resilience to parasites.

Gene and transcript expression data and genomic markers from Morada Nova sheep were investigated and associated through eQTL studies to assess their roles in the phenotype of H. contortus resistance. Some candidate genes were related to immune response and parasite resistance (ALOX12, ALOX12B, ALOX15, CD68, CXCL16, DHRS9, DNAH2, FGF14, GAST, GNLY, GP1BA, IL13, KDM6B, LRP2, MGRN1, NLRP1, RABEP1, RORC, and TNFSF13), growth and weight (ENO3, IGF1, TP53, LRRC8B, TRNAC-GCA, and WSCD1), environmental adaptation (CHD3, DTX3, DVL2, GNGT1, SLC2A4, TRNAG-CCC, and ZNF789), and microbiota regulation (ABCB11, NLGN2, and RAPGEF2). Through eQTL mapping, SNP variants were predicted to regulate the differential expression of transcripts, including some of the candidate genes, such as DTX3, FGF14, LRRC8B, RORC, and SNF789. These genes are potential molecular markers to detect and select H. contortus-resistant animals.

Transcriptomic and genomic data from hosts, expression data from H. contortus, and ASVs from the microbiota of feces and rumen were integrated to detect biological functions and pathways resulting in resistance. The genes (mitochondrial and collagen- and cuticle-related genes) expressed in H. contortus may modulate evasion from the host’s immune system, facilitating the establishment of the infection. Furthermore, sheep genes (CD8A, CD8B, HKDC1, KCNJ16, SLC26A7, and TRGC1) and biological pathways (positive regulation of immune system, response to bacterium, inflammatory response, gastric acid secretion, and mucus secretion) identified through this multi-omics approach are potential modulators of host immunity. Additionally, enriched microbiota (Bacteroides genus, Christensenellaceae R-7, [Eubacterium] coprostanoligenes group, Methanobrevibacter genera, Prevotellaceae family, and Verrucomicrobiota) may regulate host metabolism (homeostasis, glucose use, amino acid and lipid metabolism, and dietary fiber degradation), resulting in resilience to parasite infection. The parasitic genes, biological pathways, and microbiomes detected here are potential targets for the development of new therapeutics aiming to increase host resistance and parasite control in sheep flocks.

Ethics approval and consent to participate

All experimental protocols were approved by the Embrapa Pecuária Sudeste Ethics and Animal Experimentation Committee under protocol numbers 04/2017 and 01/2020. All methods are reported in accordance with ARRIVE guidelines for the reporting of animal experiments.

Consent for publication

Not applicable.

Availability of data and materials

The datasets generated and/or analyzed during the current study are available in the European Nucleotide Archive repository, under the project PRJEB67593, https://www.ebi.ac.uk/ena/browser/view/PRJEB67593.

Competing interests

The authors declare that they have no competing interests.

Funding

FAPESP (processes no. 2017/01626-1, 2019/02929-3, 2021/11842-9, 2021/02535-5, and 2022/06281-0). AMGI is recipient of a productivity fellowship from CNPq. BGNA was funded by Science Foundation Ireland and co-funded by the European Regional Development Fund through the ADAPT Centre for Digital Content Technology Resource Centre grant number 13/RC/2106_P2.

Authors' contributions

SCMN generated, analyzed and interpreted genomic and transcriptomic (GWAS, variant calling, eQTL, and RNA-seq) data. TFC analyzed and interpreted microbiota data and performed multi-omics integration. AMGI and CRO generated and interpreted transcriptomic data. BGA and LCAR generated and interpreted microbiota data. MVB, ACSC, SNE and APM obtained phenotypic data and collected samples. CG analyzed and interpreted genomic and transcriptomic data. All authors read and approved the final manuscript.

Acknowledgements

We thank the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES), Brazil—Finance Code 001 for the free access to the journals used in the literature review.

Chagas ACS, Tupy O, Santos IB, Esteves SN. Economic impact of gastrointestinal nematodes in Morada Nova sheep in Brazil. Rev Bras Parasitol Vet. 2022;31:e008722.
Amarante AFT, Bricarello PA, Rocha RA, Gennari SM. Resistance of Santa Inês, Suffolk and Ile de France sheep to naturally acquired gastrointestinal nematode infections. Vet Parasitol. 2004;120:91–106.
Van Wyk JA, Reynecke DP. Blueprint for an automated specific decision support system for countering anthelmintic resistance in Haemonchus spp. at farm level. Vet Parasitol. 2011;177:212–23.
Zvinorova PI, Halimani TE, Muchadeyi FC, Matika O, Riggio V, Dzama K. Breeding for resistance to gastrointestinal nematodes – the potential in low-input/output small ruminant production systems. Vet Parasitol. 2016;225:19–28.
Piedrafita D, Raadsma HW, Gonzalez J, Meeusen E. Increased production through parasite control: can ancient breeds of sheep teach us new lessons? Trends Parasitol. 2010;26:568–73.
McManus C, Louvandini H, Paiva SR, Oliveira AA, Azevedo HC, Melo CB. Genetic factors of sheep affecting gastrointestinal parasite infections in the Distrito Federal, Brazil. Vet Parasitol. 2009;166:308–13.
Nunes SF, Ferreira J, Paiva SR, Faria DA, Sousa JER, Soares CEA, et al. Fine genetic structure of Brazilian white Morada Nova hair sheep breed from semi-arid region. Small Rumin Res. 2022;211:106694.
Britton C, Laing R, Mcneilly T, Perez MG, Otto TD, Hildersley KA, et al. New technologies to study helminth development and host-parasite interactions. Int J Parasitol. 2023;53:393–403.
Swann J, Jamshidi N, Lewis NE, Winzeler EA. Systems analysis of host–parasite interactions. Wiley Interdiscip Rev Syst Biol Med. 2015;7:381–400.
Kristensen VN, Lingjærde OC, Russnes HG, Vollan HKM, Frigessi A, Børresen-Dale AL. Principles and methods of integrative genomic analyses in cancer. Nat Rev Cancer. 2014;14:299–313.
Paz EA, Chua EG, Hassan SU, Greeff JC, Palmer DG, Liu S, et al. Bacterial communities in the gastrointestinal tract segments of helminth-resistant and helminth-susceptible sheep. Anim Microbiome. 2022;4:23.
Mamun MAA, Sanderman M, Rayment P, Brook-Carter P, Scholes E, Kasinadhuni N, et al. Variation in gut bacterial composition is associated with Haemonchus contortus parasite infection of sheep. Anim Microbiome. 2020;2:3.
Agüero VCG, Estaban-Blanco CE, Argüello H, Valderas-García E, Andrés S, Balaña-Fouce R, et al. Microbial community in resistant and susceptible Churra sheep infected by Teladorsagia circumcincta. Sci Rep. 2022;12:17620.
Peng S, Yin J, Liu X, Jia B, Chang Z, Lu H, et al. First insights into the microbial diversity in the omasum and reticulum of bovine using Illumina sequencing. J Appl Genet. 2015;56:393–401.
El-Ashram S, Nasr IA, Abouhajer F, El-Kemary M, Huang G, Dinçel G, et al. Microbial community and ovine host response varies with early and late stages of Haemonchus contortus infection. Vet Res Commun. 2017;41:263–77.
Cortés A, Wills J, Su X, Hewitt RE, Roberston J, Scotti R, et al. Infection with the sheep gastrointestinal nematode Teladorsagia circumcincta increases luminal pathobionts. Microbiome. 2020;8:60.
Toscano JHB, Santos IB, Von Haehling MB, Giraldelo LA, Lopes LG, Silva MH, et al. Morada Nova sheep breed: resistant or resilient to Haemonchus contortus infection? Vet Parasitol. 2019;276S:1000019.
Ueno H, Gonçalves PC. Manual para diagnóstico das helmintoses de ruminantes. 4th ed. Tokyo: Japan International Cooperation Agency; 1998. pp. 14–45.
Amadeu RR, Cellon C, Olmstead JW, Garcia AAF, Resende MFR, Muñoz PR. AGHmatrix: R package to construct relationship matrices for autotetraploid and diploid species: a blueberry example. Plant Genome. 2016. 10.3835/plantgenome2016.01.0009.
Bates D, Vazquez AI. pedigreemm: pedigree-based mixed-effects models. R package version 0.3-3. 2014. https://CRAN.R-project.org/package=pedigreemm. Accessed 10 Oct 2023.
R Core Team. R: a language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing; 2022.
Xavier A, Xu S, Muir WM, Rainey KM. NAM: association studies in multiple populations. Bioinformatics. 2015;31:3862–4.
Andrews S. FastQC: a quality control tool for high throughput sequence data. 2010. https://www.bioinformatics.babraham.ac.uk/projects/fastqc/. Accessed 10 Oct 2023.
Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30:2114–20.
Kim D, Paggi JM, Park C, Bennett C, Salzberg SL. Graph-based genome alignment and genotyping with HISAT2 and HISAT-genotype. Nat Biotechnol. 2019;37:907–15.
Danecek P, Bonfield JK, Liddle J, Marshall J, Ohan V, Pollard MO, et al. Twelve years of SAMtools and BCFtools. GigaScience. 2021;10:giab008.
Ewels P, Magnusson M, Lundin S, Käller M. MultiQC: summarize analysis results for multiple tools and samples in a single report. Bioinformatics. 2016;32:3047–8.
Pertea M, Pertea GM, Antonescu CM, Chang TC, Mendell JT, Salzberg SL. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat Biotechnol. 2015;33:290–5.
Van Rossum G, Drake FL. Python 3 reference manual. Scotts Valley: CreateSpace; 2009.
Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-se data with DESeq2. Genome Biol. 2014;15:550.
Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26:139–40.
Fu J, Frazee AC, Collado-Torres L, Jaffe AE, Leek JT. ballgown: flexible, isoform-level differential expression analysis. R package version 2.28.0. 2022. 10.18129/B9.bioc.ballgown.
Stelzer G, Rosen R, Plaschkes I, Zimmerman S, Twik M, Fishilevich S, et al. The GeneCards suite: from gene data mining to disease genome sequence analyses. Curr Protoc Bioinformatics. 2016;54(130):1–33.
Brouard JS, Bissonnette N. Variant calling from RNA-seq data using the GATK joint genotyping workflow. In: Ng CKY, Piscuoglio S, editors. Variant calling. Methods Mol Biol. 2022;2493:205 – 33.
Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013;29:15–21.
Sahraeian SME, Mohiyuddin M, Sebra R, Tilgner H, Afshar PT, Au KF, et al. Gaining comprehensive biological insight into the transcriptome by performing a broad-spectrum RNA-seq analysis. Nat Commun. 2017;8:59.
Broad Institute. Picard toolkit. http://broadinstitute.github.io/picard (2018). Accessed 10 Oct 2023.
Van Der Auwera GA, O'Connor BD. Genomics in the cloud: using Docker, GATK, and WDL in Terra. 1st ed. O'Reilly Media; 2020.
Gao Y, Jiang G, Yang W, Jin W, Gong J, Xu X, et al. Animal-SNPAtlas: a comprehensive SNP database for multiple animals. Nucleic Acids Res. 2023;51:D816–26.
Cingolani P, Platts A, Wang LL, Coon M, Nguyen T, Wang L, et al. A program for annotating and predicting the effects of single nucleotide polymorphisms, SnpEff: SNPs in the genome of Drosophila melanogaster strain w1118, iso-2, iso-3. Fly. 2012;6:80–92.
Cingolani P, Patel VM, Coon M, Nguyen T, Land SJ, Ruden DM, et al. Using Drosophila melanogaster as a model for genotoxic chemical mutational studies with a new program, SnpSift. Front Genet. 2012;3:35.
Gondro C, Porto-Neto LR, Lee SH. snpQC - an R pipeline for quality control of Illumina SNP genotyping array data. Anim Genet. 2014;45:758–61.
Ferdosi MH, Kinghorn BP, Van Der Werf JHJ, Lee SH, Gondro C. hsphase: an R package for pedigree reconstruction, detection of recombination events, phasing and imputations of half-sib family groups. BMC Bioinformatics. 2014;15:172.
Niciura SCM, Benavides MV, Okino CH, Ibelli AM, Minho AP, Esteves SN, et al. Genome-wide association study for Haemonchus contortus resistance in Morada Nova sheep. Pathogens. 2022;11:939.
Genome Data Viewer. Bethesda: National Library of Medicine, National Center for Biotechnology Information. ; 2004. https://www.ncbi.nlm.nih.gov/genome/gdv/browser/genome/?id=GCF_016772045.1. Accessed 10 Oct 2023.
Shabalin AA. Matrix eQTL: ultra fast eQTL analysis via large matrix operations. Bioinformatics. 2012;28:1353–8.
Doyle SR, Tracey A, Laing R, Holroyd N, Bartley D, Bazant W, et al. Genomic and transcriptomic variation defines the chromosome-scale assembly of Haemonchus contortus, a model gastrointestinal worm. Commun Biol. 2020;3:656.
Andrade BGN, Bressani FA, Cuadrat RRC, Tizioto PC, Oliveira PSN, Mourão GB, et al. The structure of microbial populations in Nelore GIT reveals inter-dependency of methanogens in feces and rumen. J Anim Sci Biotechnol. 2020;11:6.
Bolyen E, Rideout JR, Dillon MR, Bokulich NA, Abnet CC, Al-Ghalith GA, et al. Reproducible, interactive, scalable and extensible microbiome data science using QIIME 2. Nat Biotechnol. 2019;37:852–7.
Callahan BJ, Mcmurdie PJ, Rosen MJ, Han AW, Johnson AJ, Holmes SP. DADA2: high-resolution sample inference from Illumina amplicon data. Nat Methods. 2016;13:581–3.
Quast C, Pruesse E, Yilmaz P, Gerken J, Schweer T, Yarza P, et al. The SILVA ribosomal RNA gene database project: improved data processing and web-based tools. Nucleic Acids Res. 2013;41:D590–6.
Ling W, Lu J, Zhao N, Lulla A, Plantinga AM, Fu W, et al. Batch effects removal for microbiome data via conditional quantile regression. Nat Commun. 2022;13:5418.
Russo PST, Ferreira GR, Cardozo LE, Bürger MC, Arias-Carrasco R, Maruyama SR, et al. CEMiTool: a Bioconductor package for performing comprehensive modular co-expression analyses. BMC Bioinformatics. 2018;19:56.
Zhou Y, Zhou B, Pache L, Chang M, Khodabakhshi AJ, Tanaseichuk O, et al. Metascape provides a biologist-oriented resource for the analysis of systems-level datasets. Nat Commun. 2019;10:1523.
Ribeiro ELA, Gonzálex-García E. Indigenous sheep breeds in Brazil: potential role for contributing to the sustainability of production systems. Trop Anim Health Prod. 2016;48:1305–13.
Paim TP, Paiva SR, Toledo NM, Yamaghishi MB, Carneiro PLS, Facó O, et al. Origin and population structure of Brazilian hair sheep breeds. Anim Genet. 2021;52:492–504.
McManus C, Facó O, Shiotsuki L, Rolo JLJP, Peripolli V. Pedigree analysis of Brazilian Morada Nova hair sheep. Small Rumin Res. 2019;170:37–42.
Lacerda TS, Caetano AR, Facó O, Faria DA, McManus CM, Lôbo RN, et al. Single marker assisted selection in Brazilian Morada Nova hair sheep community-based breeding program. Small Rumin Res. 2016;139:15–9.
Laurenson YCSM, Kyriazakis I, Bishop SC. Can we use genetic and genomic approaches to identify candidate animals for target selective treatment. Vet Parasitol. 2013;197:379–83.
Xiang R, Oddy VH, Archibald AL, Vercoe PE, Dalrymple BP. Epithelial, metabolic and innate immunity transcriptomic signatures differentiating the rumen from other sheep and mammalian gastrointestinal tract tissues. PeerJ. 2016;4:e1762.
Li RW, Rinaldi M, Capuco AV. Characterization of the abomasal transcriptome for mechanisms of resistance to gastrointestinal nematodes in cattle. Vet Res. 2011;42:114.
Tang J. 3 - Pepsin A. In: Barrett AJ, Rawlings ND, Woessner JF, editors. Handbook of proteolytic enzymes. 2nd ed. Academic Press; 2004. pp. 19–28.
Wohlkönig A, Huet J, Looze Y, Wintjens R. Structural relationships in the lysozyme superfamily: significant evidence for glycoside hydrolase signature motifs. PLoS ONE. 2010;5:e15388.
Roy HK. COX 3? Am J Gastroenterol. 1999;94:2343.
Dunphy JL, Balic A, Barcham GJ, Horvath AJ, Nash AD, Meeusen ENT. Isolation and characterization of a novel inducible mammalian galectin. J Biol Chem. 2000;275:32106–13.
Rowe A, Gondro C, Emery D, Sangster N. Sequential microarray to identify timing of molecular responses to Haemonchus contortus infection in sheep. Vet Parasitol. 2009;161:76–87.
D’Errico M, Rinaldis E, Blasi MF, Viti V, Falchetti M, Calcagnile A, et al. Genome-wide expression profile of sporadic gastric cancers with microsatellite instability. Eur J Cancer. 2009;45:461–9.
Sui G, Zhou S, Wang J, Canto M, Lee EE, Eshleman JR, et al. Mitochondrial DNA mutations in preneoplastic lesions of the gastrointestinal tract: A biomarker for the early detection of cancer. Mol Cancer. 2006;5:73.
Guo Z, González JF, Hernandez JN, McNeilly TN, Corripio-Miyar Y, Frew D, et al. Possible mechanisms of host resistance to Haemonchus contortus infection in sheep breeds native to the Canary Islands. Sci Rep. 2016;6:26200.
Patra G, Borthakur SK, Das S, Lalliankimi H, Lalrinkima H. Quantification and expression studies of IL-13 gene in response to experimental infection with Haemonchus contortus in Garole and Sahabadi sheep. Vet Parasitol Reg Stud Reports. 2017;8:99–103.
McRae KM, Good B, Hanrahan JP, McCabe MS, Cormican P, Sweeney T, et al. Transcriptional profiling of the ovine abomasal lymph node reveals a role for timing of the immune response in gastrointestinal nematode resistance. Vet Parasitol. 2016;224:96–108.
Liu J, Tan M, Xu X, Shen T, Zhou Z, Hunt PW, et al. From innate to adaptive immunity: abomasal transcriptomic responses of merino sheep to Haemonchus contortus infection. Mol Biochem Parasitol. 2021;246:111424.
Estrada-Reyes ZM, Ogunade IM, Pech-Cervantes AA, Terrill TH. Copy number variant-based genome wide association study reveals immune-related genes associated with parasite resistance in a heritage sheep breed from the United States. Parasite Immunol. 2022;44:e12943.
Aboshady HM, Mandonnet N, Félicité Y, Hira J, Fourcot A, Barbier C, et al. Dynamic transcriptomic changes of goat abomasal mucosa in response to Haemonchus contortus infection. Vet Res. 2020;51:44.
Knight JS, Baird DB, Hein WR, Pernthaner A. The gastrointestinal nematode Trichostrongylus colubriformis down-regulates immune gene expression in migratory cells in afferent lymph. BMC Immunol. 2010;11:51.
Zhang J, Toremurat Z, Liang Y, Cheng J, Sun Z, Huang Y, et al. Study on the association between LRRC8B gene indel and sheep body conformation traits. Genes. 2023;14:356.
Ding N, Tian D, Li X, Zhang Z, Tian F, Liu S, et al. Genetic polymorphisms of IGF1 and IGF1R genes and their effects on growth traits in Hulun Buir hair. Genes. 2022;13:666.
Zhang CI, Chunjie L, Zihu Z, Langman Z, Qianqian C, Zilong C, et al. Analysis on the desert adaptability of indigenous sheep in the southern edge of Taklimakan Desert. Sci Rep. 2022;12:12264.
Wiener P, Robert C, Ahbara A, Salavati M, Abebe A, Kebede A, et al. Whole-genome sequence data suggest environmental adaptation of Ethiopian sheep populations. Genome Biol Evol. 2021;13:evab014.
Eydivandi S, Roudbar MA, Karimi MO, Sahana G. Genomic scans for selective sweeps through haplotype homozygosity and allelic fixation in 14 indigenous sheep breeds from Middle East and South Asia. Sci Rep. 2021;11:2834.
Mengistie D, Edea Z, Tesema TS, Dejene G, Dessie T, Jemal J, et al. Genome-wide signature of positive selection in Ethiopian indigenous and European beef cattle breeds. Anim Gene. 2023;29:200151.
McKnite AM, Perez-Munoz ME, Lu L, Williams EG, Brewer S, Andreux PA, et al. Murine gut microbiota is defined by host genetics and modulates variation of metabolic traits. PLoS ONE. 2012;7:e39191.
Carracelas B, Navajas EA, Bera B, Ciappesoni G. Genome-wide association study of parasite resistance to gastrointestinal nematodes in Corriedale sheep. Genes. 2022;13:1548.
Bhuiyan AA, Li J, Wu Z, Ni P, Adetula AA, Wang H, et al. Exploring the genetic resistance to gastrointestinal nematodes infection in goat using RNA-sequencing. Int J Mol Sci. 2017;18:751.
Araujo RN, Padilha T, Zarlenga D, Sonstegard T, Connor EE, Tassel CV, et al. Use of candidate gene array to delineate gene expression patterns in cattle selected for resistance or susceptibility to intestinal nematodes. Vet Parasitol. 2009;162:106–15.
Liu J, Zhou J, Shao S, Xu X, Li CJ, Li L, et al. Differential responses of abomasal transcriptome to Haemonchus contortus infection between Haemonchus-selected and Trichostrongylus-selected merino sheep. Parasitol Int. 2022;87:102539.
Pemberton JM, Beraldi D, Craig BH, Hopkins J. Digital gene expression analysis of gastrointestinal helminth resistance in Scottish blackface lambs. Mol Ecol. 2011;20:910–9.
Li X, Jiang S, Wang X, Jia B. Intestinal transcriptomes in Kazakh sheep with different haplotypes after experimental Echinococcus granulosus infection. Parasite. 2021;28:14.
Gorfu G, Cirelli KM, Melo MB, Mayer-Barber K, Crown D, Koller BH, et al. Dual role for inflammasome sensors NLRP1 and NLRP3 in murine resistance to Toxoplasma gondii. mBio. 2014;5:e01117–13.
Lechner A, Bohnacker S, Bieren JEV. Macrophage regulation & function in helminth infection. Semin Immunol. 2021;53:101526.
Meulder FV, Coppernolle SV, Borloo J, Rinaldi M, Li RW, Chiers K, et al. Granule exocytosis of granulysin and granzyme B as a potential key mechanism in vaccine-induced immunity in cattle against the nematode Ostertagia ostertagi. Infect Immun. 2013;81:1798–809.
Terra JK, France B, Cote CK, Jenkins A, Bozue JA, Welkos SL, et al. Allelic variation on murine chromosome 11 modifies host inflammatory responses and resistance to Bacillus anthracis. PLoS Pathog. 2011;7:e1002469.
Signer-Hasler H, Burren A, Ammann P, Drögemüller C, Flury C. Runs of homozygosity and signatures of selection: a comparison among eight local Swiss sheep breeds. Anim Genet. 2019;50:512–25.
Widmann P, Reverter A, Weikard R, Suhre K, Hammon HM, Albrecht E, et al. Systems biology analysis merging phenotype, metabolomic and genomic data identifies Non-SMC Condensin I Complex, Subunit G (NCAPG) and cellular maintenance processes as major contributors to genetic variability in bovine feed efficiency. PLoS ONE. 2015;10:e0124574.
An B, Xu L, Xia J, Wang X, Miao J, Chang T, et al. Multiple association analysis of loci and candidate genes that regulate body size at three growth stages in Simmental beef cattle. BMC Genet. 2020;21:32.
King FJM, Visser C, Banga C. Genetic characterization of Mozambican Nguni cattle and their relationship with indigenous populations of South Africa. Livest Sci. 2022;264:105044.
Yang J, Li WR, Lv FH, He SG, Tian SL, Peng WF, et al. Whole-genome sequencing of native sheep provides insights into rapid adaptations to extreme environments. Mol Biol Evol. 2016;33:2576–92.
Mei C, Gui L, Hong J, Raza SHA, Aorigele C, Tian W, et al. Insights into adaptation and growth evolution: a comparative genomics study on two distinct cattle breeds from Northern and Southern China. Mol Ther Nucleic Acids. 2021;23:959–67.
Schubert K, Damink SWMO, Von Bergen M, Schaap FG. Interactions between bile salts, gut microbiota, and hepatic innate immunity. Immunol Rev. 2017;279:23–35.
Boisseau M, Dhorne-Pollet S, Bars-Cortina D, Courtot E, Serreau D, Annonay G, et al. Species interactions, stability, and resilience of the gut microbiota – helminth assemblage in horses. iScience. 2023;26:106044.
The GTEX Consortium. The GTEx Consortium atlas of genetic regulatory effects across human tissues. Science. 2020;369:1318–30.
Gantasala NP, Kumar M, Banakar P, Thakur PK, Rao U. Functional validation of genes in cereal cyst nematode, Heterodera avenae, using siRNA gene silencing. Nematodes Small Grain Cereals. 2015;12:353–6.
Kita K, Nihei C, Tomitsuka E. Parasite mitochondria as drug target: diversity and dynamics changes during the life cycle. Curr Med Chem. 2003;1:2535–48.
Murai M, Miyoshi H. Current topics on inhibitors of respiratory complex I. Biochim Biophys Acta. 2016;1857:884–91.
Nikolaou S, Gasser RB. Prospects for exploring molecular developmental process in Haemonchus contortus. Int J Parasitol. 2006;36:859–68.
Laing R, Kikuchi T, Martinelli A, Tsai IJ, Beech RN, Redman E, et al. The genome and transcriptome of Haemonchus contortus, a key model parasite for drug and vaccine discovery. Genome Biol. 2013;14:R88.
Nakandalage R, Guan LL, Malmuthuge N. Microbial interventions to improve neonatal gut health. Microorganisms. 2023;11:1328.
Iljazovic A, Roy U, Gálvez EJC, Lesker TR, Zhao B, Gronow A, et al. Perturbation of the gut microbiome by Prevotella spp. enhances host susceptibility to mucosal inflammation. Mucosal Immunol. 2021;14:113–24.
Guo P, Zhang K, Ma X, He P. Clostridium species as probiotics: potentials and challenges. J Anim Sci Biotechnol. 2020;11:24.
Hodžić A, Dheilly NM, Cabezas-Cruz A, Berry D. The helminth holobiont: a multidimensional host-parasite-microbiota interaction. Trends Parasitol. 2023;39:91–100.
Hewitson JP, Grainger JR, Maizels RM. Helminth immunoregulation: the role of parasite secreted proteins in modulating host immunity. Mol Biochem Parasitol. 2009;167:1–11.
Cooper D, Eleftherianos I. Parasitic nematode immunomodulatory strategies: recent advances and perspectives. Pathogens. 2016;5:58.
Zhao G, Yan R, Muleke CI, Sun Y, Xu L, Li X. Vaccination of goats with DNA vaccines encoding H11 and IL-2 induces partial protection against Haemonchus contortus infection. Vet J. 2012;191:94–100.
Antonacci R, Massari S, Linguiti G, Jambrenghi AC, Giannico F, Lefranc MP, et al. Evolution of the T-cell receptor (TR) loci in the adaptive immune response: the tale of the TRG locus in mammals. Genes. 2020;11:624.
Fichtner AS, Ravens S, Prinz I. Human γδ TCR repertoires in health and disease. Cells. 2020;9:800.
Calvete O, Reyes J, Valdés-Socin H, Martin P, Marazuela M, Barroso A, et al. Alterations in SLC4A2, SLC26A7 and SLC26A9 drive acid-base imbalance in gastric neuroendocrine tumors and uncover a novel mechanism for a co-occurring polyautoimmune scenario. Cells. 2021;10:3500.
Xu J, Song P, Nakamura S, Miller M, Barone S, Alper SL, et al. Deletion of the chloride transporter slc26a7 causes distal renal tubular acidosis and impairs gastric acid secretion. J Biol Chem. 2009;284:29470–9.
Alonso N, Zappia CD, Cabrera M, Davio CA, Shayo C, Monczor F, et al. Physiological implications of biased signaling at histamine H2 receptors. Front Pharmacol. 2015;6:45.
Howden CW, Hunt RH. Relationship between gastric secretion and infection. Gut. 1987;28:96–107.
Barker IK, Titchen DA. Gastric dysfunction in sheep infected with Trichostrongylus colubriformis, a nematode inhabiting the small intestine. Int J Parasitol. 1982;12:345–56.
Aboshady HM, Choury A, Montout L, Félicité Y, Godard X, Bambou JC. Metagenome reveals caprine abomasal microbiota diversity at early and late stages of Haemonchus contortus infection. Sci Rep. 2023;13:2450.
Ludvik AE, Pusec CM, Priyadarshini M, Angueira AR, Guo C, Lo A, et al. HKDC1 is a novel hexokinase involved in whole-body glucose use. Endocrinology. 2016;157:3452–61.
Azziz G, Giménez M, Carballo C, Espino N, Barlocco N, Batista S. Characterization of the fecal microbiota of Pampa Rocha pigs, a genetic resource endemic to eastern Uruguay. Heliyon. 2023;9:e16643.
Martinez Boggio G, Meynadier A, Daunis-I-Estadella P, Marie-Etancelin C. Compositional analysis of ruminal bacteria from ewes selected for somatic cell score and milk persistency. PLoS ONE. 2021;16:e0254874.
Sbardellati DL, Fischer A, Cox MS, Li W, Kalscheur KF, Suen G. The bovine epimural microbiota displays compositional and structural heterogeneity across different ruminal locations. J Dairy Sci. 2020;103:3636–47.
Wang H, Xia P, Lu Z, Su Y, Zhu W. Metabolome-microbiome responses of growing pigs induced by time-restricted feeding. Front Vet Sci. 2021;8:681202.
Wang X, Zhang Z, Li B, Hao W, Yin W, Ai S, et al. Depicting fecal microbiota characteristic in yak, cattle, yak-cattle hybrid and Tibetan sheep in different eco-regions of Qinghai-Tibetan plateau. Microbiol Spectr. 2022;10:e0002122.
Waters JL, Ley RE. The human gut bacteria Christensenellaceae are widespread, heritable, and associated with health. BMC Biol. 2019;17:83.
Cheng J, Hu J, Geng F, Nie S. Bacteroides utilization for dietary polysaccharides and their beneficial effects on gut health. Food Sci Hum Wellness. 2022;11:1101–10.
Low A, Lee JKY, Gounot JS, Ravikrishnan A, Ding Y, Saw WY, et al. Mutual exclusion of Methanobrevibacter species in the human gut microbiota facilitates directed cultivation of a Candidatus Methanobrevibacter intestini representative. Microbiol Spectr. 2022;10:e00849–22.
Samuel BS, Gordon JI. A humanized gnotobiotic mouse model of host-archaeal-bacterial mutualism. Proc Natl Acad Sci USA. 2006;103:10011–6.
Ohkusa T, Nishikawa Y, Sato N. Gastrointestinal disorders and intestinal bacteria: Advances in research and applications in therapy. Front Med. 2023;9:935676.
Moon CD, Carvalho L, Kirk MR, Mcculloch AF, Kittelmann S, Young W, et al. Effects of long-acting, broad spectra anthelmintic treatments on the rumen microbial community compositions of grazing sheep. Sci Rep. 2021;11:3836.
Chen H, Mozzicafreddo M, Pierella E, Carletti V, Piersanti A, Ali SM, et al. Dissection of the gut microbiota in mothers and children with chronic Trichuris trichiura infection in Pemba Island, Tanzania. Parasit Vectors. 2021;14:62.
Abosse JS, Terefe G, Teshale BM. Comparative study on pathological changes in sheep and goats experimentally infected with Haemonchus contortus. Surg Exp Pathol. 2022;5:14.
Gryaznova M, Dvoretskaya Y, Burakova I, Syromyatnikov M, Popov E, Kokina A, et al. Dynamics of changes in the gut microbiota of healthy mice fed with lactic acid bacteria and bifidobacteria. Microorganisms. 2022;10:1020.
Karpinska A, Ryan D, Germaine K, Dowling D, Forrestal P, Kakouli-Duarte T. Soil microbial and nematode community response to the field application of recycled bio-based fertilisers in Irish grassland. Sustainability. 2021;13:12342.
Xiang H, Fang Y, Tan Z, Zhong R. Haemonchus contortus infection alters gastrointestinal microbial community composition, protein digestion and amino acid allocations in lambs. Front Microbiol. 2022;12:797746.
Li RW, Li W, Sun J, Yu P, Baldwin RL, Urban JF. The effect of helminth infection on the microbial composition and structure of the caprine abomasal microbiome. Sci Rep. 2016;6:20606.
Lee SC, Tang MS, Lim YA, Choy SH, Kurtz ZD, Cox LM, et al. Helminth colonization is associated with increased diversity of the gut microbiota. PLoS Negl Trop Dis. 2014;8:e2880.
Houlden A, Hayes KS, Bancroft AJ, Worthington JJ, Wang P, Grencis RK, et al. Chronic Trichuris muris infection in C57BL/6 mice causes significant changes in host microbiota and metabolome: effects reversed by pathogen clearance. PLoS ONE. 2015;10:e0125945.
Gálvez EJC, Iljazovic A, Amend L, Lesker TR, Renault T, Thiemann S, et al. Distinct polysaccharide utilization determines interspecies competition between intestinal Prevotella spp. Cell Host Microbe. 2020;28:838–52.

No competing interests reported.

Download PDF

Version 1

posted

You are reading this latest preprint version

Multi-omics elucidated parasite-host-microbiota interactions and resistance to Haemonchus contortus in sheep

Status:

Version 1

Abstract

Background:

Results:

Conclusions:

Figures

Background

Methods

Sample collection and phenotypes

RNA-seq

Variant calling and annotation of RNA-seq data

GWAS

eQTL analysis

Analysis of unmapped RNA-seq reads

Microbiota from sheep feces and rumen content

Multi-omics integration

Results

RNA-seq

Variant calling and annotation of RNA-seq data

GWAS

eQTL analysis

Gene co-expression network analysis

Transcript co-expression network analysis

Discussion

Conclusions

Declarations

References

Additional Declarations

Supplementary Files

Status:

Version 1