Genome-wide association study of lung lesions and pleurisy in New Zealand lambs.

Pneumonia is an important issue for sheep production, leading to reduced growth rate and a predisposition to pleurisy. The objective of this study was to identify loci associated with pneumonic lesions and pleurisy in New Zealand progeny test lambs. The lungs from 3,572 progeny-test lambs were scored for presence and severity of pneumonic lesions and pleurisy at slaughter. Animals were genotyped using the Illumina Ovine Infinium HD SNP BeadChip (606,006 markers). The heritability of lung lesion score and pleurisy were calculated using the genomic relationship matrix, and genome-wide association analyses were conducted using EMMAX and haplotype trend regression. At slaughter, 35% of lambs had pneumonic lesions, with 9% showing lesions on more than half of any individual lobe. The number of lambs recorded as having pleurisy by the processing plants was 9%. Heritability estimates for pneumonic lesions and pleurisy scores adjusted for heteroscedasticity (CPSa and PLEURa) were 0.16 (± 0.03) and 0.05 (± 0.02), respectively. Five single-nucleotide polymorphisms (SNPs) were significantly associated with pneumonic lesions at the genome-wide level, and additional 37 SNPs were suggestively significant. Four SNPs were significantly associated with pleurisy, with an additional 11 SNPs reaching the suggestive level of significance. There were no regions that overlapped between the 2 traits. Multiple SNPs were in regions that contained genes involved in either the DNA damage response or the innate immune response, including several that had previously been reported to have associations with respiratory disease. Both EMMAX and HTR analyses of pleurisy data showed a significant peak on chromosome 2, located downstream from the transcription factor SP3. SP3 activates or suppresses the expression of numerous genes, including several genes with known functions in the immune system. This study identified several SNPs associated with genes involved in both the innate immune response and the response to DNA damage that are associated with pneumonic lesions and pleurisy in lambs at slaughter. Additionally, the identification in sheep of several SNPs within genes that have previously been associated with the respiratory system in cattle, pigs, rats, and mice indicates that there may be common pathways that underlie the response to invasion by respiratory pathogens in multiple species.

INTRODUCTION

Chronic nonprogressive pneumonia is the most common form of ovine pneumonia in New Zealand, and is an important issue for sheep production, leading to reduced growth rate (Kirton et al., 1976; Alley, 1987; Goodwin et al., 2004; McRae et al., 2016) and a predisposition to pleurisy (Alley, 2002). There is well-documented evidence for between-animal variation in the ability of livestock to resist multiple diseases of economic importance, including respiratory disease (Bishop and Morris, 2007; Davies et al., 2009). Previous work has established that the heritability of pneumonic lesions at slaughter in New Zealand mixed breed progeny tested lambs is 0.07 ± 0.02 (McRae et al., 2016), which is comparable to estimates of respiratory disease in cattle (Snowder, 2009). With growing pressure to reduce the use of antibiotics and drugs in agriculture, these heritable differences mean that improvement of animal health through genetic selection for enhanced resistance can be used as a complementary approach to current methods for disease control (Goddard, 2012). More fundamentally, genomics, through tools such as genome-wide association studies (GWAS), can also be used to further increase our understanding of the genetic mechanisms underlying the host response to disease, and compare these mechanisms between breeds or species. Discovering regions of the genome associated with resistance or susceptibility may also lead to the development of new diagnostic tools and alternative treatments. The aim of this study was to utilize genotype data to identify regions of the ovine genome associated with pneumonic lesions and pleurisy in New Zealand lambs.

MATERIALS AND METHODS

Animals were managed in accordance with the provisions of the New Zealand Animal Welfare Act 1999, and the New Zealand Codes of Welfare developed under sections 68–79 of the Act.

Animals

The lungs from a total of 3,572 ewe and ram lambs from 4 flocks were scored for the presence and severity of pneumonic lesions. Lambs were from 3 South Island (Flocks A, B, and C) and one North Island (Flock D) progeny test flocks. All flocks were fixed-date slaughters, which took place when lambs were between 4 and 8 mo of age. Dams were composites of the main dual-purpose sheep breeds used in New Zealand, including Romney, Coopworth, Perendale, and Texel. Sires were a mixture of dual-purpose and terminal sire composites.

Phenotypic Measurements

The methodology for scoring pneumonic lesions has been previously described (Baird et al., 2012; McRae et al., 2016). Briefly, lungs were scored at chain speed postslaughter at the processing plant. The “consolidated pneumonia score” (CPS) system has a range from 0 to 2, where 0 = no lesions present; 1 = any individual lobe with up to 50% of the lobe affected and 2 = any individual lobe with greater than 50% of the lobe affected. Pneumonic lesions were defined as compacted, dark purple-red areas of the lung that were firm to touch. Information on lamb carcasses that were identified as having pleurisy and detained for trimming was obtained from the processing plants.

