Genetic susceptibility to feline infectious peritonitis in Birman cats.

Genetic factors are presumed to influence the incidence of feline infectious peritonitis (FIP), especially among pedigreed cats. However, proof for the existence of such factors has been limited and mainly anecdotal. Therefore, we sought evidence for genetic susceptibility to FIP using feline high density single nucleotide polymorphism (SNP) arrays in a genome-wide association study (GWAS). Birman cats were chosen for GWAS because they are highly inbred and suffer a high incidence of FIP. DNA from 38 Birman cats that died of FIP and 161 healthy cats from breeders in Denmark and USA were selected for genotyping using 63K SNPs distributed across the feline genome. Danish and American Birman cats were closely related and the populations were therefore combined and analyzed in two manners: (1) all cases (FIP) vs. all controls (healthy) regardless of age, and (2) cases 1½ years of age and younger (most susceptible) vs. controls 2 years of age and older (most resistant). GWAS of the second cohort was most productive in identifying significant genome-wide associations between case and control cats. Four peaks of association with FIP susceptibility were identified, with two being identified on both analyses. Five candidate genes ELMO1, RRAGA, TNFSF10, ERAP1 and ERAP2, all relevant to what is known about FIP virus pathogenesis, were identified but no single association was fully concordant with the disease phenotype. Difficulties in doing GWAS in cats and interrogating complex genetic traits were discussed.

Introduction

Feline infectious peritonitis (FIP) is caused by internally arising 3c and S gene mutants of a ubiquitous and largely non-pathogenic feline enteric coronavirus (FECV) (Chang et al., 2012, Pedersen et al., 2012). The resultant FIP viruses (FIPVs) cause disease mainly among young cats from pedigreed catteries, shelters, and other dense multi-cat environments. FIP prevalence is around 1–5% in these populations, but mortality among affected cats is virtually 100% once clinical signs appear (Pedersen, 2009). This suggests that either the FECV-to-FIPV mutation rate is very low or that most infections are abortive. The first scenario is unlikely based on two different studies. Chang et al. (2012) calculated that thousands of FECV progeny carrying the S protein mutation occur during any single infection. Another study indicated that the FECV-to-FIPV mutation and FIP occurs in at least 20% of older cats following FECV exposure, but only in individuals that were immunocompromised (Poland et al., 1996). Therefore, additional extrinsic (environment) or intrinsic (host) factors must play an important role in actual disease occurrence.

Extrinsic factors that increase the level of exposure to FECV (Foley et al., 1997) or that depress the immune system at the time of initial FECV infection (reviewed Pedersen, 2009, Poland et al., 1996) affect FIP prevalence. Age at the time of initial FECV exposure, which usually occurs around 9 weeks of age (Pedersen et al., 2008), may also be a critical factor, because the adaptive immune system of kittens is still maturing during this period (Pedersen, 1987). Kittens from dense multi-cat environments such as catteries and shelters are often under heavy stress due to the effects of overpopulation, diet changes, early weaning and exposure to a myriad of other common kitten-hood infections (Pedersen, 2009).

Susceptibility to FIPV is also affected by factors that are intrinsic to the host. Confirmed cases of FIP did not exist before the 1950s (reviewed Pedersen, 2009). As such, cats have not had a long period of time to evolve significant immunity to the virus. Pedigreed cats are more likely to develop FIP than random-bred cats and certain breeds are more likely to succumb to FIP than others (Pesteanu-Somogyi et al., 2006, Rohrbach et al., 2001, Worthing et al., 2012). One study of Persian catteries and pedigrees indicated that susceptibility to FIP was at least 50% heritable (Foley and Pedersen, 1996). Identifying a genetic association with FIP susceptibility may be possible using a newly available moderate density feline SNP array with a highly inbred cat population. Studies of disease associations with inbred dog breeds have identified a large number of single gene mutations, as well as a complex genetic trait, using <100 cases and controls and relatively few (20–172K) SNP markers (reviewed in Ostrander, 2012). There has been much less experience with GWAS in cats than dogs, although many phenotypic and disease traits of cats have been identified by candidate gene approaches (Bighignoli et al., 2007, Filler et al., 2012, Gandolfi et al., 2010, Grahn et al., 2012) and a recent GWAS study discovered the causative mutation for hypokalemia in Burmese using just 35 cases and 25 controls (Gandolfi et al., 2012).

