Genomic scan for identifying candidate genes for paratuberculosis resistance in sheep

Additional keywords: disease susceptibility, GWAS, Ovine SNP50K BeadChip.

Paratuberculosis, or Johne’s, disease is a chronic granulomatous enteritis, which affects ruminants, caused by Mycobacterium avium subsp. Paratuberculosis (MAP). It impacts not only on agriculture economy but also on public health. In fact, although the relationship between MAP and Crohn’s disease in humans is still subject of debate (Balfour Sartor 2005) migrant studies showed that the incidence of the disease, in people moving from a low Crohn’s disease and Johne’s disease incidence area to a high incidence area, subsequently rises to that of the host population (Hermon-Taylor 2009). Moreover Olsen et al. (2009) suggested a possible role of MAP in the inflammation seen in Crohn’s disease. There are no effective treatments against MAP; control programs are very complex and expensive for farmers, and also the vaccine did not show to be sufficiently protective (Reddacliff et al. 2006). However, although MAP is widespread in the environment, the low rate of infected individuals, even in the same flock, suggested the presence of genetic factors influencing disease resistance.

Some functional candidate genes have been associated with MAP in cattle, and with increasing risk of Crohn’s disease in humans. Toll-like receptors (TLR) showed polymorphisms associated with susceptibility to MAP infection (Mucha et al. 2009; Pinedo et al. 2009). Their role in the pathogenesis of Johne’s disease has been identified based on gene expression in sheep (Nalubamba et al. 2008; Taylor et al. 2008). TLR, being members of mammalian pattern-recognition receptors, play a primary role in the recognition of pathogen-associated molecular patterns in bacteria, viruses, protozoa and fungi (Vasselon and Detmers 2002; Kaisho and Akira 2006). Several studies have shown that mutations in TLR genes may reduce the ability to recognise pathogen-associated molecular patterns and hence interfere with innate immune activation. TLR6 may be a potential marker of exposure to MAP and could be used to identify sheep resistant to MAP infection (Plain et al. 2010). Mutations found in TLR1 and 2 may increase the susceptibility to ovine mycobacterial infection (Bhide et al. 2009). NOD2 and CARD15 genes were associated with increasing risk of Crohn’s disease in humans (Hugot et al. 2001; Hugot 2006); however, Sezzi et al. (2007) could not find any significant associations between the polymorphisms detected in the NOD2/CARD15 gene and the disease.

In order to identify genomic regions associated with Johne’s disease paratuberculosis, whole-genome association studies (GWAS) using the Illumina Bovine BeadChip, have been performed in cattle (Settles et al. 2009; Minozzi et al. 2010; Pant et al. 2010; Kirkpatrick et al. 2011; Zanella et al. 2011), but, in sheep, the studies are still scarce.

The recent availability of the Ovine SNP50k BeadChip (Kijas et al. 2012) containing more than 50 000 single nucleotide polymorphism (SNP), distributed along the whole genome, with knowledge of the position on each specific chromosome, recently started to give impulse to the research also in sheep breeding and sheep genomics, with a small number of examples of detected selection sweeps (Moradi et al. 2012; Moioli et al. 2013). With the aim to identify putative candidate genes that play a role in the immune responses, and potentially affect paratuberculosis disease resistance, we genotyped with the Illumina SNP50K BeadChip two groups of dairy sheep of the Sarda breed, selected among the most divergent after diagnosis assessment using enzyme-linked immunosorbent assay (ELISA) test. The genomic regions either surrounding or encoding the markers, that showed a significant effect on the results of the ELISA test, were then checked for the presence of annotated genes either in sheep or other mammals.

Because the complicated nature of the long, slowly developing MAP infections and the lack of good reference tests might likely introduce selection bias in traditional test evaluations (Nielsen and Toft 2008), the trial was performed in a flock of 1000 dairy sheep of the Sarda breed, under the monitoring of the National Veterinary Services Laboratories for several years, that was known to have sheep infected with MAP. Seven-hundred and fifty-nine adult sheep were screened on serum through ELISA. A commercial kit (IDVet, Montpellier, France) was used to diagnose paratuberculosis with the serum antibody ELISA test. Sera were treated according to the manufacturer’s protocol. Animals were considered either serologically positive or negative based on the optical density ratio between the sample and the positive control, as follows:

With the purpose of characterising selective sweep in sheep for resistance to MAP infection, we applied the selective genotyping strategy, which consists of genotyping the extreme divergent animals for a trait, in this case the potential disease resistance. Within the screened 759 sheep, 50 positive and 50 negative sheep were selected, the kit considering positivity S/P values >70. A second ELISA assessment was performed on serum of positive individuals so to confirm the previous results. Average S/P value of the positive sheep was 145.02 ± 35.34; min. = 76.00; max. = 224.00; the corresponding value of the negative sheep was 0.44 ± 16.67; min. = –20.89; max. = 50.78. To further support resistance to MAP, all selected negative sheep were over 6 years of age, under the assumption that older animals were to be considered resistant for not having showed the infection for several consecutive years. The trial was performed in a commercial flock, which keeps no pedigree of the sheep, except the registration number of the ram; breeding groups are composed of ~30 sheep for each ram, and rams are maintained for breeding for about three consecutive years. Therefore, in order to prevent that the between-individuals relatedness affected the results, the sheep of each of the two groups were selected among the daughters of the same rams. The number of the positive-negative daughters of each ram was reported in Table 1.

Table 1.  Distribution of the rams in the two groups of sheep

By selecting highly divergent sheep in the same flock and of the same breed, the linkage disequilibrium (LD) between the SNP and the alleles of the true responsible genes affecting the trait is likely to have been maintained (Goddard and Hayes 2009; Wientjes et al. 2013), this being the required premise for a GWAS aiming to detect the genes influencing the target trait.

DNA of the 100 sheep of the two groups, extracted from blood sample by using Qiagen QIAamp DNA blood mini/midi kit (Qiagen, San Diego, CA, USA), was genotyped using the Ovine SNP50 BeadChip manufactured by Illumina (San Diego, CA, USA). Genotyping was performed by LGS (Cremona, Italy). Raw data were analysed using the GenomeStudio Genotyping Module version 1.7 (Illumina) by applying a no-call threshold of 0.15. Moreover, the markers not satisfying the following filtering parameters were excluded: SNP call rate ≤ 90%; SNP minor allele frequency ≤ 5%; out of Hardy–Weinberg equilibrium at P < 0.01.

Previous works on the genetic diversity of sheep breeds (Kijas et al. 2012; Ciani et al. 2014) made evident that the Ovine SNP50k BeadChip clearly separates the breeds and the geographical breeding areas on the basis of the markers’ allele frequencies. Because the sheep of the present trial belonged to the same breed, the Sarda, and the same farm, a high frequency of markers with similar allele frequencies between the positive and the negative group had to be expected. For this reason, the markers where the frequency of the reference allele was <0.10 different between positive and negative sheep were considered non-informative and were excluded from the statistical analysis.

