An analysis of effects of heterozygosity in dairy cattle for bovine tuberculosis resistance.

Genetic selection of cattle more resistant to bovine tuberculosis (bTB) may offer a complementary control strategy. Hypothesising underlying non-additive genetic variation, we present an approach using genome-wide high density markers to identify genomic loci with dominance effects on bTB resistance and to test previously published regions with heterozygote advantage in bTB. Our data comprised 1151 Holstein-Friesian cows from Northern Ireland, confirmed bTB cases and controls, genotyped with the 700K Illumina BeadChip. Genome-wide markers were tested for associations between heterozygosity and bTB status using marker-based relationships. Results were tested for robustness against genetic structure, and the genotypic frequencies of a significant locus were tested for departures from Hardy-Weinberg equilibrium. Genomic regions identified in our study and in previous publications were tested for dominance effects. Genotypic effects were estimated through ASReml mixed models. A SNP (rs43032684) on chromosome 6 was significant at the chromosome-wide level, explaining 1.7% of the phenotypic variance. In the controls, there were fewer heterozygotes for rs43032684 (P < 0.01) with the genotypic values suggesting that heterozygosity confers a heterozygote disadvantage. The region surrounding rs43032684 had a significant dominance effect (P < 0.01). SNP rs43032684 resides within a pseudogene with a parental gene involved in macrophage response to infection and within a copy-number-variation region previously associated with nematode resistance. No dominance effect was found for the region on chromosome 11, as indicated by a previous candidate region bTB study. These findings require further validation with large-scale data.

Introduction

Bovine tuberculosis (bTB) is an infectious zoonotic disease caused by Mycobacterium bovis and has devastating consequences for the UK cattle industry. Previous studies have demonstrated genetic variation in resistance to bTB (Bermingham et al. 2009; Brotherstone et al. 2010) and have shown the feasibility of genomic selection for bTB resistance using information from markers dispersed throughout the genome (Tsairidou et al. 2014). An alternative approach for selecting more resistant individuals is through identifying individual QTL that explain significant amounts of the genetic variance and exploiting this knowledge in marker‐assisted selection, which was examined for bTB resistance by Bermingham et al. (2014). Selection for disease resistance based on single QTL has previously been implemented successfully in aquaculture and livestock species: e.g. resistance to infectious pancreatic necrosis in Atlantic salmon (Houston et al. 2010) and resistance to Escherichia coli in pigs (Jørgensen et al. 2003; National Committee for Pig Production 2004).

Such studies typically assume that the QTL have additive effects. However, some of the phenotypic variation in disease resistance could also be explained by QTL with dominance effects. The increasing and decreasing expression of dominance is commonly assumed to underlie hybrid vigour and inbreeding depression respectively (Falconer & Mackay 1997; Oldenbroek 2007). Assuming a biallelic locus, over‐ or under‐dominance occurs when the heterozygote genotype is superior or inferior (respectively) to both homozygotes, and partial dominance occurs when the heterozygote is intermediate. For fitness, with over‐dominance, the selection is balancing and both variants would be expected to be maintained in a large population, whereas in other cases the population would be expected to move towards homozygosity for one of the alleles although it has reduced fitness. Over‐dominance is more likely than partial dominance to influence the decision on whether or not to fix a QTL with a favourable additive effect. Over‐dominance has been documented in disease traits, e.g. sickle cell anaemia in humans (Allison 1964) and crooked tail syndrome in Belgian Blue cattle (Fasquelle et al. 2009). Under‐dominance has also been documented, e.g. for the Rhesus blood group system in humans (Haldane 1941) and in establishing chromosomal rearrangements (White 1973) but, to our knowledge, not for resistance to infectious disease.

Therefore, the objective of this study was to detect non‐additive genetic variation associated with bTB resistance by searching for individual QTL exhibiting dominance effects. For this purpose, we used a genome‐wide association study (GWAS). This method has been widely used in human and animal studies to identify associations between genetic polymorphisms and traits through the use of large‐scale SNP data (Andersson 2009). However, the standard GWAS often fails to explore non‐additive genetic variation, and we present an adaptation directed towards this aim. A previous publication provided suggestive evidence of over‐dominance protecting heterozygotes from culling as a result of the bTB diagnostic SICCT test (Amos et al. 2013), which would be an important consideration shaping bTB eradication policy if it were confirmed. Therefore, an additional objective was to test these previous findings on a much larger dataset. Our findings are compared with results from standard genome scans, and the genomic regions identified are explored for candidate genes.

Materials and methods

Population