This was the first attempt to elucidate genes associated with a complex disease susceptibility trait in cats using GWAS, which is a far more difficult task than identifying causative genes for simple Mendelian traits. However, given sufficient numbers of cases and controls, it may be possible to identify regions of the genome containing relevant candidate genes. The Birman breed is ideal for such a study, given its high incidence of FIP (Pesteanu-Somogyi et al., 2006), relatively small population size and high linkage disequilibrium (Alhaddad et al., 2013).

Materials and methods

Cats used in study

Buccal swabs of 199 Birman cats were submitted by breeders, veterinarians and owners from Denmark and USA (Table 1 ), along with three generation pedigrees. The diagnosis of FIP was confirmed by attending veterinarians based on clinical signs, blood tests, analysis of ascites or pleural fluids, and in some cases by gross or complete necropsies. The male to female ratio of cats in the study was 1.3:1 and ~70% of cats dying of FIP were 1½ years of age or younger (Fig. 1 ). No cats greater than 5½ years of age in this study developed FIP. Sixty four percent of cats died with the non-effusive form of FIP and 36% with effusive disease. The gender and age distribution was typical for FIP, while the predominance of non-effusive over effusive FIP was atypical (Pedersen, 2009).

Table 1

Numbers of cases and controls from Denmark and the United States used in GWAS.

	Samples used in GWAS	Samples after GWAS filtration
Analysis-1
Cases (Denmark)	23	19
Cases (US)	15	15
Total cases	38	34
Controls (Denmark)	121	114
Controls (US)	40	37
Total controls	161	151
Analysis-2
Cases (Denmark)	8	8
Cases (US)	12	12
Total cases	20	20
Controls (Denmark)	63	63
Controls (US)	25	25
Total controls	88	88

Fig. 1

Ages of Birman cats used in the study. Total number of cats with known ages = 176 (cases = 30; controls = 146).

DNA extraction and genotyping

DNA from buccal swabs was extracted following the manufacturer's recommendations (Puregene DNA buccal cell kit; Gentra Systems, Inc., Minneapolis, MN). DNA samples were genotyped with the Illumina Infinium HD feline 63K single nucleotide polymorphism (SNP) BeadChip® Array (Illumina Inc., San Diego, CA) by GeneSeek Inc., Lincoln, NE. A large number of SNPs (n  = 6890) were not assigned to a chromosome or position at the time the array was developed and were therefore packaged on an artificial chromosome i.e., “chr.Un.”. The correct chromosome and position assignments for the 63K SNPs for this study was based upon the ICGSC Felis_catus 6.2/felCat5 (September 2011) assembly on the UCSC genome browser.

Association analysis

The 199 Birman DNA samples were genotyped and genome-wide association (GWAS) studies were performed with PLINK (Purcell et al., 2007). All markers were subject to stringent quality-control criteria to assure high-quality data. Fourteen samples were excluded due to low genotyping (MIND > 0.1, which excluded individuals with more than 10% missing genotype). The total genotyping rate in remaining individuals was 0.99. Some SNPs (n  = 811) failed missingness test (GENO > 0.1, which excluded SNPs with 10% missing genotyping rate). Some SNPs (n  = 25,484) failed the frequency test (MAF < 0.05, exclusion of SNPs with less than 5% minor allele frequency). After frequency and genotyping pruning, 36,862 SNPs remained for analysis. Population stratification in case and control populations from each region were assessed by multidimensional scaling (MDS) and relatedness was confirmed by Pi-Hat test using the PLINK whole genome association analysis toolset (Purcell et al., 2007). After filtering, GWAS data from 185 cats (cases = 34, controls = 151) remained for further analysis (Table 1). Allele frequencies were compared between the cases and controls, and the strength of associations between specific SNPs and disease was shown in terms of P(raw).