Data cleaning consisted of removal of records with 1) missing values, and 2) contemporary groups (CG) containing less than 5 observations. CG was defined as flock, birth year, sex, weaning mob, and slaughter date; animals needed to have all of these in common to be considered in the same CG. Weaning mob was obtained from Sheep Improvement Limited (SIL), the New Zealand sheep genetic evaluation database. To adjust for heteroscedasticity, CPS (initially scored as 0, 1, or 2) was scaled using the formula CPSa = CPS/SQRT[CPSm*(2-CPSm)], where m is the mean incidence rate within the CG where phenotypic score is being adjusted. Pleurisy (initially coded as 0 or 1) values were also transformed using the formula PLEURa = PLEUR/SQRT[PLEURm*(1-PLEURm)].

Genotypes and Quality Control

Genomic DNA was extracted from ear tissue samples collected from lambs at tailing, using a high-throughput DNA extraction method (Clarke et al., 2014). Animals were genotyped with the Illumina Ovine Infinium HD SNP BeadChip (606,006 markers) according to the manufacturer’s protocol. Genome coordinates of each single-nucleotide polymorphism (SNP) were based on the OARv3.1 ovine genome assembly (Jiang et al., 2014). Quality control checks excluded markers that appeared nonautosomal (including pseudoautosomal), had a call rate below 90%, and/or had a minor allele frequency (MAF) ≤0.01. Individuals were excluded from the analysis if there was more than 5% genotyping failure. After quality control measures, 3,546 phenotyped animals were available, with 537,117 of the initial 606,006 SNPs utilized for analysis.

Heritability

Variance components were estimated using restricted maximum likelihood (REML) procedures fitting an animal model in ASReml (Gilmour et al., 2015), with the genomic relationship matrix (GRM) estimated in GenABEL (Aulchenko et al., 2007) using HD genotypes. Heritabilities were obtained by running a univariate analysis on the respective traits. Data analysis models for both CPSa and PLEURa, included CG as a fixed effect (McRae et al., 2016).

Genome-wide Association Analyses

Pneumonic lesion and pleurisy data were analyzed using values adjusted for heteroscedasticity (CPSa and PLEURa). Pneumonic lesion data was also analyzed by only including animals with no lesions and those with severe lesions [CPSa (0&2)]. Genome-wide association analyses were performed using SNP & Variation Suite v8.4.0 (Golden Helix, Inc., Bozeman, MT, www.goldenhelix.com) using 2 of the following approaches: 1) Efficient Mixed-Model Association eXpedited (EMMAX) using identity-by-state (IBS), and 2) haplotype trend regression (HTR) with a 3-SNP sliding window. Analyses were performed on adjusted values, with CG fitted as a covariate. Genome browse software was used to visualize results with an added track of Ovis aries genes from Ensembl 84. After Bonferonni correction, thresholds were 9.31 × 10−8 and 1.86 × 10−6 for genome-wide significance (P < 0.05) and suggestive significance (P < 0.1), respectively.

RESULTS

Incidence of Pneumonia and Pleurisy

In total, 3,572 lungs were scored for pneumonic lesions from lambs born between 2013 and 2015 (Table 1). Of these, 1,234 (35%) had lesions, with 329 (9%) showing lesions on more than 50% of any individual lobe (CPS of 2). The number of lambs recorded as having pleurisy by the processing plants was 310 (9%). The incidence of pneumonia was significantly higher in 2014-born lambs than those born in 2013 or 2015 (P < 0.001). Of the 310 animals recorded as having pleurisy, 118 (38%) had a CPS of 0, 71 (23%) had a score of 1, 60 (19%) had a score of 1, and 61 (20%) were unable to be scored due to the lungs being retained in the carcass.

Table 1.

Incidence of pneumonic lesions and pleurisy by flock and year of birth

Flock	Year born	Lungs scored	CPSa > 0	CPSa = 2	Pleurisy
A	2013	292	52 (18%)	20 (7%)	25 (9%)
	2014	483	194 (40%)	57 (12%)	70 (14%)
	2015	467	98 (21%)	29 (6%)	77 (16%)
B	2014	766	334 (44%)	80 (10%)	59 (8%)
	2015	292	46 (16%)	11 (4%)	6 (2%)
C	2015	56	14 (25%)	3 (5%)	13 (23%)
D	2014	1,216	496 (41%)	129 (11%)	60 (5%)
	Total	3,572	1,234 (35%)	329 (9%)	310 (9%)

