Genome-wide association study to identify genetic loci associated with gastrointestinal nematode resistance in Katahdin sheep.

Resistance to gastrointestinal nematodes has previously been shown to be a moderately heritable trait in some breeds of sheep, but the mechanisms of resistance are not well understood. Selection for resistance currently relies upon faecal egg counts (FEC), blood packed cell volumes and FAMACHA visual indicator scores of anaemia. Identifying genomic markers associated with disease resistance would potentially improve the selection process and provide a more reliable means of classifying and understanding the biology behind resistant and susceptible sheep. A GWAS was conducted to identify possible genetic loci associated with resistance to Haemonchus contortus in Katahdin sheep. Forty animals were selected from the top and bottom 10% of estimated breeding values for FEC from a total pool of 641 sires and ram lambs. Samples were genotyped using Applied Biosystems™ Axiom™ Ovine Genotyping Array (50K) consisting of 51 572 SNPs. Following quality control, 46 268 SNPs were included in subsequent analyses. Analyses were conducted using a linear regression model in plink v1.90 and a single-locus mixed model in snp and variation suite. Genome-wide significance was determined by a Bonferroni correction for multiple testing. Using linear regression, loci on chromosomes 2, 3, 16, 23 and 24 were significantly associated at the genome level with FEC estimated breeding values, and we identified a region on chromosome 2 that was significant using both statistical analyses. We suggest a potential role for the gene DIS3L2 for gastrointestinal nematode resistance in Katahdin sheep, although further research is needed to validate these findings.

Of the gastrointestinal nematodes (GIN) that infect small ruminants, the abomasal parasite Haemonchus contortus arguably represents the greatest economic concern. H. contortus possesses the highest prevalence of anthelmintic resistance and is the most abundant GIN (Fleming et al. 2006). Female worms can produce 5000–15 000 eggs daily, resulting in rapid accumulation of infective larvae on pastures (Emery et al. 2016). Fourth‐stage larvae and adult nematodes consume blood through the hosts’ abomasal mucosa (Emery et al. 2016). Individual worms can remove up to 30 µl of blood per day, which can cause fatal anaemia in young or immunocompromised animals (Zajac 2006; Emery et al. 2016). H. contortus infections account for significant production losses, and concerns regarding treatment costs and anthelmintic resistance have encouraged the development of other methods of GIN control (Saddiqi et al. 2011). Current selection strategies use phenotypic markers such as faecal egg count (FEC), antibody assays, packed cell volume (PCV) and FAMACHA visual anaemia score to identify animals that are more resistant to H. contortus (Burke & Miller, 2008; Shaw et al. 2012; Aguerre et al. 2018). Phenotypic selection can be labour intensive and costly, and accuracy of selection depends upon many factors that may be difficult to control, such as variations in natural helminth infection and environmental load depending upon the season (Woolaston & Baker, 1996; Uriarte et al. 2003; Jackson & Miller, 2006).

It is hypothesised that inherited GIN resistance is polygenic and related to the immune system (Saddiqi et al. 2011; Atlija et al. 2016; Aguerre et al. 2018). Protection against GIN has been associated with the T helper 2 (Th2) immune response (Moncada et al. 2003) characterised by secretion of interleukin‐10 (IL‐10) and other cytokines which promote recruitment of eosinophils, basophils and mast cells to control infection and mediate helminth expulsion (Begley & Nicola, 1999; Hussaarts et al. 2014; McRae et al. 2015). The Th2 response has also been described as a mediator for acute wound healing during helminth infection (Chen et al. 2012). IL‐10 expression promotes maintenance of the Th2 response by inhibiting the development of Th1 cells and thereby preventing the expression of proinflammatory cytokines such as interferon gamma (IFN‐γ) (Bigley & Nicola 1999; Couper et al. 2008).