Results

The first step in this study was to select case and control populations for GWAS based on available information. Multidimensional scaling of case and control cats from both the USA and Denmark indicated that the two populations were closely related but geographically distinguishable, but with no stratification between cases and controls within the same countries (Fig. 2 ). Given these facts, and assuming that FIP susceptibility traits have been under neutral selection since 1925 when the breed was formally recognized, a decision was made to first compare all case and control cats regardless of country of origin (Table 1). A second analysis involved stratifying cases and controls by limiting cases to cats 1½ years of age or younger (the most susceptible) and the controls to cats 2 years of age or older (the most resistant) (Fig. 1 and Table 1). Genomic inflation factor (GIF) values for both analysis groups were 1.35 and 1.43, respectively. The highest ranking SNPs in both analyses are shown in Table 2 . All potential candidate genes involved with immune responses within 500 kb of each SNP were identified, and if no candidate genes were identified the search was extended to a distance of 1 Mb.

Fig. 2

MDS plot showing the distribution of four groups of 185 Birman cats (DK cases = 19, DK controls = 114, US cases = 15, US controls = 37) based on genome-wide IBS distance and geographic location. X-axis represents first MDS principle component and Y-axis represents second MDS principle component.

Table 2

Highest ranking SNPs from analysis-1 (34 cases vs. 151 controls) and analysis-2 (20 cases vs. 88 controls) based on unadjusted P(raw) and corrected P(Bonf) probabilities, and candidate genes known to be strongly involved in immune responses.

Chr	SNP	SNP position	P(raw)	P(Bonf)	Candidate gene
Analysis-1
A2	A2.191286425	126618108	1.3 × 10−6	0.048	ELMO1
C2	C2.107344890	93184134	7.8 × 10−6	0.288	TNFSF10
D4	D4.37641174	45868722	6.6 × 10−6	0.243	RRAGA
Analysis-2
A1	A1.196617776	154265118	2.5 × 10−7	0.009	ERAP1,2
A1	A1.206840008	164728174	4.1 × 10−7	0.015	ERAP1,2
A1	Un.59861682	155715831	4.3 × 10−7	0.016	ERAP1,2
E2	E2.65509996	54165589	1.2 × 10−6	0.043
A3	A3.140889392	129597254	3.7 × 10−6	0.137
A1	Un3.4893672	160942986	8.2 × 10−6	0.259	ERAP1,2
A2	A2.191286425	126618108	1.4 × 10−5	0.399	ELMO1
C2	C2. 107021201	92907113	1.2 × 10−5	0.347	TNFSF10