For each marker, the allelic substitution effect was calculated by regressing the S/P value on the number of copies of the reference allele, using the GLM procedure of SAS software (SAS Institute 2007). To control the chance of any false positive among those markers, the False Discovery Rate (FDR) correction was applied using the Multtest procedure of SAS software (SAS Institute 2007). It was arbitrarily decided that the allelic substitution effect of the potential markers for disease resistance should have FDR P < 0.09. For these markers, the position on the OARv3.1 genome assembly was assessed (http://www.livestockgenomics.csiro.au/sheep/). The presence of an annotated gene close to the marker was then checked in the National Center for Biotechnology Information (NCBI) Ovis aries genome database (http://www.ncbi.nlm.nih.gov/genome/?term=Ovis+aries, verified 23 November 2014), by exploring ~500 kb upstream and downstream the marker of the OARv3.1 region, as suggested by Zare et al. (2014). The explored region was larger than the haplotype block of LD estimated in humans, on one side because the size of these regions are still subject of debate (Reich et al. 2001); on the other side because modern cattle breeds are characterised by large LD blocks, likely caused by evolutionary forces such as genetic drift, admixture, selection and small effective population size, which are common in livestock (Odani et al. 2006; Khatkar et al. 2007).

The position of the marker either in the gene, or in the flanking regions, was assessed by performing a standard nucleotide Basic Local Alignment Search Tool (Altschul et al. 1990) of the OARv3.1 region encompassing the annotated gene with the deposited annotated gene sequence in the NCBI database.

There is increasing evidence that the onset of many human diseases is due to mutations in the intronic regions of genes (Alshatwi et al. 2012). Because such mutations may cause alterations in regulatory regions and splicing process (Doss and Sethumadhavan 2009), to verify whether any of these mutations occurred within the putative binding sites for transcription factors, in silico analysis was performed using the Match software (http://www.gene-regulation.com, verified 23 November 2014).

After data editing, 5091 markers were available for the analysis; their distribution on the ovine chromosomes was reported in Table 2. Thirty-two markers satisfied the significance threshold, having Raw P ≤ 0.0006 and FDR P < 0.9. Fig. 1 represents the Manhattan plot of the F-test values obtained, for each marker, from the GLM procedure of SAS software (SAS Institute 2007). The horizontal line separates the 32 significant markers. For these markers, the allelic substitution effect of the reference allele was reported in Table 3; in this Table, markers were sorted according to the statistical significance of the effect. The 32 markers fell either within the sequence, or in the 5ʹ-3ʹ UTR, of a described gene. The in silico analysis (http://gene-regulation.com, verified 23 November 2014) indicated that two markers were located in putative binding site of transcription factors: OAR8_270360.1 was in the binding site of the myeloid zinc finger 1 (MZF1), reported as one regulator of transcriptional events during hemopoietic development (Hromas et al. 1991); OAR24_23721386.1 was in the binding site of the CCAAT Enhancer Binding Protein β (C/EBPbeta), which is required for the gene transcription of the nuclear factor of activated T-cells (Yang and Chow 2003).

Table 2.  Distribution on the ovine chromosomes of the markers used in the present study

Fig. 1.  Manhattan plot of the F-test values obtained, for each marker, from the GLM procedure. The horizontal line separates the 32 significant markers.

Table 3.  Allelic substitution effect of the reference allele, F-value and probability values for the 32 markers with FDR P < 0.09, with the position of the marker on OARv3.1

The search of the annotated genes was performed for each of the 32 markers of Table 3; the genes and the position in the gene of these markers as well the accession number of an annotated mRNA of the gene, when available, were reported in Table 4. In this Table the markers were sorted according to the chromosome and the position on the chromosome, to allow for a simplified evaluation of a potential LD of some of them, which might corroborate the effect of a putative candidate gene on the analysed performance. The Table reports also the gene function, as in the databases http://ghr.nlm.nih.gov/gene (verified 23 November 2014) and http://www.uniprot.org/uniprot (verified 23 November 2014); according to their function, the genes were classified in a pie graph based on the gene ontology (Fig. 2).

Table 4.  Gene symbol, name and function, accession number and position in the gene of the markers of Table 3, ordered by chromosomes
n.a., no deposited mRNA sequence

Fig. 2.  Major ontological and functional classification of the genes in proximity of the markers significantly affecting paratuberculosis positivity (FDR P < 0.09).

With the purpose to detect genes that might play a role in disease resistance, this study was performed in a dairy flock that was known to have sheep infected with MAP; sheep were screened on serum through ELISA. Serology provides a cost-effective alternative to organism detection-based diagnostic methods for bovine paratuberculosis, although assay specificity is less than that for faecal culture (Collins 2002). Sensitivity and specificity estimates from studies of serum antibody ELISA for detection of affected, infectious and infected sheep were reported by Nielsen and Toft (2008) and Gumber et al. (2006). Both authors refer values for specificity ranging from 95% to 100% and suggest that appropriate animal sampling for flock diagnoses can be performed to meet specific levels of confidence. To achieve the best possible conditions, the selective genotyping strategy was here applied, by selecting two groups of 50 positive and 50 negative sheep, out of 759 screened animals, representing the highest divergent animals for the target trait. On the sheep of the two divergent groups, we performed a GWAS with the Illumina SNP50K BeadChip with the purpose to detect genes that might have a role in disease resistance.

In order to discern whether the genes in proximity of the 32 markers were either directly or indirectly involved in Johne’s disease, we examined the literature that included experimental trials performed with non-infected versus Mycobacterium bovis-infected cattle. These works are based on two different approaches, the first is a gene expression screenings with specific DNA microarrays containing immunity genes; the second is the identification of loci associated with tolerance to the disease through a GWAS. Because single genetic markers of susceptibility genes with small effects have limited predictive value, a genetic risk profile requires extensive knowledge of gene–gene interactions (Janssens et al. 2008); therefore, to expedite the search for genes associated with polygenic diseases, the complementary approaches of whole-genome scans and microarray gene profiling is used to identify and validate clusters of relevant genes (Gibbons et al. 2004). Coussens et al. (2005) followed the first approach and evaluated gene-expression profiles of peripheral blood mononuclear cells from two groups of Holstein cows, consisting of four animals naturally infected with MAP and four healthy uninfected controls from a commercial herd. These authors performed microarray hybridisations using a commercial microarray hybridisation station (BOTL-3 bovine cDNA microarrays, CAFG – Michigan State University) containing expressed sequence tag clones representing genes known to function in immune response, and evaluated function and ontology of the differentially expressed genes; they showed that 40% of the genes could be classified into the following ontological and functional categories: cytokines and receptors, signal transduction, nuclear and transcription factors, cell cycle control, signal transduction, cytoskeletal; a further 30% of the differentially expressed genes were involved in apoptotic processes, cell growth and migration, adhesion, inhibitors and metabolism. Similarly to what was reported in the work by Coussens et al. (2005), in Fig. 2 the here detected genes were classified in a pie graph based on the gene ontology. The majority of the genes, also in the present study, fell in the category of the activator of transcription, the growth factor activators and the signal transduction (50%; Fig. 2); another important category (37%) included the genes involved in metabolic processes and ion transport.

Also Meade et al. (2007), Magee et al. (2012) and Nalpas et al. (2013) have performed gene expression screenings, with DNA microarrays containing immunity genes, on peripheral blood mononuclear cells from non-infected versus M. bovis-infected cattle. Meade et al. (2007) reported the lower relative expression of key innate immune genes in the infected animals. Magee et al. (2012) revealed a differential expression for genes involved in: (1) the inflammatory response; (2) cell signalling pathways, including TLR and intracellular pathogen recognition receptors; and (3) apoptosis. Nalpas et al. (2013) demonstrated the differential expression of immune, apoptotic and cell signalling genes in the non-infected versus infected animals. The results of the reported studies actually showed that M. bovis infection is associated with the repression of host immunity gene expression.

The second approach to identify loci associated with tolerance to the disease is the GWAS. Similarly to the present study, Zanella et al. (2011) performed a GWAS in Holstein cattle with the Illumina Bovine SNP50 BeadChip, and found that only one gene, the GNA12 gene (guanine nucleotide binding protein α 12) could be considered a positional candidate associated to tolerance to Johne’s disease, so that the previous findings, emphasising the important role of CARD15/NOD2 and TLR4 genes, were not corroborated. Zare et al. (2014) performed a case-control GWAS to identify genomic regions underlying susceptibility to MAP infection in Jersey cattle, but using a less dense SNP panel (7K SNP). They reported that the significant markers were in proximity of the following candidate genes: major histocompatibility complex, TCF19, FAT10, HIVEP1, CCDC17, ZNF684, UBE2 L3, UBE2K, FAM5C, FAM109A. These genes included transcription regulators and ubiquitins. Kirkpatrick et al. (2011) performed a GWAS with the SNP50 BeadChip on MAP-infected Holstein cows and compared SNP allele frequencies with specific reference populations, their aim being the identification of a set of SNP that could be used as a predictor for Holstein cattle susceptibility to MAP infection. Of the genes detected though GWAS in humans and mice, inferred as potential candidates of Crohn’s disease and colitis, Kirkpatrick et al. (2011) found that only the PTGER4 gene was in proximity of one of the SNP of the proposed set.

In the present study, 5 of the 30 detected genes had been previously reported to play a direct role in the immune system: SEMA3D, CD109, PCP4, PRDM2 and ITFG2. SEMA3D is a member of the semaphorins that have been demonstrated to be important for immune response by reducing the activation of T-cells through its cell-surface receptors, including members of the neurophilin and plexin families (Moretti et al. 2006). Emerging evidence also indicated that additional semaphorins and related molecules seem to function in the reciprocal stimulation of T-cells and antigen-presenting cells.

The CD109 gene, belonging to the complement gene family, has been studied in depth in humans, where the differential expression of the gene was reported for different human carcinoma, compared with the control (Hashimoto et al. 2004; Zhang et al. 2005). A direct effect of the marker in proximity of this gene, on the differential expression of the CD109 gene in resistant versus non-resistant sheep, might furthermore be inferred, falling marker OAR8_270360.1 in the binding site of the MZF1 transcription factor.

The PCP4 gene is involved in the activation of B-cells: in the murine B-cell co-receptor complex, Jacobson et al. (2009) demonstrated an increased expression of PCP4, in wild type mice, compared with complement-deficient animals; the same authors correlated the increased expression of PCP4 with B-cell maturation into end stage phenotypes.

The PRDM2 gene belongs to the positive regulatory domain I containing genes, with zinc fingers (PRDM) family, which has recently invoked considerable interest as it has been implicated in fundamental aspects of cellular differentiation and exhibits expanding ties to human diseases (Fog et al. 2012). It is interesting to note that the PRMD1 gene was considered by Minozzi et al. (2010) a strong functional candidate of MAP resistance; these authors have, in fact, performed a GWAS in Holstein cows, using the Illumina Bovine SNP50 BeadChip, and have identified SNP on chromosomes 12, 8, 9, 11, and 27 with significant association with MAP antibody response, but reported that the only gene that could be considered a strong functional candidate was PRMD1, a transcription repressor that acts on the β interferon gene expression and affects the maturation of B-lymphocytes into antibody secreting cells.

The ITFG2 gene, belonging to the integrin α family, was reported to be involved in cell-mediated responses by Zhuang et al. (2008).

With the exception of the PRDM members, proposed as potential candidate of MAP resistance in sheep, in the present study, and in cattle, by Minozzi et al. (2010), each GWAS has identified different genes. It is therefore evident that the complexity of the immune system cannot be summarised in few genes with major effect.

On the contrary, in this study it was evident that the genes responsible for basic cellular processes play a role in MAP resistance, corroborating the findings of Coussens et al. (2005), Meade et al. (2007) and Nalpas et al. (2013).

Most of the existing studies on immune response to MAP infection, both in humans and in livestock, used the candidate gene approach, which is based on a pre-selection of the genes that are hypothesised to influence natural resistance, and following analysis of these genes, either of the polymorphisms or the gene expression, in the positive and control groups. This approach precludes the detection of new genes, potentially involved in the immune response, simply because their role was not previously observed or hypothesised. For this reason, in the present study we aspired to identify regions associated with resistance to paratuberculosis in sheep. Although the present study should be considered preliminary, because the GWAS was performed on a limited number of individuals, it describes the first genome-wide characterisation of selective sweeps in sheep addressing disease resistance. Although animal health has been identified as a priority, breeding for disease resistance is still a scarcely explored area and there is a clear need of research, particularly for sheep.

In this study, 32 genomic regions were made evident for encoding polymorphic markers with significant effect on the S/P value indicating paratuberculosis positivity. These regions harboured already described ovine genes; for five of these genes, a direct involvement in the immune response was reported in the literature.

The explorative value of this study might be corroborated by using the obtained results, on one side, to deeply analyse each gene, so to detect the polymorphism that alter gene expression, likely being in LD with the marker. On the other side, specific quantitative real-time PCR trials of gene expression for these genes could be designed on experimental resistant versus non-resistant populations, so to find markers to be used in the selection for disease resistance.