The dataset comprised 1151 Holstein‐Friesian cows of unknown pedigree from 165 herds in Northern Ireland; 592 of these were confirmed bTB cases and the remaining 559 were controls with multiple negative tuberculin test results. The available information included the cow's date of birth, testing dates and reasons for testing (e.g. routine surveillance testing and control tests within breakdowns); see Tsairidou et al. (2014) for further details. All cases were positive for the single intradermal comparative cervical test (SICCT) and had bTB lesions confirmed by carcass inspection, culture and molecular tests (Bermingham et al. 2014). Controls were age‐matched with cases and selected from herds with higher disease prevalence to increase their probability of exposure to the pathogen (Bishop 2012).

Quality control

All individuals were genotyped with the 700K BovineHD Illumina BeadChip. Quality control was conducted within the genabel package (r/2.15.2) (Aulchenko et al. 2007b). Markers with less than 0.05% minor allele frequencies and less than 0.90% call rates were excluded. Animals were excluded when having less than 0.30% call rate, and one putative cow was discarded as it was a male (odds > 1000). This left 573 102 markers from 1150 individuals for calculating identity‐by‐state (IBS) relationships (see below). For the SNP of interest identified in the Results, the genotype call graph was examined for ambiguities, but none were found. The outcomes are given in Fig. S1, which shows that genotypic classes were found to form three distinct clusters, indicative of reliable genotyping.

Genetic structure

The underlying genetic structure was examined using classical multidimensional scaling (CMDS) based upon pairwise distances between individuals calculated from the genome‐wide IBS information using the r/2.15.2 function ‘cmdscale’ (R Development Core Team 2008). CMDS revealed a secondary distinct cluster of 40 individuals, 39 of these from a single herd that included crossbreds, as had been observed previously (Bermingham et al. 2014). Therefore to address this, the analyses described below were repeated with and without this cluster.

Heterozygote advantage GWAS

An adaptation of the GWAS was developed to identify non‐additive genetic variation. Firstly, genotypes were recoded as a binary heterozygosity indicator: so d  = 1 if individual i was a heterozygote at locus j and d  = 0 otherwise. The ‘polygenic’ function of genabel was used to obtain residuals under a polygenic model (Aulchenko et al. 2007a,b) for use with the ‘mmscore’ function to test for associations between the d and the bTB phenotypes. The polygenic model was: where y is the binary bTB status (0: control; 1: case); m is the overall mean and 1 a vector of 1's; b is the vector of fixed effects with incidence matrix X, comprising the age of the individual as a covariate, year of testing, season of testing, the reason for testing and the two principal components calculated from the CMDS analysis on the genomic relationship matrix G (where data from the secondary cluster was included); a is the vector of polygenic breeding values with a assumed to follow a multivariate normal distribution, a~MVN (0, σ a 2 G), where σ a 2 is the additive genetic variance; and e is the residual error, assumed to be distributed MVN (0, σ e 2 I).

The model was fitted within genabel, which calculates G following Amin et al. (2007): where g is the genomic relationship between animals i and j; n is the number of loci used for estimating relationships; x is the count of alternative alleles (0, 1 or 2) of individual i at SNP locus k, coded as 0, 1 and 2, where the reference allele is arbitrarily chosen; p is the frequency of the reference allele in the data; H is the expected heterozygosity at locus k assuming Hardy‐Weinberg equilibrium (HWE); and H is the observed heterozygosity for animal i (0 or 1).

The model applied to each locus k by the ‘mmscore’ function was: where m is the overall mean with 1 a vector of 1's; d is the vector of heterozygosity indicators for locus k and b is the associated regression coefficient; and e * are residuals assumed to be distributed MVN (0, νe 2 I). To provide a check on consistency with the additive screen published by Bermingham et al. (2014), d was replaced by x in Model (Model 2), where x is the vector of counts for locus k described above.

A detailed analysis of departures from HWE was carried out for significant SNPs. The models and results are presented in Appendix S1.

Significance thresholds applied

The genome‐wide degree of inflation (λ) for Model (Model 2) using the heterozygosity indicators d was 1.00 (SE = 1.02 × 10−5), indicating no serious inflation due to hidden structure in the data, although accounted for in significances quoted. Bonferroni corrections were applied to obtain genome‐ and chromosome‐wide significance thresholds, derived as −log10(0.05/n) and −log10(0.05/n ) respectively (where n is the number of SNPs on chromosome c). Particular attention was given to regions identified by Driscoll et al. (2011) and associated with a protective effect against reacting to the bTB diagnostic SICCT test (Amos et al. 2013), namely microsatellites INRA111 (BTA11: 40 311 694–40 311 817 bp) and BMS2753 (BTA9: 76 800 661–76 800 769 bp). Here, the 20 SNPs on each side of the two microsatellites were assessed using a significance threshold of −log10(0.05/80) = 3.2.