When the 34 FIP cases and 151 controls were analyzed (analysis-1), the strongest associations (−log10   P(raw) ≥ 5.0) were found on chr.A2, chr.C2 and chr.D4 (Table 2, Fig. 3a). Given a Bonferroni correction of −log10   P(raw) ≥ 5.8 for genome-wide significance, only the highest ranked SNP on chr.A2 attained significance (p  = 0.048) (Table 2). The Manhattan plot of analysis-2 cats demonstrated 3 peaks with 6 SNPs having −log10   P(raw) ≥ 5.0 on chr.A1, chr.A3, and chr.E2 (Table 2, Fig. 3b). Only the three highest ranking SNPs on chr.A1 and highest ranked SNP on chr.E2 had genome-wide significance after Bonferroni correction (p  = 0.009–0.043) (Table 2). Two SNPs on chr.A2 (log10   P(raw) = 4.86) and chr.C2 (−log10   P(raw) = 4.93) were just below the cutoff at −log10   P(raw) ≥ 5.0 (Table 2). The SNP on chr.A2 was identical to the highest ranking SNP identified in analysis-1, further bolstering its significance. The chr.C2 SNP identified in analysis-2 was separated by only 300 kb from the second ranking SNP on analysis-1. Although neither SNP was significant upon Bonferroni correction, their proximity and high ranking by both analyses was taken as reason to also interrogate this chr.C2 region for candidate genes. GWAS of analysis-1 and -2 cats identified five potential candidate genes on four different chromosomes: engulfment and cell motility protein 1 (ELMO1), Ras-related GTP binding A (RRAGA), tumor necrosis factor superfamily 10 (TNFSF10), and endoplasmic reticulum aminopeptidase 1 and 2 (ERAP1 and ERAP2) (Table 2). No candidate genes were identified on chr.A3 and chr.E2, both of which demonstrated strong associations in analysis-2. Only a single high ranking SNP was associated with ELMO1 (Fig. 4 ) and RRAGA (Fig. 5 ), while several high ranking SNPs were associated with TNFSF10 (Fig. 6 ). Nineteen of 30 highest ranking SNPs in analysis-2, including the top three ranked SNPs, were on chr.A1 (data not shown). Seventeen of the 19 SNPs formed a 10.5 Mb region flanked by the three highest ranking SNPs (Fig. 7 ). This region contained 31 or more genes, two of which (ERAP1 and ERAP2) were strong candidates at positions chr.A1.158,772,633 and chr.A1.158,978,956 bp, respectively. These two genes were near the center of the 10.5 Mb region.

Fig. 3

Genome-wide association analysis of combined and stratified populations: (a) analysis-1 containing 34 cases and 151 controls from both the USA and Denmark; (b) analysis-2 containing 20 cases (1½ years and younger) and 88 controls (2 years or older). Bonferonni correction significant at −log10 P(raw) = 5.8. “Chromosome Un” contains a small number of unmapped SNPs. Peaks > 5.0 − log10 P(raw) were investigated for candidate genes.

Fig. 4

Manhattan plot of chr.A2 based on GWAS of 34 FIP and 151 healthy cats (analysis-1). Peak SNP is A2.191286425 near candidate gene ELMO1.

Fig. 5

Manhattan plot of Chr.D4 based on GWAS of 34 FIP and 151 healthy cats (analysis-1). Peak SNP is D4.37641174 near candidate gene RRAGA.

Fig. 6

(a) Manhattan plot of chr.C2 based on GWAS of 34 FIP and 151 healthy cats (analysis-1). Peak SNP is C2.107344890 near candidate gene TNFSF10 (TRAIL). (b) Manhattan plot of chr.C2 based on GWAS of 20 FIP and 88 healthy cats (analysis-2). Peak SNP is C2.107021201 near candidate gene TNFSF10 (TRAIL).

Fig. 7

Manhattan plot of chr.A1 based on GWAS of 20 FIP and 88 healthy cats (analysis-2). Seventeen of 30 highest ranking SNPs in GWAS were on chr.A1 and within a 10.5 Mb region flanked by the three highest ranking SNPs.

Discussion

This was the first attempt to identify genes associated with FIP susceptibility; a trait presumed to be genetically complex. Its complexity was confirmed by the multiplicity of strong associations between case and control cats observed in GWAS. However, none of these associations was totally concordant with the FIP susceptibility phenotype. Therefore, the associations need to be further confirmed with GWAS of additional case and control cats and sequencing of candidate genes for polymorphisms.