Katahdin sheep are an economically important breed in the United States, ranking within the top six of registered breeds for the last several years (Morgan 2016). The Katahdin breed was developed through crosses of British wool breeds with the St Croix hair sheep (Wildeus 1997). The St Croix breed has been noted for its immune response and relative resistance to GIN (Burke & Miller 2004; Bowdridge et al. 2015). Some mature Katahdin ewes have been shown to have GIN resistance traits similar to St Croix upon natural infection (Burke & Miller 2002). Recent work estimated the heritability (h 2) of FEC in Katahdin lambs at 60 (h = 0.18–0.26) and 120 (h 2 = 0.23–0.46) days of age (Ngere et al. 2018), and periparturient ewes at lambing (h 2 = 0.29–0.41) and 30 days postpartum (h 2 = 0.17–0.31) (Notter et al. 2018). These generally moderate heritabilities for FEC suggest that parasite resistance can be improved through genetic selection.

In the current study animals were selected from a pool of 641 Katahdin ram lambs and sires enrolled in the National Sheep Improvement Program database. Selection was based upon high and low EBV for weaning FEC. Katahdin breeders collected stool samples from each animal’s rectum and FEC was quantified by a certified parasitology laboratory (LSU & V Tech Vet School). FEC data were submitted to the National Sheep Improvement Program for EBV prediction by LAMBPLAN (MLA 2004) to reduce variability owing to non‐genetic factors (Brown & Tier 2003; Ferguson 2016). To generate a more diverse sample set for genotyping, less related animals within each EBV category were chosen using a measure of genetic relatedness based on pedigree (Lewis et al. 2005; Kuehn et al. 2008) (Fig. S1). Using a multivariate approach, more related individuals were clustered. Individuals from discrete clusters were then selected. In total, 33 ram lambs and seven sires from eight different US farms were selected for inclusion in the GWAS, with 16 of these animals possessing high FEC EBV and 24 animals possessing low FEC EBV. The two categories captured approximately the top and bottom 10% of animals based on FEC EBV, a phenotypic indicator of GIN resistance.

DNA was extracted from blood samples using the phenol–chloroform method as previously described (Sambrook et al. 1989). The Applied Biosystems™ Axiom™ Ovine Genotyping Array (50K), which included 51 572 SNPs, was used for genotyping against the oar v4.0 reference genome assembly. Quality control was performed in plink v1.07, first excluding non‐autosomal SNPs and SNPs with a call rate less than 0.90, and then a MAF less than 0.01. Following quality control, 4885 SNPs were identified as duplicate markers. These duplicate SNPs were filtered by Fisher’s linear discriminate (FLD) genotype cluster quality score (Johnson & Wichern 2002). The marker with the best FLD genotype cluster score for each of the duplicated SNPs was retained. Following quality control and FLD filtering, 46 268 SNPs were included in analyses.