Regional and genomic additive and dominance effects

The significance of the SNPs identified using Model (Model 2) and those flanking the microsatellites identified by Driscoll et al. (2011) was further tested using mixed linear models with both genome‐wide and regional additive and dominance effects fitted as random effects. The regional effects were defined by windows of 40 SNPs, 20 on each side of the locus of interest. The full model was: where y is the binary bTB status (0: control, 1: case); m is the overall mean; β is the vector of fixed effects as described in (1) with incidence matrix X; a is the vector of additive polygenic effects with a~MVN(0, σ a 2 G); d is the vector of polygenic dominance effects with d~MVN(0, σ d 2 D); α is the vector of regional additive effects with α~MVN(0, va 2 G ), δ is the vector of regional dominance effects with δ~MVN(0, vd 2 D ); and e is the vector of residual errors with e ~ MVN (0, I σ e 2). The incidence matrices for all the random effects was I in this data. The (co)variance matrix G was as described above for Model (Model 1), and G was obtained using the same procedure using only the 40 SNPs defining the window. The (co)variance matrices of dominance effects D and D were calculated as: where d is the genomic relationship between individuals i and j; m is the number of SNPs included, m = n for D and 40 for D ; z is the heterozygosity indicator for locus k of individual i; p is the frequency of the reference allele of locus k in the data; and q is frequency of the alternative allele. Likelihood ratio tests were conducted to identify significant regional effects with the genome‐wide terms a and d assumed to be nuisance factors and included in all models. In addition, alternative models were tested for the region containing the SNP with heterozygote disadvantage after excluding this SNP from the calculation of the local matrices.

Results

Additive GWAS model

Results from the standard GWAS are shown in Appendix S2. The SNPs identified in the present study were consistent with findings previously published by Bermingham et al. (2014) using the same data but different statistical procedures.

GWAS for heterozygote advantage

The GWAS for heterozygote advantage identified a SNP on BTA6 (rs43032684) which was significant at the chromosome‐wide level with −log10 P = 6.21 (Table 1, Fig. 1). The rs43032684 locus on BTA6 was found to be consistently significant in both sample sets including and excluding the secondary cluster identified by CMDS (Table 1). The second most significant SNP identified, on BTA25 (rs109960101), did not reach chromosome‐wide significance (−log10 P = 4.91).

Table 1

Association P‐values for the SNP on BTA6 identified in the GWAS for heterozygote advantage on all the animals and after removing the animals in the minor cluster identified through classical multidimensional scaling (CMDS), with Bonferroni‐corrected significance thresholds

Analysis	rs43032684 –log10(P‐value)
All animals	6.21
After removel of CMDS‐identified cluster	6.09
Genome‐wide threshold	7.09
Chromosome‐wide threshold	5.78

Figure 1

Heterozygote disadvantage GWAS Manhattan plot showing SNP associations with the bTB status. The green line represents the chromosome‐wide threshold for BTA6, and the black line represents the genome‐wide threshold.

None of the 80 SNP loci in the two regions identified by Driscoll et al. (2011) and Amos et al. (2013) achieved the lower statistical significance threshold. The most significant SNP locus on chromosomes BTA9 and BTA11 was rs110974556 on BTA9 with chromosomal position 40 944 275 bp, which is ~36 Mb away from the published microsatellite BMS2753, but this did not achieve chromosomal significance.

Genotypic frequencies and deviations from HWE