aCPS = Consolidated Pneumonia Score, where 0 = no lesions present; 1 = individual lobes with up to 50% of the lobe affected and 2 = individual lobes with greater than 50% of the lobe affected.

The heritability estimated for CPSa was 0.16 (± 0.03) and for PLEURa was 0.05 (± 0.02). This is slightly higher than the previously published estimates of 0.07 ± 0.02 and 0.02 ± 0.01, respectively (McRae et al., 2016). This is likely to be due to the use of a GRM rather than recorded pedigree information in estimating the heritability in the current analysis; the 1,216 lambs from Flock D were included in both studies, however only sire information is recorded for these animals, therefore using a GRM rather than pedigree is a more accurate estimation of relatedness. The genetic correlation between the 2 traits was 0.58 (± 0.16), and the phenotypic correlation was (0.15 ± 0.02), which was in line with previous estimates.

When adjusted pneumonic lesion (CPSa) information from all animals was included, there were no SNPs that passed the threshold for suggestive significance in the EMMAX analysis (Fig. 1A). In the HTR analysis, however, 4 regions, on chromosomes 3, 6, 8, and 13, were significant (Fig. 1B), with a further 31 SNPs passing the level for suggestive significance (Table 2). When only the extreme animals were included (i.e., animals with no lesions compared to those with severe lesions; CPS of 0 vs. CPS of 2), several SNPs in each analysis were of suggestive significance, although none reached the genome-wide level of significance (Fig. 1C and D). The top 3 SNPs in the EMMAX analysis were all intronic variants, in the LSAMP, PPIL6, and KCNMA1 genes (Table 2). The top SNPs in the HTR analysis included the same SNP in EYA4 that reached the suggestive significance level in the CPSa HTR analysis, along with 2 missense variants in exon 2 of ATAD5. Both EMMAX and HTR analyses of pleurisy data showed a significant peak on chromosome 2 (Fig. 2). Additionally, there were multiple suggestively significant intergenic SNPs on chromosomes 8 and 11 (Table 2).

Figure 1.

Manhattan plot of genome-wide association analysis for consolidated pneumonia score (CPS) in New Zealand lambs. Analyses were performed on all animals (A and B), or only including animals with scores of 0 or 2 (C and D). Genome-wide association analyses were conducted using 2 approaches: 1) Efficient Mixed-Model Association eXpedited (EMMAX) using identity-by-state (IBS) (A and C), and 2) haplotype trend regression (HTR) with a 3-SNP sliding window (B and D). Analyses were performed on pneumonic lesion scores after adjustment for heteroscedasticity, with contemporary group fitted as a covariate.

Table 2.

SNPs suggestively and significantly associated with consolidated pneumonia score (CPS) and pleurisy in New Zealand lambs

Testa	Analysisb	P-valuec	Chr	Position	RSID	Gened	Gene named	Variant consequence (impact)d
CPSa	HTR	1.72E-06	1	198656040	rs421794454	ENSOARG00000020512	RFC4	Intron variant
	HTR	1.02E-06	1	198674451	rs416081302	Multiple small nucleolar RNAs		Upstream gene variant
	HTR	8.00E-07	2	158991503	rs407273673
	HTR	1.15E-06	2	158998132	rs417102378
	HTR	8.05E-07	3	197683811	rs426850802
	HTR	8.74E-07	3	197705247	rs405096150
	HTR	1.72E-07	3	197720146	rs430716198
	HTR	1.34E-07	3	197720936	rs412869687
	HTR	3.80E-07	3	197824787	rs407726225
	HTR	4.75E-08*	3	197825391	rs424070250
	HTR	1.06E-06	4	100364303	rs403394816
	HTR	6.45E-08*	6	77843695	rs399606595
	HTR	1.63E-09*	8	7733798	rs429357466
	HTR	7.65E-09*	8	7743164	rs400905064
	HTR	4.29E-07	8	12061122	rs425423371
	HTR	1.01E-06	8	12061431	rs402511423
	HTR	1.85E-06	8	12149272	rs423436094
	HTR	6.85E-07	8	46753989	rs422310670
	HTR	4.64E-07	8	58610312	rs418966278	ENSOARG00000014564	EYA4	Intron variant
	HTR	1.10E-06	8	65817474	rs419752214
	HTR	2.43E-07	8	65827376	rs399425501
	HTR	3.76E-07	8	65834390	rs429368446
	HTR	6.06E-07	9	56587375	rs424700173
	HTR	8.75E-07	12	12552747	rs408273790
	HTR	1.01E-06	13	8840940	rs400804234
	HTR	3.36E-08*	13	8848881	rs417728121
	HTR	2.36E-07	13	8850719	rs423025524
	HTR	4.67E-07	13	8855334	rs416260513
	HTR	1.49E-06	13	8860196	rs420406541
	HTR	3.06E-07	13	8864106	rs428620400
	HTR	4.21E-07	13	8871693	rs430812458
	HTR	9.09E-07	13	8872416	rs404379883
	HTR	8.25E-07	18	51969325	rs419274927	ENSOARG00000026456	Novel lincRNA	Noncoding transcript variant
	HTR	1.26E-06	20	18894279	rs401389671
	HTR	1.37E-06	20	23420251	rs417356235
	HTR	1.10E-06	20	23423229	rs399485900