Genome‐wide association was first conducted for FEC EBV using a linear regression (LR) model using plink v1.90 software (Chang et al. 2015). The package 'qqman' in r version 3.5.1 was used to create the Manhattan plot (Turner 2014) (Fig. 1a).A second analysis was performed using a single‐locus mixed‐model (SLMM) through snp and variation suite™ version 8.7.2 (Golden Helix, Inc., http://www.goldenhelix.com) (Fig. 1b). The SLMM fits a kinship matrix to correct for cryptic relatedness as a random effect and offers more stringent control of the false discovery rate (Segura et al. 2012; Brzyski et al. 2017). Genome‐wide significance was defined by a Bonferroni correction for multiple testing [P ≤ −log10 (1.25−6)]. Genome‐wide suggestive was defined by [P ≤ −log10 (1 × 10−5)] (Fig. 1a,b).

Figure 1

Significant SNPs identified through linear regression and single‐locus mixed‐model GWASs. Genome‐wide significance is defined by a Bonferroni correction for multiple testing [P ≤ −log10 (1.25 × 10−6); red line] and genome‐wide suggestive is defined by [P ≤ −log10 (1 × 10−5); blue line]. (a) Left: Manhattan plot for the linear regression model; right: quantile–quantile plot for the linear regression model. (The QQ plot shows an early departure of observed from expected P‐values, suggesting that this method does not sufficiently account for P‐value inflation and may indicate the presence of type I error.) (b) Left: Manhattan plot for the single locus mixed‐model; right: quantile–quantile plot for the single locus mixed‐model. (c) Table displaying results for both models. LR, Linear regression; SLMM, single‐locus mixed‐model; MAF, minor allele frequency. Markers identified in previous studies are denoted as: *Periasamy et al. (2014), ⁺Atlija et al. (2016), #Sallé et al. (2012).

Using LR analysis, loci on chromosomes 2, 3, 16, 23 and 24 were significantly associated at the genome level with FEC EBV (Fig. 1a). Interestingly, a locus on chromosome 3 was within the gene complement C3a receptor 1 (C3AR1) and the locus on chromosome 16 was 87 kb upstream of the gene integrin subunit alpha 2 (ITGA2). C3AR1 has been reported to be differentially expressed in more susceptible vs resistant sheep (Ahmed et al. 2015) and has been associated with the Th1 immune response (Ghannam et al. 2014). Integrin subunit alpha genes have been found to be upregulated in more resistant animals under GIN infection (Zhang et al. 2019). The most significant SNP was located 2507 bp upstream of the gene DIS3 like 3′–5′ exoribonuclease 2 (DIS3L2). The MAF of this SNP was over‐represented in resistant (0.479) in comparison with susceptible (0.094) sheep (Fig. 1c). An additional SNP was located within the second intron of DIS3L2; this SNP reached significance in the LR model and was most significant using the SLMM model (Fig. 1b). Conversely, the MAF of this SNP was over‐represented in susceptible (0.719) in comparison with resistant (0.208) sheep (Fig. 1c).

DIS3L2 is the cytoplasmic exoribonuclease required for the decay of uridylated pre‐let‐7 and repression of let‐7 (lethal‐7) miRNAs in the Lin28A pathway (Chang et al. 2013; Ustianenko et al. 2013). Let‐7 has been associated with the immune response to parasite infection through the direct regulation of toll‐like receptor 4 expression, and researchers found that in vitro suppression of let‐7i miRNAs in human cholangiocyte cells corresponded with decreased Cryptosporidium parvum parasite burden (Chen et al. 2007). Let‐7 directly affects IL‐10 expression through binding to the IL10 3′UTR; both let‐7a and let‐7d family members have been associated with repression of IL‐10 (Schulte et al. 2011; Swaminathan et al. 2012). An absence of IL‐10 has been shown to result in increased levels of IFN‐γ and delayed expulsion of Trichinella spiralis in mice (Helmby & Grencis, 2003). Upregulation of IFN‐γ has been described in the abomasum and abomasal lymph nodes of sheep that were classified as susceptible to H. contortus infection (Zaros et al. 2014). Reduction of IFN‐γ expression may enhance the immune response to GIN, as this would favour the Th2 cell subset and antibody‐associated immune mechanisms (Coltman et al. 2001).

This preliminary study suggests that the gene DIS3L2 may have a role in GIN resistance in Katahdin sheep. We theorise that polymorphisms within or associated with DIS3L2 moderate its function to allow preferential degradation of the oligoU tail but not degradation of pre‐let‐7 itself. Instead, pre‐let‐7 may enter the maturation pathway following removal of the oligoU tail (Fig. 2). Considering the small sample size and uncorrected inflation factor in the LR model (λ = 1.51) (Fig. 1a), additional work is needed to validate these preliminary findings. Validation may be accomplished through increasing the sample size or expanding the study to include other populations of hair sheep. Future research is required to understand the immune mechanisms that differentiate GIN resistance from susceptibility.

Figure 2

Theorised DIS3L2 and let‐7 pathway mediating gastrointestinal nematode (GIN) resistance. In the presence of LIN28A and TUT4/TUT7 the pre‐let‐7 miRNA is uridylated and marked for degradation by DIS3L2. We propose that SNPs associated with DIS3L2 alter its function such that uridylation of pre‐let‐7 becomes a reversible event. DIS3L2 removes the oligoU tail in this model to allow pre‐let‐7 to continue into the maturation pathway with TUT4/TUT7 and Dicer. Mature let‐7 miRNAs bind to the IL10 3′UTR and prevent/reduce IL‐10 expression, ultimately allowing the Th1 immune response to be favored over the Th2, resulting in decreased resistance to GIN infection.

Conflict of interest

The authors have no conflict of interest to declare.

Figure S1 Eigenvector plot showing animal relatedness and phenotypic distribution.

Click here for additional data file.