The genotypic frequencies for rs43032684 on BTA6 are shown in Table 2. Departure from HWE was tested as described in Appendix S1, and maximum likelihood estimates for p and α for each line are shown in Table S2 for the four models tested. Likelihood ratio tests showed that the full model, with both p and α differing between lines, was a significantly better fit than any reduced model. Controls had a greater frequency of allele A (p  = 0.728 cf. p'  = 0.658; P < 0.01) and comparatively fewer heterozygotes, even after accounting for differences in p (α = 0.24 cf. α' = −0.02; P < 0.001) (Fig. S4). The cases showed no significant departure from HWE (α' = 0), whereas the controls showed excess homozygosity compared to HWE (χ 2 = 18.274, 1 d.f.; P < 0.001).

Table 2

Genotypic frequencies for rs43032684 on BTA6 for the cases and the controls

rs43032684	Cases	Controls	Fraction of total
A/A	0.21	0.25	0.46
A/G	0.22	0.13	0.35
G/G	0.05	0.05	0.11
n	0.48	0.44	–
Missing	0.03	0.05	0.08

Regional and genomic additive and dominance variance

The region containing rs43032684, with apparent heterozygote disadvantage, contributed significant dominance variance (P < 0.01) (Table 3). The dominance variance associated with this region remained significant after removing rs43032684 (also P < 0.01), suggesting that the result was not due to SNP‐specific genotyping errors. No dominance variance was detected for either of the regions of microsatellites INRA111 and BMS2753 (Table 3). For the region containing INRA111, no additive genetic variance was detected. For BMS2753, Driscoll et al. (2011) had reported an unspecified genotype effect, and a small additive variance component (P = 0.05) was detected in this region.

Table 3

Likelihood ratio test (LRT) for dominance and additive effects for the regions of interest. The full models contain the additive genomic matrix, the dominance genomic matrix, the local (40‐SNP) additive and the local dominance matrices; the reduced models are without the local dominance matrix. The alternative full model for the SNP with heterozygote disadvantage is for the same region but without the SNP of interest in the calculation of the local matrices

Model	BTA	Marker	Region	LogL	LRT	P‐value	d.f.
Full (1)	6	rs43032684	10 178 841–10 307 829	218.808
Reduced (1)	6	rs43032684	10 178 841–10 307 829	213.603	10.410	<0.01	1
Full (2)	6	rs43032684	10 178 841–10 307 829	218.214
Reduced (2)	6	rs43032684	10 178 841–10 307 829	213.570	9.288		1
Full	9	BMS2753	76 744 316–76 878 698	214.024
Reduced	9	BMS2753	76 744 316–76 878 698	214.024
Full	11	INRA111	40 196 811–40 503 568	211.732
Reduced	11	INRA111	40 196 811–40 503 568	211.732

LogL, log‐likelihood of the corresponding model; LRT = 2[LogL(full Model) – LogL(reduced Model)].

Region exploration

The SNP identified through the heterozygote disadvantage GWAS (chromosomal position on BTA6: 10 245 091 bp) was found to be adjacent to several candidate genes including the translocation associated membrane protein 1‐like 1 gene (TRAM1L1, chromosomal position on BTA6: 9 455 206–9 456 342 bp). A closer look at the surrounding region revealed that the SNP resides within two partly overlapping copy number variation (CNV) regions which were previously associated with resistance to gastrointestinal nematodes in cattle (Ensembl, UMD 3.1; BTA6: 9 524 225–10 295 238 bp and 9 914 188–10 346 486 bp) (Hou et al. 2011, 2012).

Furthermore, the SNP showing heterozygote disadvantage has two alleles (G/A) and resides within the peroxiredoxin‐6 pseudogene (LOC784039, BTA6: 10 223 720–10 246 027 bp, transcript length: 840 bps). The pseudogene contains two exons, and the SNP resides downstream of exon 1. Pseudogenes are homologues of functional parental genes, from which they have been derived through reverse mRNA transcription (processed pseudogenes) or through gene duplication (non‐processed pseudogenes) (Gerstein & Zheng 2006). In Ensembl, blasting the sequence of peroxiredoxin‐6‐like pseudogene against that of the peroxiredoxin‐6 gene (BTA16) confirmed that peroxiredoxin‐6 is the parental gene of peroxiredoxin‐6‐like pseudogene and that the latter is likely to be a processed pseudogene, as the parental gene resides on a different chromosome (Vanin 1985; Hirotsune et al. 2003; Kandouz et al. 2004; Bischof et al. 2006).

Discussion

In this study we examined the hypothesis of a heterozygote (dis)advantage and a role for dominance in genetic resistance to bTB, which was proposed by Amos et al. (2013), and used a case–control population to test published candidate regions and other regions using high‐density SNP data. No evidence was found, neither a SNP effect nor dominance variance, to support the region previously published by Amos et al. (2013), which had been interpreted as showing a heterozygote advantage for resistance. In contrast, the findings of this study suggest that a region on BTA6 (rs43032684) may show dominance variance in bTB resistance, the SNP analysis suggesting a heterozygote disadvantage with heterozygotes more likely to be diseased. However, this association did not achieve genome‐wide significance and so emphasises a need for further validation. The additive effects identified through standard GWAS were consistent with findings previously reported by Bermingham et al. (2014) using the same data, although the two studies used different quality control and statistical procedures.

The putative heterozygote advantage reported by Amos et al. (2013) concerned a candidate microsatellite (INRA111), previously published by Driscoll et al. (2011) in a small mixed breed population. Two approaches to validation were adopted: one analogous to classical GWAS testing for associations with individual SNP loci within the candidate region and the second using regional variance components, which more closely models the haplotype structure of the regions and hence may be more appropriate for validating microsatellites, which were not included in the data. The lack of support for these candidates was despite the targeted testing procedures and strengthens the likelihood that the previous results were false positives. It is possible that the associations described by Amos et al. (2013) are not common to all breeds; however, the Holstein–Friesian gene pool studied here comprised a substantial proportion of the population investigated by Amos et al. (2013).

The region of the rs43032684 SNP on BTA6 was shown to harbour dominance variance in this study, but in contrast to the hypothesis of Driscoll et al. (2011) and Amos et al. (2013), the SNP identified appeared to show heterozygote disadvantage where by heterozygotes appeared to be over‐represented in cases compared to controls. The initial analysis, regressing the phenotype on the binary heterozygosity indicator, identified regions in which heterozygosity differs between cases and controls, but differences in heterozygosity can arise from changes in allele frequency (i.e. an additive effect). However, the analysis of genotype frequencies showed that the differences observed were due primarily to differences in the deviation from HWE between the cases and controls, which in effect, is consistent with some degree of heterozygote (dis)advantage. Nevertheless, the value of α was large and positive in the controls, and it would be expected that the controls would exhibit a degree of departure from HWE that is closer to the whole population than to the cases. This is because the prevalence of bTB is low, so the selection intensity in sampling the controls is smaller than that in sampling the cases. In the absence of selection, the value of α for random mating of female and male parents is expected to be negative but small (Robertson 1961) of the order of the rate of inbreeding, as in the cases, although the mating sub‐populations (such as farms) would tend to shift α towards positive values but to the extent observed for rs43032684. A random sample of 1500 SNP loci from chromosomes showing no evidence of associations with bTB resistance had a mean α = 0.0057 (SE = 0.0014) in the cases and α = 0.005 (SE = 0.0014) in controls and a mean difference between cases and controls of 0.0007 (SE = 0.0018). Therefore, the evidence presented for this region requires interpretation with caution and undoubtedly requires further validation.

The rs43032684 SNP resides within a pseudogene, but pseudogenes may retain some tissue‐specific and/or state‐specific functionality (Hirotsune et al. 2003; Kandouz et al. 2004). The parental gene product is peroxiredoxin‐6, an enzyme localised in lysosomes and active in lipid degradation pathways in various types of cells including alveolar macrophages (Sorokina et al. 2009; Chatterjee et al. 2011). Mycobacteria over‐ride the immune response by releasing lipid mediator molecules within the macrophages, which accumulate in the membrane of the phagosome, prohibiting its fusion with the lysosome, which would result in destroying the mycobacteria (Anes et al. 2003; Koul et al. 2004; Raman et al. 2010), therefore impairing clearance of infection (Gammack et al. 2004; Raman et al. 2010). Despite this apparent functional relevance, expression in monocytes could not be confirmed using PCR (K. Jensen, personal communication, January 18, 2013).

Conclusions

We present here a genome‐wide association method for capturing non‐additive genetic variation, which allows for the identification of associations that are not detected via standard GWAS. In our study, we found no evidence for a region previously reported to have heterozygote advantage that would hinder culling of bTB infected cattle in the current eradication regime. In contrast, our methods identified a novel candidate region on BTA6 associated with bTB resistance, where locus heterozygosity was linked to increased susceptibility to bTB in cattle, i.e. a heterozygote disadvantage. The SNP was located in a pseudogene within a CNV region associated with nematode resistance in cattle; however, despite the apparent functional relevance, functionality was unable to be confirmed. Further studies are required to confirm these findings and validate the QTL on large scale data.

Figure S1 Genotype calls scoring graph for sequencing quality control for rs43032684, where the samples are displayed in three distinct shaded areas based on their genotype calls.

Click here for additional data file.

Appendix S1 Departures from Hardy‐Weinberg equilibrium.

Table S1 Estimates for rs43032684 on BTA6 of the frequency of allele A (pA), and the degree of departure from Hardy‐Weinberg equilibrium (α) for cases and controls with different constraints.

Figure S2 Contour plots for the cases and controls showing departures from HWE.

Click here for additional data file.

Appendix S2 Standard GWAS.

Table S2 Chromosome‐wide significant SNPs identified from standard GWAS and corresponding P‐values.

Figure S3 Manhattan plot from standard GWAS showing significance of SNP associations based on their P‐values.

Figure S4 Q‐Q plot showing observed compared to expected x2 values under the null hypothesis of no association.

Click here for additional data file.