CPSa (0&2)	EMMAX	7.91E-07	1	178857472	rs410655004	ENSOARG00000019641	LSAMP	Intron variant
	EMMAX	1.42E-06	8	27846793	rs430024463	ENSOARG00000010461	PPIL6	Intron variant
	EMMAX	1.29E-06	25	32872254	rs422854508	ENSOARG00000009163	KCNMA1	Intron variant
	HTR	9.49E-07	4	110104884	rs425466808
	HTR	1.39E-06	8	7743164	rs400905064
	HTR	2.23E-07	8	46753989	rs422310670
	HTR	1.43E-06	8	58610312	rs418966278	ENSOARG00000014564	EYA4	Intron variant
	HTR	8.91E-07	11	17720321	rs419581914	ENSOARG00000012322	ATAD5	Missense variant (moderate)
	HTR	1.79E-06	11	17720415	rs400520703	ENSOARG00000012322	ATAD5	Missense variant (moderate)
PLEURa	EMMAX	3.10E-09*	2	134984962	rs398681238
	EMMAX	2.32E-08*	2	134985148	rs424471052
	EMMAX	1.57E-06	2	135163372	rs415671617
	EMMAX	1.48E-06	2	242723012	rs421193149
	EMMAX	1.58E-06	11	15261261	rs420254502
	EMMAX	3.29E-07	11	15262540	rs409974296
	EMMAX	1.80E-07	11	15265356	rs417033802
	HTR	1.38E-07	2	134976058	rs404285802	ENSOARG00000000469	SP3	Downstream gene variant
	HTR	7.66E-08*	2	134979525	rs428634189
	HTR	1.53E-08*	2	134984962	rs398681238
	HTR	6.59E-07	2	134985148	rs424471052
	HTR	1.88E-07	2	134998369	rs414115266
	HTR	1.82E-06	2	135006264	rs412779979
	HTR	1.72E-06	8	13863996	rs414046873
	HTR	9.86E-07	8	88651287	rs412134993
	HTR	4.03E-07	8	88659717	rs398705894

aAnalyses were performed on consolidated pneumonia and pleurisy scores after adjustment for heteroscedasticity (CPSa and PLEURa, respectively). For CPSa data, analyses were performed using all animals, or only including animals with scores of 0 or 2 [CPSa (0&2)].

bGenome-wide association analyses were conducted using 2 approaches: 1) Efficient Mixed-Model Association eXpedited (EMMAX) using identity-by-state (IBS), and 2) haplotype trend regression (HTR) with a 3-SNP sliding window. Contemporary group (sex, birth year, flock, weaning mob, and kill date) was fitted as a covariate in all analyses.

cAfter Bonferonni correction, thresholds were 9.31 × 10−8 and 1.86 × 10−6 for genome-wide significance (P < 0.05*) and suggestive significance (P < 0.1), respectively.

dGene names and variant consequences are based on Ensembl Release 84.

Figure 2.

Manhattan plot of genome-wide association analysis for pleurisy in New Zealand lambs. Genome-wide association analyses were conducted using 2 approaches: 1) Efficient Mixed-Model Association eXpedited (EMMAX) using identity-by-state (IBS) (A), and 2) haplotype trend regression (HTR) with a 3-SNP sliding window (B). Analyses were performed on pleurisy scores after adjustment for heteroscedasticity, with contemporary group fitted as a covariate.

DISCUSSION

In sheep, as with other ruminants, respiratory disease such as pneumonia is etiologically complex, resulting from a complex interaction between multiple infectious agents and the host, which is often compromised by physical and physiological stress. GWAS help provide an understanding of the genes and pathways involved in the response to disease. GWAS in both dairy (Neibergs et al., 2014) and beef (Keele et al., 2015) cattle have identified multiple loci associated with bovine respiratory disease complex (BRDC). Neibergs et al. (2014) discovered candidate loci involved in viral susceptibility, viral entry into cells, and modulation of inflammation in a case–control analysis of preweaned Holstein calves. A GWAS of lung lesions in beef cattle identified SNPs near candidate genes involved in functions such as tissue repair and regeneration, cell proliferation, apoptosis, and immunity (Keele et al., 2015).