There are findings in studies of this type that either diminish or increase confidence in a particular candidate gene. Candidate genes on chr.A2 (ELMO1) and chr.C2 (TNFSF10) were identified by two different GWAS configurations. However, only a single high ranking SNP was associated with RRAGA and ELMO1, which was somewhat discomforting. In contrast, a number of adjacent high ranking SNPs were associated with TNFSF10, ERAP1 and ERAP2, giving more confidence to these associations. The distance between high ranking SNPs and potential candidate genes was another consideration, and the preference is for high ranking SNPs to lie within a candidate gene's coding region. This is infrequently the case, so the scanning regions are often set at 500 kb or more. However, confirmed candidate genes for cat phenotypic and disease traits are often found 1 Mb or more distant from their identifying SNPs (Gandolfi et al., 2012). Candidate genes TNFSF10 and RRAGA were 69–359 kb away from their respective highest ranking SNPs on chr.C2 and chr.D4, while ELMO1 was 814 kb from the highest ranking SNP on chr.A2. ERAP1 and ERAP2, highly relevant candidate genes, were situated in the center of a 10.5 Mb region flanked by the three highest ranking SNPs and strongly associated with the disease phenotype. Although distant, the SNP structure of this region and the relevance of these two genes to FIP immunopathogenesis made them obvious candidates for further study.

Analysis of GWAS data was also complicated by case/control numbers and population substructure apparent in MDS and from GIF values. Population substructure exists in all breeds of cats because of their short generation time, large numbers of offspring, and intense selection pressures applied by individual breeders. These factors were highlighted by genetic comparisons of Danish and American Birman cats, which evolved from a common gene pool and were segregated for only a few decades by geography. Although highly related as a breed, Danish and American cats were distinguishable by MDS, with the latter being more homogenous than the former. This could only have resulted from subtle phenotypic preferences among breeders in the two countries.

Numbers of cases and controls is another confounding factor in identifying complex genetic traits, even within relatively inbred populations. The removal or addition of relatively few animals in the current study caused strong peaks of association to lose significance or disappear altogether, weaker peaks to gain significance, and new peaks to appear. This effect was evident between GWAS of analysis-1 (all cases vs. all controls, regardless of age or geography) and analysis-2 (young cats with FIP vs. older healthy cats, regardless of geography).

The actual disease and exposure status of case and control is another covariate in a study of this type. There was a high degree of confidence in the FIP diagnoses, because they were made by veterinarians knowledgeable about FIP and who often worked closely with breeders on the problem. The exposure status of control cats was much more difficult to determine, even though FECV is enzootic in virtually all catteries (Pedersen, 2009) and the FECV to FIPV mutation rate is relatively high (Poland et al., 1996). The logic of stratifying the larger combined American and Danish populations into a smaller analysis group according to age was based on a single fact–exposure and susceptibility to FIPV was only certain for those cats that died of FIP. Extending the age for control cats enhanced the chance for FIPV exposure and therefore increased the odds that healthy cats surviving beyond this time were resistant. Stratifying the GWAS population in this manner did produce stronger associations, bolstered confidence in candidate genes such as ELMO1 and TNFSF10 and identified promising new candidate genes ERAP1 and ERAP2.

