INTRODUCTION
Progressive retinal atrophy (PRA) is a comprehensive term encompassing several inherited retinal diseases which all progress over time and eventually leads to vision loss in affected dogs. PRA in dogs has been reported to be genetically similar to the retinitis pigmentosa (RP) disease occurs in human[1]. Similar to RP, PRA shows genetically heterogeneous and phenotypic similarities [2]. The first significant clinical sign of this disease is night blindness, followed by total eyesight loss. [1]. In PRA disease, various ophthalmoscopic complication mostly takes place in eye such as retinal vascular degeneration, retinal detachment, loss of pigmentation, hyperreflectivity in tapetal area, and optic atrophy [3].
PRA disease has been first recognized during the early 20th century after an impulsive retinal deterioration due to PRA was reported in setter family dogs [4–6]. At present, several forms of PRA have been described in more than 100 different breeds of dogs [1], although many of those exhibit similar ophthalmologic signs; the difference between breeds is only in the rate of progression, aetiology, and age. The majority of hereditary PRA diseases in dog are autosomal recessive, while some are X-linked forms or autosomal dominant [1,5,7]. So far, researchers have identified as many as 31 mutations in 24 genes, as underlying cause for canine retinal disease [8]. The methods to identify susceptible genetic variants and genes discovery in canine have seen considerable improvement over the past 20 years [8,9]. Unfortunately, the real treatment for PRA disease in canine still does not exist. Hence, identification of PRA associated genetic variations using deferent approach like Single nucleotide polymorphism (SNP) microarrays and next generation sequencing, is essential for finding disease associated candidate genes.
Genome-wide association studies (GWAS) is considered worldwide as an influential tool for identifying causative loci and genes associated with genetic disease. Nevertheless, conducting robust GWAS analysis with the aim of detailed biological understanding is challenging due to certain factors like large number of samples, stringent statistical thresholds to account for multiple testing, small effect sizes of variants resulting in very few reaching the stringent threshold levels, and linkage disequilibrium [10–12]. GWAS cannot describe the fact that pathways and cluster of genes are connected with each other to implement a specific molecular function [13]. Additionally, GWAS possibly will not able to fully capture the effect of multi-allelic quantitative trait locus because of markers having two allelic variants [14]. Hence, the identification of genetic loci through GWAS only might offer inadequate insight into the polygenic disease like PRA.
Recommended resolution for addressing those issues mentioned earlier is to perform GWAS accompanied by gene-set and pathway analyses. The integration of GWAS, gene-set analysis, pathway and network analyses might help to overcome the limitations of GWAS and provide insights into the biological pathways affecting the particular trait. Such approaches have been used in humans [15], and dairy cattle [13,14]. Though several GWAS have been performed for PRA disease in canine [1,5,7,16], those studies have resulted in only few markers and have not exploited the use of pathway analyses for understanding the underlying biology of PRA. Therefore, aim of the present study was to carry out GWAS analysis for finding out the sensitive loci and genes associated to PRA disease. Additionally, we have identified regulatory gene ontology (GO) terms and pathways which are associated to canine PRA disease for understanding the biological functions regulating this disease.
MATERIALS AND METHODS
Whole blood sample was collected from a total of 129 dogs (101 ophthalmoscopically normal dogs, 28 dogs diagnosed with PRA) of 18 different breeds (Table 1) and stabilized by ethylenediaminetetraacetic acid. We followed an appropriate instruction articulated by The Institutional Animal Care and Use Committee (IACUC) of National Institute of Animal Science (NIAS). These dogs were diagnosed individually by veterinary ophthalmologists from the division of animal disease and health, NIAS, Korea. Genomic DNA was isolated from blood sample followed by wizard® genomic DNA purification kit Protocol (Promega, Madison, WI, USA).
Genome-wide association analysis was performed in 28 cases and 101 controls using the Illumina Canine HD Bead-Chip, which genotypes 173,662 SNPs evenly spaced across the genome. The GWAS analysis was performed with PLINK software (version 1.9). Prior to the analysis, quality control (QC) filtering with a minor allele frequency < 5%, low genotyping call rate < 90%, missing genotype calls > 10%, Hardy–Weinberg equilibrium at p < 1E-06 were carried out. The final genotyping call rate was 98.5%. Following QC filtering, 135,553 SNPs and 128 animals remained for further analyses. The PRA phenotype used in this analysis was discrete. A case–control GWAS was performed to test the association between disease trait and allele in each SNP with the -- assoc option of PLINK software. We generated a total of 20 PC; the eigenvalues of all the PCs were fit as covariance to account for population stratification. The following statistical model equation was used in this GWAS study.
Where, y is experiential phenotype, x and z are known design matrices, b is a vector of fixed effects (sex, age, and breed), u is the additive genetic effect for each SNP and e is vector of residual terms. Manhattan and quantile-quantile (Q-Q) plots were generated with CMplot [17]. Bonferroni correction was deemed as too stringent; therefore, the significance level was set at nominal p< 0.01 to identify putative markers.
To identify the pathways and genes enriched for the PRA disease in dogs, we performed GO and pathway analysis following the process discussed by Dadousis et al. [13]. The SNPs (nominal p < 0.01) from GWAS analysis were used for the analysis. The markers were annotated to the nearby genes using the SnpEff software version 4.3 [18]. Then, the genes overrepresentation in every GO terms and pathway were tested using a Fisher’s exact test, which was conducted with the web-based functional enrichment analysis tool, the Database for Annotation, Visualization, and Integrated Discovery (DAVID, https://david.abcc.ncifcrf.gov/) [19]. The GO [20] and Kyoto Encyclopedia of Genes and Genomes (KEGG) [21] databases were used to describe the GO terms and pathways related to the enriched genes. The GO terms/KEGG category pathways enriched by more than 10 and less than 1,000 genes were examined in order to narrow down the functional categories [22].
RESULTS
After QC, 4,108, 3,418 and 30,583 markers were removed due to missing genotype data, HWE (p < 0.000001) and MAF < 5%. A total of 135,553 markers were finally passed the QC and genotyped in the GWAS analysis. We performed case–control based association analysis with PRA. Manhattan plots were constructed based on the association p-values to find major loci associated with PRA disease. The physical location of SNPs in chromosome was plotted against −log10p-values (Fig. 1A). Even though the number of cases was not sufficient for effective GWAS analysis, our findings provide some insight into major loci underlying PRA. However, no significant loci were found from the association analysis. A Q-Q plot of observed p-values was plotted against expected p-values for PRA is shown in Fig. 1B. There were no significant deviation (λ value = 1.01) which indicates absence of significant SNPs. The lack of highly significant loci might be due the inclusion of large number of breed (n = 18) and lower sample number.
GWAS was supplemented by a functional enrichment analysis for detecting significant GOs and pathways associated to PRA. Among the total 135,553 SNPs used in this study, 91,779 SNPs were mapped to the genes. Based on reduced GWAS threshold of (p < 0.01) 1,114 SNPs were deemed to be associated with PRA (Supplementary Table S1). These SNPs were annotated to assign 640 unique genes, (Supplementary Table S1). Subsequently, the gene list was uploaded to DAVID for GO enrichment and pathway analysis. A total of 76 GO terms and 11 KEGG pathways were found to be enriched (Supplementary Table S2 and Table S3). Out of the total enriched GO terms, 35 significantly enriched (p < 0.05) GO terms with at least 4 genes are presented here. There were 26 biological processes related to cell adhesion, cell migration, protein kinase C, regulation of cell signaling, phosphorylation, ion transport, cell growth, cell proliferation, axon guidance, DNA replication, apoptosis and startle response; 4 cellular component GO terms including focal adhesion, endoplasmic reticulum, basal lamina, and early endosome membrane protein, and 4 molecular function GO including calcium ion binding, adenosine triphosphate (ATP) binding, protein kinase activity, potassium channel regulator activity, and heterodimerization activity (Table 2). ATP binding (GO:0005524) was the most significant with the highest number of genes (49) among all GO terms. Meanwhile, homophilic cell adhesion via plasma membrane adhesion molecules (GO:0007156) was the most significantly enriched biological process for PRA trait (p = 8.32E-06).
5 KEGG pathways were also significantly (p < 0.05) enriched for PRA trait including genes involved in focal adhesion, axon guidance, cyclic guanosine monophosphate (cGMP)-protein kinase G (PKG) signaling pathway, gastric acid secretion, and platelet activation (Table 3). Focal adhesion pathway was the most significantly enriched pathway with PRA associated genes. Besides, it was enriched with the highest number of genes (14) from the gene list uploaded for pathway analysis.
DISCUSSION
The present study investigated functional enrichment based association analysis, a post-GWAS method to understand PRA disease biology in dogs. These approaches were not previously reported for revealing the functional associations among genes associated to PRA disease. The pathway-based association analyses effectively check the cluster of gene enriched in every pathway for providing the knowledge on the pathway-level divergence [10,23]. Moreover, it offers us different understandings about genes as well as their interaction with each other. Thus, these approaches have been quite successful in acquiring new information regarding a particular disease trait.
In this study, no locus was identified at the genome wide and suggestive significant levels in our GWAS data set. Therefore, nominal p < 0.01 from GWAS analysis was used to subset for SNPs for further functional enrichment analysis. In our gene-set enrichment analysis, the most significant biological process (GO:0007156) enriched for the PRA trait was “Homophilic cell adhesion via plasma membrane adhesion molecules”. This biological process is enriched with 13 significant genes (CDH12, CDH8, CDH13, CRTAM, CDH9, PTPRM, FAT3, ROBO1, DSC3, CDHR3, ROBO2, PCDH15, and CDH10) and all of them except two (ROBO2 and PTPRM) encode proteins that belong to cadherin superfamily. Interestingly, some of these genes, specially belong to cadherin family are also enriched in calcium ion binding molecular function (GO:0005509). It is widely known that cadherins, a class of transmermbrane protein depends on calcium ion binding to function, cell-cell adhesion [24]. Notably, these cadherin proteins have been significantly associated with PRA because of their involvement in retinal detachment and proliferative vitreoretinopathy [24,25]. The retinal separation was reported to influence the neural retinal atrophy by cutting choroidal blood supply to the outer segments of retinal pigment epithelium [26,27]. Moreover, cadherin proteins are found to be an important adhesion molecule in retinal pigmental epithelial cells [28]. In addition, another adhesion related biological process for PRA disease was GO:0005925 term “focal adhesion” shared one gene (CDH13) with other GO term “homophilic cell adhesion via plasma membrane adhesion molecules”. The gene (CDH13) encodes cadherin13 protein, which was reported to be extremely associated with cancer disease including retina cancer, lung cancer and breast cancer [29]. The disease relevance of cadherin13 studies found that the expression of this protein was invariable in retinoblastoma eye specimen of human than healthy eye [29]. Moreover, the association of CDH13 with Ca (2+)-dependent cell-cell adhesion molecule suggested distinctive connection in retinal atrophy [30].
Apoptosis or programmed cell death process is implicated in several neurodegenerative diseases, including the inherited retinal degenerations knows as RP diseases [31,32]. Marigo [33] identified two apoptotic pathways in undergoing cell death photoreceptors of autosomal recessive RP mouse model . They even suggested the use of anti-apoptotic strategies against such apoptosis may open new perception about innovative therapeutic strategy for treatment of retinal deterioration [33]. In this study we found the GO term apoptotic process (GO:0006915) was to be significantly enriched with PRA associated genes. Among the genes included in this specific GO, deleted in colorectal carcinoma (DCC) gene was significantly enriched in another biological process, axon guidance (GO:0007411). The DCC gene has a key role in the dual functionality of netrin-1 in axon guidance biological process [34]. Moreover, this axon guidance was detected to be associated with axonal degradation of retinal ganglion cells, which is directly linked to PRA disease [35].
The maximum number of genes was found to be enriched in GO term “ATP binding” (GO:000524). Specially, retina-specific gene ABCA4 and Na+/K+ transporting gene ATP1A4 is of particular interest. The ABCA4 gene encodes ATP-binding cassette, subfamily A, member 4, protein which was described to be an essential protein for human retinal degeneration disease. Moreover, this gene was found as an causal gene associated to generalized PRA disease in different breeds of dogs [36]. Besides, ATP1A, gene comes under Na+/K+ transporting ATPase alpha 1 family, has been a causative gene for autosomal recessive cone-rod dystrophy (CORD8) disease [37]. Moreover, this Na+/K+ ATPase was exhibited recently to be important in photoreceptor cell [38]. These above descriptions suggest an indirect link to the ATP binding term and PRA disease.
Another important GO related to “endoplasmic reticulum” was significantly enriched for PRA. Endoplasmic reticulum (ER) has vital contribution in membrane protein folding [39]. In cells, ER generally goes under stress upon failing to match the protein folding capacity demand [39,40]. Interestingly, this ER stress was found to be elevated in multiple animal model of RP disease [41]. The key role of ER in retinal neurodegeneration has already been described [40–43]. ER leads the degradation of pigment in photoreceptor cells which ultimately causes blindness [39]. Now, it can be hypothesized that GO term “endoplasmic reticulum” is strongly associated with the retinal disease including PRA.
Among the 5 significant KEGG pathway enriched for PRA, Focal adhesion and Axon guidance were found to be common in GO and KEGG categories. The gene sets enriched for Focal adhesion and Axon guidance varies between GO and pathway analysis. Focal adhesion pathway (Cfa04510) was the most significantly enriched pathway for the PRA trait. Focal adhesion kinase, member of the focal adhesion pathway, has been associated in retinal angiogenesis, which is a major feature of eye diseases that can lead to blindness [44–46]. Hence, a causal relationship among focal adhesion and retinal degeneration can be hypothesized (part of the PRA disease). In addition, among the included genes in focal adhesion pathway, AKT2 is of particular interest. This gene was common in another two pathways in our result such as cGMP-PKG signaling pathway and platelet activation pathway. This gene is closely related with serine/threonine protein kinases (AKT1, AKT2, and AKT3) and was reported to have contribution in light stress-induced retinal photoreceptor degeneration [47]. Furthermore, axon guidance pathway was significantly enriched with 10 genes including ABLIM1, DCC, EPHA4, CXCR4, ROBO1, NTN4, SEMA3C, ROBO2, EPHB1, and RASA1. The relationship between axon guidance pathway and degeneration of retina (causing progressive vision loss), was described above.
Our pathway analysis revealed that a vital signaling pathway such as cGMP-PKG signalling pathway (Cfa04022) was amid the enriched pathways in our study. The cGMP- dependent (PKG) signaling pathway is known to be intricately linked to the hereditary retinal degeneration such as RP [48]. The importance of cGMP signaling pathway for neurodegenerative disease retinal degeneration was recently reviewed [49]. A previous study showed that cGMP-dependent PKG activation is a hallmark in photoreceptor degeneration in human homologus mouse [48]. Moreover, the same authors suggested to target PKG in cGMP dependent PKG pathway as a unique therapy for retinal degeneration disease [48].
Additionally, the gastric acid secretion pathway was significantly enriched for PRA trait. The ATP1A1 (ATPase Na+/K+ Transporting Subunit Alpha 1) gene encodes the protein which was reported to be connected to photoreceptor cell, which is a part of retina [38]. There is reported relation between Gastric acid secretion pathway and PRA trait. Last of all, the platelet activation (Cfa04611) was significantly enriched pathway for PRA trait. Among the most of the genes in this pathway were shared with focal adhesion pathway and cGMP-PKG signaling pathway. Moreover, there is an evidence of association between platelet activating factor and pathogenesis of age-related eye disease [50]. Hence, it suggests an unintended link between platelet activation pathway and PRA trait.
Seemingly, aside from ontologies and pathway that could be directly allied with PRA disease (e.g., GO term apoptotic process and cGMP-PKG signaling pathway); ontologies and pathway not strictly associated to the traits of our interest were also identified (negative regulation of T cell proliferation, platelet activation). Moreover, it is obvious that some of our result may need further validation. It is not noteworthy to mention here that the GO and KEGG pathways of canine which are accessible for public are remain scarce and not adequately elucidated. Therefore, an RNA-Seq analysis may be more useful for identifying critical pathways and genes influencing PRA.
CONCLUSION
The present study is the first to report about GWAS analysis together pathway and gene-set enrichment analysis for PRA in dogs. Our results showed new insight into gene sets and pathways associated to PRA suggesting that this disease is the influence of combined consequence of a variety of genes work together in a pathway. Out of many GO term for PRA, in particular calcium ion binding, endoplasmic reticulum, homophilic cell adhesion via plasma membrane adhesion molecules, apoptotic process were found to be strongly related to PRA trait. Among the pathways, two pathways (focal adhesion and cGMP-PKG signaling) were observed to be closely related to PRA. Thus, the genes associated in highlighted GOs and pathways identified may be valuable candidate for further studies regarding the underlying cause of PRA disease. However, additional gene expression studies are needed to identify rare genes and regulatory SNPs which play significant role behind the PRA disease.