The majority of SNPs associated with pneumonic lesions in this study were in intergenic regions of the sheep genome. Intergenic variants within RFC4, EYA4, and a novel lincRNA were suggestively associated with pneumonic lesions when including all the data and variants within LSAMP, PPIL6, and KCNMA1 reached suggestive significance when only including the extreme animals. Additionally, 2 missense variants in exon 2 of ATAD5 also reached the suggestively significant level in the analysis of the extreme animals.

EYA4, ATAD5, and RFC4 all have roles in the response to DNA damage. Eyes Absent (EYA) proteins are implicated in a diverse range of processes, including DNA damage repair and innate immunity (Tadjuidje and Hegde, 2013). EYA4 has been shown to enhance the innate immune response to viruses through stimulating the interferon regulatory factor 3 (IRF3)-mediated transcription of IFN-β and CXCL10 in response to undigested DNA (Okabe et al., 2009). EYA4 has been associated with familial lung cancer risk (Wilson et al., 2014), and a SNP located within 15 kb of EYA4 has been significantly associated with lung lesions in commercial beef cattle (Keele et al., 2015).

The replication factor C (RFC) complex, composed of subunits RFC1-5, also plays an essential role in DNA replication and repair in eukaryotes (Kim and MacNeill, 2003). Additionally, several RFC-like complexes (RLC), made up of RFC2-5 and an alternative subunit that replaces RFC1, have been reported, including ATAD5-RLC (Ben-Aroya et al., 2003). Atad5+/− mice show high levels of genomic instability (Bell et al., 2011), and delayed DNA replication and cell division, leading to an altered adaptive immune response though reduced immunoglobulin class switching (Zanotti et al., 2015). The identification of 2 suggestively significant SNPs within genes that form ATAD5-RLC highlights the potential importance of this complex in the host response to respiratory challenge.

Although not their primary role, both PPIL6 and KCNMA1 have previously been associated with the respiratory system. A QTL containing the cyclophilin-like PPIL6 was associated with the variability of immune response in a crossbred swine population postinfluenza vaccination (Zanella et al., 2015). The potassium channel gene KCNMA1 was expressed at significantly higher levels in the lungs of asthmatic rats compared to those of control rats (Yin et al., 2008), and is differentially methylated during normal development in the mouse and human lung (Cuna et al., 2015).

The significant peak on chromosome 2 associated with pleurisy was detected using 2 independent methods. This peak is located downstream from the transcription factor SP3, which is involved in the activation or suppression the expression of numerous genes, including the interferon regulatory factor IRF3 and IL-10, an anti-inflammatory cytokine (Tone et al., 2000). Of interest is that Sp3 knockout mice die at birth of respiratory failure, although only minor structural abnormalities are observed in the lungs (Bouwman et al., 2000). As mentioned above, IRF3 is involved in the innate response to viral infection (Xu et al., 2012), and several bovine viral pathogens including bovine herpesvirus 1 (BHV-1) and bovine diarrhoea virus (BVDV) target IRF3 activity, halting the interferon response (Srikumaran et al., 2007). As with other farmed ruminants, in sheep, pneumonia is etiologically complex. While Mannheimia haemolytica is considered to be the predominant agent responsible for lung damage, multiple viruses (Davies et al., 1977; Davies et al., 1982; Davies and Jones, 1985) have also been shown to play a role through compromising the respiratory system, allowing secondary invasion by bacteria (Brogden et al., 1998). An enhanced immune response to viruses could therefore result in a reduced chance of developing lung damage.

As previously discussed, pneumonia can arise through a combination of a variety of environmental and pathogenic factors. Despite the complex nature of this disease, previous research in both sheep and cattle has shown that there is an underlying genetic component in the variation observed between animals in their susceptibility to pneumonia (Snowder, 2009; McRae et al., 2016). This indicates that it is possible to select for animals with the ability to withstand and/or recover from infection that can be the result of multiple causative factors. This study identified several SNPs associated with genes involved in both the innate immune response and the response to DNA damage that are associated with pneumonic lesions and pleurisy in lambs at slaughter. Additionally, the identification in sheep of several SNPs within genes that had previously been reported to be involved in the respiratory system in cattle, pigs, rats, and mice indicates that there may be common genetic pathways underlying the response to respiratory disease in multiple mammalian species.