The potential candidate genes identified in this study have been strongly implicated in immune pathways concerning cellular migration, phagocytosis, programmed cell death (apoptosis), and virus/host cell interactions. ELMO1 (engulfment and cell motility 1) is involved in cytoskeletal rearrangements required for cell motility and phagocytosis of apoptotic cells (Gumienny et al., 2001). ELMO1 has been associated with the inhibition of HIV-1 transcription in monocyte/macrophages in the face of apoptotic events occurring in the T cell compartment (Gekonge et al., 2006). TNFSF10, also known as TNF-related apoptosis-inducing ligand (TRAIL), plays an important role in induction of cell death by apoptosis and is indispensable for a normal functioning immune system (Falschlehner et al., 2009). TRAIL can be expressed by various cells of the immune system, including T cells, natural killer (NK) cells, natural killer T cells, dendritic cells and macrophages. Type I interferon and interferon regulatory factors also regulate TNF-related apoptosis-inducing ligand (TRAIL) in HIV-1-infected macrophages (Herold et al., 2009). Highly virulent strains of avian influenza virus induce high levels of TRAIL in macrophages, which in turn may contribute to the disease pathogenesis (Zhou et al., 2006). Ebola virus infection is known to greatly disrupt immunity by infecting phagocytic cells and releasing cytokines that induce apoptosis in bystander lymphocytes (Hensley et al., 2002). RRAGA encodes a protein involved in the RCC1/Ran-GTPase pathway, which plays a direct role in a TNF-alpha signaling pathway leading to apoptosis. RRAGA may also cause cell death by acting as a cellular target for adenovirus E3-14.7K, an inhibitor of TNF-alpha functions (Li et al., 1997). ERAP1 and ERAP2 (endoplasmic reticulum aminopeptidase 1 and 2) encode aminopeptidases involved with cleavage of viral and bacterial proteins into smaller fragments within the endoplasmic reticulum and transport to the cell surface (Falk and Rötzschke, 2002). At the surface, these proteins can then cleave cytokine receptors and reduce their ability to transmit signals into the cell, thus affecting inflammatory responses. Small peptides produced by cleavage of bacterial or viral proteins can also attach to MHC class I molecules in the endoplasmic reticulum and be exported to the cell surface and displayed to the immune system. If the MHC class I/peptide complex is recognized as foreign by the immune system, the infected cell will undergo apoptosis. Polymorphisms in ERAP1, 2 as well as MHC class I have been associated with natural selection among human populations for resistance to intracellular pathogens such as HIV-1 (Cagliani et al., 2010).

The fact that these potential candidate genes were involved in some manner with processes such as immune regulation and apoptosis makes them highly relevant to what is known about FIPV pathogenesis (reviewed Pedersen, 2009). However, two crucial questions regarding FIP pathogenesis remain to be explained: (1) why does FIPV target a type of macrophage that has a peculiar tropism for venules in the omentum, serosal surfaces, meninges, and uveal tract, and (2) why do these macrophages fail to undergo pathogen-induced apoptosis either by intracellular or extracellular signaling? Targeting of a specific macrophage subset is most likely associated with the way the virus and host cell membrane interact, while polymorphisms within any or all of the five candidate genes identified in the present study are more likely to be associated with why infected macrophages remain persistently infected. All five genes relate to a programmed cell death in the face of viral infection through inhibition of viral replication, cytokine interactions, or MHC class I/peptide complex induced apoptosis.

Although there were several moderately strong SNP/genome regions associated with FIP susceptibility, 100% concordance with the FIP phenotype was not detected for any individual SNP or region. This could be expected in situations where multiple genes contribute only small to moderate degrees of risk either independently or collectively (epistasis) and are not inherited uniformly by every individual. Complex traits must be also interpreted with clear guidelines in mind. Three such guidelines for dissecting complex genetic traits have been advocated by Lander and Kruglyak (1995) and include: (1) scan entire genomes with a dense collection of markers, (2) calculate a linkage statistic at each marker position, and (3) identify regions where the association deviates significantly from what would be expected from random assortment. However, this recipe does not take into account population size, degree of heritability, and numerous non-genetic covariates (Thornton-Wells et al., 2004). Lander and Kruglyak (1995) also stressed the importance of reporting negative or weak findings when the evidence clearly supports a negative or weak association, while confirmatory linkage requires the same result to be found in a second study, preferably done by an independent group of investigators and at an adjusted probability of ≤0.01.

It is hoped that the present study can be expanded to the point where genotypes linked to FIP susceptibility can be confirmed. Confirmatory studies may require greater numbers of case and control cats, much denser SNP arrays, and even whole genome sequencing. Confirmed FIP susceptibility genotypes could be used to breed cats more resistant to FIP. Until this time is reached, the best genetic advice that can be given to pedigreed cat breeders is not to breed healthy toms and queens that produce offspring dying from FIP, or to breed offspring of parents that died subsequently of FIP.