Natural Genetic Variation Screen in Drosophila Identifies Wnt Signaling, Mitochondrial Metabolism, and Redox Homeostasis Genes as Modifiers of Apoptosis.

Apoptosis is the primary cause of degeneration in a number of neuronal, muscular, and metabolic disorders. These diseases are subject to a great deal of phenotypic heterogeneity in patient populations, primarily due to differences in genetic variation between individuals. This creates a barrier to effective diagnosis and treatment. Understanding how genetic variation influences apoptosis could lead to the development of new therapeutics and better personalized treatment approaches. In this study, we examine the impact of the natural genetic variation in the Drosophila Genetic Reference Panel (DGRP) on two models of apoptosis-induced retinal degeneration: overexpression of p53 or reaper (rpr). We identify a number of known apoptotic, neural, and developmental genes as candidate modifiers of degeneration. We also use Gene Set Enrichment Analysis (GSEA) to identify pathways that harbor genetic variation that impact these apoptosis models, including Wnt signaling, mitochondrial metabolism, and redox homeostasis. Finally, we demonstrate that many of these candidates have a functional effect on apoptosis and degeneration. These studies provide a number of avenues for modifying genes and pathways of apoptosis-related disease.

Phenotypic heterogeneity is the driving force behind the Precision Medicine Initiative (Scriver and Waters 1999; Nadeau 2001; Queitsch ; Gallati 2014). Patients suffering from the same genetic disorders can carry identical causal mutations but often display wildly variable phenotypes and symptom severity. A large part of this variation is due to inter-individual differences in genetic background, including silent cryptic genetic variation that is revealed upon disease or stress (Queitsch ; Chow 2016). Understanding the role of this variation and the genes or pathways which modify disease will lead to improved personalized therapeutic predictions, strategies, and diagnostics.

One process implicated in many genetic disorders is programmed cell death or apoptosis (Elmore 2007; Sano and Reed 2013; Kurtishi ). During normal development and tissue turnover, cells can receive both internal and external signals that trigger a programmed response which eventually results in the death of the cell (Elmore 2007). Because cell death is essential to cellular, tissue, and organismal homeostasis, disruption of apoptosis pathways can be catastrophic. Inhibition of apoptosis is an important step in transformation and cancer, while excess apoptosis, often activated by chronic cellular stress, is a primary cause of degeneration in different neuronal, retinal, muscular-skeletal, and metabolic diseases (Mattson 2000; Elmore 2007; Ouyang ). As a result, an important area of therapeutic development is focused on targeting apoptosis without disrupting normal tissue homeostasis (Elmore 2007). Our previous work demonstrated that hereditary variation in apoptotic genes is associated with phenotypic variation in a model of retinal degeneration, suggesting that modifiers of apoptosis could serve as drug targets in degenerative diseases (Chow ).

Model organism tools, such as the Drosophila Genetic Reference Panel (DGRP), enable the study of the impact of natural genetic variation on diseases and related pathways. The DGRP is a collection of ∼200 isogenic strains derived from a wild population, such that each strain represents one wild-derived genome (Mackay ). The variation in the DGRP is well tolerated under healthy, non-disease conditions and allows for the identification of genetic polymorphisms that are associated with phenotypic variation in models of human disease (Chow and Reiter 2017). Importantly, the availability of full-genome sequence for these strains allows for genome-wide association analyses that link quantitative phenotypes with genetic variation and modifier genes.

In this study, we report the results of natural variation screens of - () and -induced apoptosis (Figure 1). Overexpression of either of these genes leads to massive apoptotic activation (Hay ; Jin ). While there is a great deal of overlap between these pathways, they can each activate apoptosis independently. p53 is stabilized in response to DNA damage and initiates apoptosis by transcriptionally activating the inhibitor of apoptosis (IAP) inhibitors , , and (Mollereau and Ma 2014) (Figure 1). P53 can also increase apoptosis by activating JNK signaling and stabilizing the IAP inhibitor Hid (Shklover ). is activated transcriptionally by either p53 or the JNK signaling cascade, which is induced downstream of oxidative, ER, and other cellular stresses (Kanda and Miura 2004; Shlevkov and Morata 2012) (Figure 1). We designed this study to identify genetic modifiers of general apoptosis induced by any cellular pathway, including modifiers that are specific to stress-induced, p53-independent pathways or specific to canonical p53-dependent pathways.

Figure 1

Activation of apoptosis through and -associated pathways. Apoptosis is primarily initiated through either p53 or Jun-induced (JNK) transcriptional activation of the Inhibitor of Apoptosis (IAP, in Drosophila Diap) inhibitors , and . While p53 is primarily activated by DNA damage and disruption of the cell cycle, JNK signaling is activated downstream of cellular stress, such as endoplasmic reticulum (ER) stress, through Ire1 and Cdk5. ER stress occurs when misfolding proteins, like the rhodopsin mutant Rh1, accumulate in the ER (Chow ). Expression of , , and leads to inhibition of Diap1, releasing the inhibition on initiator caspases and allows for the activation of effector caspases and downstream apoptosis. Models used in this or previous studies of retinal degeneration in the DGRP are indicated in white.

We observed substantial phenotypic variation across the DGRP for both - and -induced apoptosis. Using genome-wide association analysis, we identified a number of modifying pathways and genes, several of which have known roles in cell death pathways, neuronal development, neuromuscular diseases, and cancer. Using systems biology approaches, we also identified Wnt signaling, mitochondrial redox homeostasis, and protein ubiquitination/degradation as possible modifiers of apoptosis. Finally, we confirmed that reduction in the expression of many of these candidate modifier genes significantly alters degeneration. Our findings highlight several exciting new areas of study for apoptotic modifiers, as well as a role for stress-induced cell death in the regulation of degenerative disorders.

Methods

Fly stocks and maintenance

Flies were raised at room temperature (∼20-22°) on a diet based on the Bloomington Stock Center standard medium with malt. This diet consists of cornmeal, yeast, malt, and corn syrup and followed the recipe and proportions of the Bloomington cornmeal diet without the addition of soy flour. In this study, p53 is driven by the GAL4/UAS system and the strain contains a GMR-GAL4 transgene and a UAS-p53 transgene (GMR > ). expression is driven directly by the GMR promoter sequences from a single transgene (GMR-rpr). The strains containing GMR-GAL4 and UAS-p53 or GMR-rpr on the second chromosome have been previously described (Hay ; Jin ). In both cases, the apoptotic gene ( or ) is overexpressed specifically in the developing Drosophila eye under the control of the GMR promotor. This leads to excess apoptosis in the eye imaginal disc and ultimately a small, degenerate eye once the adult ecloses. These are referred to as the apoptotic models throughout the manuscript. 204 strains from the DGRP were used for the GMR > study (Table S1) and 202 were used for the GMR-rpr study (Table S2). In both cases virgin females carrying one of the apoptosis models were crossed to males of the DGRP strains. F1 progeny carrying GMR > or GMR-rpr were collected and scored for eye size. As excess apoptosis leads to the degenerate eye phenotype, eye size was used as a proxy for cell death and variation in eye size as a proxy for variation in the degree to which apoptosis was activated during development. The following RNAi and control strains are from the Bloomington Stock Center: swim RNAi (55961), RNAi (57560), RNAi (32967), αMan1a RNAi (64944), RNAi (62153), RNAi (53345), RNAi (38998), RNAi (34320), RNAi (34012), RNAi (44483), RNAi (33645), shab RNAi (55682), RNAi (64671), Cyt-c-P RNAi (64898), Ir40A (57566), RNAi (61934), control attP40 (36304), and control attP2 (36303). To test RNAi in the GMR-rpr model, a GMR-GAL4 transgene was introduced into the strain through genetic crosses. While this effectively changes the genetic background of the GMR-rpr model, we control for this as much as possible by crossing this line to the attP40 and attP2 lines to generate the no RNAi controls.

Eye size imaging

For eye images, adult females of the necessary genotypes were sorted by phenotype and sex under CO2 anesthesia. As cell death in these models occurs during development, eye size is stable through adulthood. Flies were aged at least 2 days to ensure they are fully mature, but not more than 7 days to prevent accumulated lethality due to a leaky GAL4 driver or UAS promotor, then flash frozen on dry ice. Left eyes were imaged for all measurements to maintain consistency. 10-15 eyes per strain were imaged using a Leica EC3 camera (Leica Microsystems) mounted on a FlyStuff Trinocular Microscope (Genesee Scientific) at 3X magnification. Camera settings were as follows: brightness-70%, gamma-0.70, saturation-140.00, EC3-14140014, configuration-last used, capture format-2048x1536, live format-1024x768, shading-none, sharpening-low. Eye area was measured in ImageJ as previously described (Chow ). Two-dimensional eye area was measured as a proxy for eye size for each individual of each line. The outlines of the eyes were carefully traced using the freehand selection tool on ImageJ. We then used ImageJ to calculate the area in pixels enclosed by this selection and recorded that measurement as eye size in that individual. Average eye size for each line or transgenic strain was determined from 10-15 individual measurements.

Phenotypic analysis and genome-wide association

For each trans-heterozygous DGRP/GMR > or GMR-rpr line, eyes from 10-15 individual females were imaged and measured. For each model, a one-way ANOVA (R software) was used to determine if there was an effect of genetic background/strain on eye size. A genome-wide association (GWA) analyses was performed to identify genomic variants that are significantly associated with differences in eye size. Mean eye area was used for the GWA. GWA was performed as previously described (Chow ). DGRP genotypes were downloaded from the website, http://dgrp.gnets.ncsu.edu/. Variants were filtered for minor allele frequency (≥ 0.05), and non-biallelic sites were removed. A total of 1,967,719 variants for and 1,962,205 variants for were included in the analysis. Mean eye size for 204 F1 DGRP/GMR > strains (representing 2953 flies) or 202 F1 DGRP/GMR-rpr strains (representing 2987 flies) were regressed on each SNP. To account for cryptic relatedness (He ; Huang ), GEMMA (v. 0.94) (Zhou and Stephens 2012) was used to both estimate a centered genetic relatedness matrix and perform association tests using the following linear mixed model (LMM):where, as described and adapted from Zhou and Stephens 2012, y is the n-vector of mean eye sizes for the n lines, α is the intercept, x is the n-vector of marker genotypes, β is the effect size of the marker. u is a n x n matrix of random effects with a multivariate normal distribution (MVN_n) that depends on λ, the ratio between the two variance components, τ^(-1), the variance of residuals errors, and where the covariance matrix is informed by K, the calculated n x n marker-based relatedness matrix. K accounts for all pairwise non-random sharing of genetic material among lines. ϵ, is a n-vector of residual errors, with a multivariate normal distribution that depends on τ^(-1) and I_n, the identity matrix. Quantile-quantile (qq) plots demonstrate an appropriate fit to the LMM (Figure S1). Genes were identified from SNP coordinates using the BDGP R54/dm3 genome build. A SNP was assigned to a gene if it was +/− 1 kb from a gene body.

Correlation Analysis

A Pearson Correlation test was performed to compare mean eye size between F1 DGRP/GMR > strains and F1 DGRP/GMR-rpr strains of the same DGRP backgrounds (i.e., the F1 of the RAL 227 crossed to GMR > compared to the F1 of the RAL 227 crossed to GMR-rpr). The same test was used to compare the F1 DGRP/GMR-rpr and F1 DGRP/GMR > Rh1 (Chow ), as well as F1 DGRP/GMR > and F1 DGRP/GMR > Rh1. Statistics were calculated using using R software.

RNAi Validation

Virgin females of either the GMR-GAL4; GMR-rpr or GMR > genotype were crossed to males carrying UAS-RNAi constructs targeting candidate modifier genes. Eye size of F1 progeny expressing both the apoptotic model and the RNAi construct was measured as described above. The eyes of 10-15 females were imaged and measured. Eye size from RNAi-carrying strains were compared directly to genetically matched attP40 or attP2 controls using a Dunnett’s multiple comparisons test.

Bioinformatics Analysis

Genetic polymorphisms within an annotated gene were associated with that particular gene. Polymorphisms located outside of an annotated gene were associated with the closest annotated gene within 1 Kb. Polymorphisms located more than 1 Kb away from an annotated gene were labeled “intergenic” and not further considered. Information about candidate genes and their human orthologs was gathered from a number of databases including Flymine, Flybase, OMIM, and NCBI. Genetic interaction maps were generated using the GeneMANIA plugin on Cytoscape (version 3.6.1) (Shannon ; Montojo ). Gene Set Enrichment Analysis (GSEA) was run on all polymorphisms within 1 kb of a gene to generate a rank-list of genes based on their enrichment for significantly associated polymorphisms (see figshare for code). For GSEA analysis, polymorphisms within 1 kb of more than 1 gene were assigned to one gene based on an a priori list of exon, UTR, intron, and upstream or downstream. Genes were then assigned to GO categories. Using the entire ranked gene list instead of limiting the list by corrected p-value cut-off (Dyer ), GSEA determines whether members of a given set of genes belonging to a specific GO category are randomly distributed throughout the ranked list, are found primarily at the top of the ranked list, or are found primarily at the bottom of the ranked list. GO categories enriched at the top of the list functionally describe the phenotype of the gene set. Calculation of enrichment score was performed as described (Subramanian ). Categories with ES scores > 0 (enriched for associated genes with low p-values, ES scores ≤ 0 are unenriched), gene number > 3 (only see pathways represented by multiple variants), and corrected p-values <0.05 (significant associations) were included in the final output.

Reagent and Data Availability

Strains and stocks are available upon request. Genomic sequence for the DGRP is available at http://dgrp.gnets.ncsu.edu/. Code and related guides for the GSEA analysis and supplemental material available at figshare: https://doi.org/10.25387/g3.9808379.

Results and Discussion

rpr- and p53-induced apoptosis is dependent on genetic background

We used the Drosophila eye to model apoptosis. Expression of either or in the ommatidial array of the developing eye imaginal disc results in massive cell death and smaller, rough adult eyes (Hay ; Jin ). The model is induced by direct drive of the GMR promoter (GMR-rpr) on a second chromosome balancer. The model is induced using the GAL4/UAS system, where GMR-GAL4 drives expression of UAS-p53 (GMR > ). Importantly, in both of these models, adult eye size is an easily scorable, quantitative proxy for levels of apoptosis. The lines described serve as the donor strains (GMR > p53/CyO or GMR-rpr,CyO/sna) that we crossed to each DGRP strain. Females from the donor strains were crossed with males of each of 204 or 202 DGRP strains to generate F1 progeny that overexpressed or , respectively, in the eye disc. The progeny received 50% of their genome from the maternal donor strain and 50% from the paternal DGRP strain. Therefore, we are measuring the dominant effect of the DGRP background on the or retinal phenotype. This cross design is similar to a study of ER stress-induced degeneration (Chow ) and a model of Parkinson’s Disease (Lavoy ) we previously reported. The GMR promoter has been used in several previous publications (Chow ; He ) to drive overexpression of a UAS-transgene in the developing Drosophila eye imaginal disc. In all of these cases, the top candidates have been substantially different and highly dependent on the UAS-transgene. This suggests that the protein being overexpressed is more influential on the phenotype than the GMR-GAL4 transgene. We examined eye size in the F1 progeny to determine the average eye size in individual genetic backgrounds (Figure 2A-D). As previous work has demonstrated a lack of correlation between variation in body size and variation in eye size in a model of developmental eye degeneration, there was no need to correct for overall body size (Chow ).

Figure 2

Apoptosis levels vary across genetic background in and models of apoptosis-induced degeneration. Apoptosis induced by overexpression of (A) or (B), as measured by adult eye size, varies across different genetic backgrounds. DGRP strains are arranged along the X-axis from smallest to largest. Each point represents the median eye size of an individual DGRP strain, while the error bars represent the standard deviation of eye size for that strain. Representative images of GMR > eyes (C) or GMR-rpr eyes (D) in different DGRP backgrounds demonstrate the phenotypic variation quantified in panels A and B.

We first tested the effect of sex on apoptosis in a pilot study. We measured eye area in at least ten females and ten males from eight different DGRP strains crossed to either the or model. Eye size is positively correlated between males and females (Figure S2). Because variation is greater in females (Figure S2), we elected to focus on female eye size for the remainder of our analysis.

We found a significant effect of genetic background on eye size in the GMR > model (P < 2.2 × 10−16) (Figure 2A,C, Table S1). Average eye size measured in pixels on ImageJ ranged from 10542 pixels (RAL812) to 17835 pixels (RAL374) (Table S1). Similarly, we found a significant effect of genetic background on eye size in the GMR-rpr model (P < 2.2 × 10−16), with median eye size ranging from 7957 pixels (RAL83) to 16884 pixels (RAL304) (Figure 2B,D, Table S2). For both the GMR > and the GMR-rpr models, the variation in eye size within individual DGRP strains is substantially smaller than the variation observed between DGRP strains (Figure 2A-B, Table S1-S2).

We noted that the range in average eye size for the GMR-rpr model (8927 pixels) is greater than that seen in the GMR > model (7293 pixels). This could be due to the greater involvement of in a variety of stress-induced, -independent apoptotic pathways (Shlevkov and Morata 2012). Alternatively, it is possible that variation in -associated pathways is simply less well-tolerated than in -associated pathways. It is also possible that the DGRP simply carries more variation affecting the GMR-rpr model than GMR > model. Another possibility is that the differences in genetic background between these two models is driving differences in the way they respond to modifier variation in the DGRP. Finally, it is possible that the difference in range we observe are not statistically significant and is simply a stochastic artifact.

We observed qualitative differences between the apoptotic models, with flies expressing the GMR > model displaying a teardrop-shaped eye (Figure 2C) and flies expressing the GMR-rpr model displaying a rounder eye (Figure 2D). These qualitative shapes were not subject to effects of genetic variation. The differences in eye shape noted between GMR > and GMR-rpr, however, could be indicative of differences in the mechanisms by which apoptosis and degeneration progress in these two models. Alternatively, this could be evidence of the technical differences in the two models, since is driven by the GAL4/UAS system and is driven directly by the GMR promotor. We saw no accumulation of necrotic tissue in strains experiencing severe degeneration, nor did we note obvious differences in pigmentation (Figure 2C,D). Eyes from all strains maintained the rough-eye phenotype that is characteristic of or -induced degeneration, indicating that while modifying variation may reduce the amount of cell death in the eye imaginal disc, it cannot fully rescue the degenerative phenotype.

Eye size in rpr-expressing DGRP lines correlates with eye size in other models of degeneration

Because canonical p53 signaling activates the expression of , we expected high correlation in apoptosis levels and eye size between these models (Shlevkov and Morata 2012; Mollereau and Ma 2014). Indeed, there is a significant positive correlation in eye size between DGRP strains expressing GMR > and GMR-rpr (r = 0.19, P = 0.0071) (Figure 3A). In a previous study, we examined the impact of genetic variation on a model of retinitis pigmentosa (RP) and ER stress-induced apoptosis (Chow ). In that study, we found that the degeneration induced by overexpression of a misfolded protein (Rh1) in the developing eye imaginal disc is modified by a number of genes involved in apoptosis (Chow ). This is to be expected, as the primary cause of degeneration in this model is JNK-hid/grim/-mediated cell death (Figure 1) (Kang ). Consistent with this mechanism of Rh1-induced degeneration, we found a significant correlation in eye size between the Rh1 and models (r = 0.25, P = 0.001, Figure 3B). In contrast, we see no correlation between the Rh1 and models of apoptosis (r = 0.12, P = 0.13) (Figure 3C). These results suggest that there is shared genetic architecture between Rh1 and -mediated apoptosis and degeneration that is independent from that shared between and .

Figure 3

Eye size is correlated between GMR-rpr and both GMR > and GMR > Rh1 models of degeneration. Correlation in mean eye size between the GMR-rpr, GMR > , and GMR > Rh1 models across the DGRP. A. Eye size is significantly correlated in the same DGRP strains expressing GMR-rpr and GMR > (r = 0.19, P = 0.0071). B. Eye size is significantly correlated in the same DGRP strains expressing GMR-rpr and GMR > Rh1 (r = 0.25, P = 0.001). C. Eye size is not correlated in the same DGRP strains expressing GMR > and GMR > Rh1 (r = 0.12, P = 0.13). * P < 0.05, ** P < 0.005.

rpr-induced degeneration is modified by apoptosis, Wnt signaling, and mitochondrial metabolism

Genome-wide association analysis:

To identify the genes driving this phenotypic variability, we performed a genome-wide association analysis to identify genetic polymorphisms that impact the severity of degeneration in the GMR > and GMR-rpr models of apoptosis. We used mean eye size as a quantitative phenotype to test for association with polymorphisms in the DGRP. With the large number of variants (1,967,719 for and 1,962,205 for ) tested for the apoptosis models in ∼200 lines, the analyses were not sufficiently powered for associations to survive multiple testing correction at a false discovery rate of 0.2 (p-value < 1x10−7). Previous work from our lab and others, however, have demonstrated the relevance of weak signals obtained from DGRP studies. Specifically, it has been observed that association of candidate genes can be replicated in different lab environments and in different populations (Everman ; Pitchers ). Ultimately, our goal is not to treat these associations as definitive, but to use these variants to nominate candidate modifier genes or pathways for subsequent functional validation. Despite the weak signals, this approach has been used to successfully identify bona fide modifier genes in previous work (Chow , 2015, 2016; Palu and Chow 2018, 2019; Lavoy ; Ahlers ). Here again, this genome-wide association analysis yields a number of genes that potentially underlie the variable expressivity of degenerative phenotypes induced by these models of apoptosis. Using an arbitrary p-value cutoff of <1x10−04, we identified 128 significantly associated polymorphisms for the GMR-rpr model (Table S3). We only considered polymorphisms that fall within an annotated gene or +/− 1 kb of an annotated gene. Polymorphisms located more than 1 kb from an annotated gene were not considered in our analysis. Sixteen polymorphisms lie outside of these parameters and were not considered further. Of the remaining 112 polymorphisms, ten are located in an intergenic region (+/− 1kb), 14 are located in UTRs, 69 are located in introns, and 19 are located in protein-coding sequences. All 19 polymorphisms in coding regions are synonymous variants. These 112 gene-associated polymorphisms lie in 82 candidate genes (Table S3, S4). Sixty-six of the candidate genes have direct human orthologs (Table S4). A more stringent p-value cutoff (<1x10−5) yields only 20 polymorphisms, 16 of which lie in 14 candidate genes (12 with human orthologs) (Table S3, S4). Because the more stringent cutoff yielded few candidates, we focused the majority of our analysis on the 82 candidate genes identified at P < 1x10−04. A less stringent cutoff allows for enrichment analyses. It will also result in more false positives, but, again, we are mainly focused on the genes, and not the variants, as the variants themselves may not be evolutionarily conserved, but the candidate genes themselves may be conserved as modifier genes. Functional analyses of the candidate genes will point to bona fide modifiers.

For the GMR > model, we identified 24 polymorphisms at a p-value cutoff of <1x10−04 (Table S5). Eight of these polymorphisms lie outside of genes and were not considered further. Of the remaining 16 polymorphisms, one is located in a UTR, 15 are located in introns, and eight are intergenic. The 16 gene-associated polymorphisms lie in 13 candidate genes (Table S5, S6). Thirteen of the associated polymorphisms have a p-value of <1x10−05. Five of these are intergenic, while the remaining six are in six candidate genes. Interestingly, there is no overlap between the GMR > candidate polymorphisms or genes and those identified using the GMR-rpr model of apoptosis (Table S3-S6). The only overlap in modifier genes is between GMR > Rh1 and GMR > (Table S6) (Chow ). They share candidate modifier genes , a disulfide oxidoreductase (FlyBase Curators 2004), and , a cell surface immunoglobulin involved in synapse organization (Gaudet ). It is unclear what the significance of this overlap might be.

We conclude from our initial analysis that the top candidates for our models of degeneration are likely specific to the method by which we induce that degeneration. However, it is also possible that our study is underpowered and we are simply missing overlapping candidate modifiers shared between the two models. Future validation of the candidates will involve specific tests examining the impact of a modifier gene on both apoptotic models.

Even at the lower p-value threshold, there are very few significant associations for GMR > , and even fewer in close proximity to a gene. This is likely because the DGRP population, simply due to chance, does not have a great deal of variation affecting the p53 pathway. Because the DGRP was generated from a single population, the entire spectrum of possible variation is simply not represented. It is plausible that screening through an alternate population might yield different and more interesting results for GMR > modifiers. With the results reported here, it would be difficult to obtain meaningful results from analyses focused on shared gene ontology, modifier interactions, and shared pathway functions. We therefore elected to focus the remaining analysis on the GMR-rpr model.

Modifier genes:

Because the model is the most direct inducer of apoptosis we have tested, we expected to see apoptotic functions for many of the candidate genes identified in our GWAS. The top hit was the gene (), a ubiquitin specific protease (USP) orthologous to human USP53 and USP54 (Table S4). We identified nine intronic SNPs in through our association analysis. Previous studies show that loss of in the developing eye results in a mild rough eye phenotype, albeit a much less dramatic one than that seen upon overexpression of (Wolff and Ready 1991; Copeland ). While this previous study reported no genetic interaction between and , this was assessed based on qualitative changes as opposed to quantitative differences in eye size (Copeland ). Our GWAS data suggests that such a genetic interaction may play an important role in -induced degeneration.

Ec is one of several apoptotic genes identified in this analysis. In fact, 16/82 (∼20%) of the candidate genes have known functions in apoptosis-related pathways, all of which have conserved human orthologs (Table S4). One of these is , a Drosophila paralog of (human orthologs: BIRC2 and BIRC3) (Hay ). The Diap proteins normally inhibit caspase activation and prevent apoptosis. Expression of the rpr/grim/hid proteins inhibits Diap1 and Diap2, allowing apoptosis to proceed. Increased expression or activity of Diap2 reduces the impact of overexpression, thereby reducing apoptosis (Hay ). Conversely, reduced expression of may not have a strong impact on -associated degeneration, as is the major functional paralog in this pathway. The identification of a gene directly involved in the pathway demonstrates the efficacy of our GWAS.

Two candidates, and (ERCC3 and ERCC2) (Table S4), have human orthologs mutated in Xeroderma pigmentosum, an inherited genetic condition where defects in DNA excision repair result in melanomas and eventually death (Kraemer and DiGiovanna 2016). These are subunits of the TFIIH helicase complex that are involved in excision repair after UV damage (Koken ; Mounkes ; Reynaud ). Besides and , we identified 4 additional genes whose human orthologs are directly involved in cancer: (OPCML), Fum4 (FH), (TMEM259), and (BLNK). Mutations in these genes have been associated with ovarian cancer (OPCML) (Sellar ), renal cancer (FH) (The Multiple Leiomyoma Consortium 2002; Pollard ), and various carcinomas (TMEM259) (Chen ). The roles of these genes in cancer are likely due to functions in apoptotic initiation or cell cycle regulation. Other candidates are activated downstream of p53, such as (ADGRB3) and (BAIAP3) (Shiratsuchi , 1998). This suggests that feedback signaling through p53 can increase -induced apoptosis and degeneration.

24/82 candidate genes (∼30%) are involved in neuronal function or implicated in neurological disease. Twenty-three have conserved human orthologs (Table S4). Human orthologs of Form3 (INF2) and (KARS) can both be mutated in different forms of the degenerative peripheral neuropathy Charcot-Marie-Tooth disease (McLaughlin ; Boyer ), while (KCNC3) and (CWF19L1) are associated with spinocerebellar ataxia (Waters ; Burns ). Mutation in the Rab3-interacting scaffold protein encoded by (RIMS1) can cause a retinal degenerative disease that is similar to retinitis pigmentosa (Johnson ), which was the focus of the Rh1 study (Chow ). Identification of genes with roles in different neuronal and muscular degenerative diseases suggests that these modifiers could be important in a variety of apoptosis-associated diseases.

Network analysis:

To understand if there are functional relationships between GMR-rpr modifiers, we examined interactions among the 82 candidate genes. Genetic, physical, and predicted interactions were compiled and visualized using Cytoscape software (Shannon ; Montojo ). Fourteen of the 82 candidate genes were found as nodes in these interaction networks, as was itself (Figure 4A). We identified several interesting clusters of candidate genes, including those with functions in apoptosis, development, and protein ubiquitination.

Figure 4

modifiers are enriched for neuronal function, Wnt signaling, and metabolic pathways. A. modifier network, as plotted by the GeneMANIA plugin in Cytoscape (Shannon ; Montojo ). Significant candidate modifiers are indicated in red, with physical interactions shown in green, genetic interactions shown in blue, and predicted interactions shown in gray. Thicker lines indicate stronger evidence for the association. B. Top 20 significant ontological categories as identified by GSEA. Categories are arranged from most significant on top to least significant along the y-axis. P-values are indicated by red-to-blue gradient, with red the lowest p-values and blue the highest P-values. Gene number identified in each category is indicated along the y-axis.

As expected, given the large number of candidates with apoptotic roles, we found an apoptosis cluster of interactions between modifiers with functions associated with cell cycle regulation and cell death (Figure 4A). A number of these genes, including and (FLNA), have either direct or indirect interactions with itself. As noted above, interacts both physically and genetically with (Figure 1) (Hay ). shows indirect genetic interactions with through its physical association with the presenillin (psn) protein (Guo ) (Figure 4A). This interaction is conserved in humans, and mutations that alter this interaction are associated with Alzheimer’s Disease (Lu ).

Among the apoptosis cluster are also regulators of developmental apoptosis in our network, including the protease (Copeland ) and the neuronal cell adhesion protein encoded by (KIRREL3) (Bao ) (Figure 4A). Indirect genetic interactions were identified between these genes, which are commonly involved in development of the Drosophila eye imaginal disc and accompanying regulated apoptosis. The chromatin-binding HmgD/Z (HMGB2) proteins are expressed at high levels in the larval CNS, suggesting that they are important for the developmental regulation of neuronal gene expression (Churchill ; Gaudet ; Brown ). They indirectly interact with through the closely related dsp1, which encodes another paralog of human HMGB2 (Figure 4A). Dsp1 recruits members of the repressive polycomb complex to chromatin. It is possible that these genetic interactions indicate a role for the HMGB2 proteins in regulating expression and, as a result, developmental regulation of cell death and tissue turnover. Our apoptotic model is expressed in a developmental tissue, suggesting that some of the variation in eye size observed across the DGRP could be due to changes in the response of developmental processes to the abnormal activation of apoptosis. Such regulators of developmental apoptosis could be excellent candidates for therapeutic targeting in degenerative diseases.

We also identified a number of predicted interactions in a cluster of modifier genes involved in protein ubiquitination (Figure 4A). Among the top candidate genes are , , Su(dx) (ITCH), and (RNF7), all of which have important roles in protein degradation through ubiquitination and the proteasome degradation pathway. Su(dx), like , encodes a ubiquitin ligase (Gaudet ). Our network analysis highlights a predicted interaction between Su(dx) and the Rab GTPase-interacting protein , another candidate gene (Laflamme ) (Figure 4A). This regulator of vesicular fusion is predicted to interact with a number of additional ubiquitin ligases as well (Figure 4A). Degradation of proteins through the proteasome is an important mechanism for maintaining cellular homeostasis under a variety of cellular stresses (Sano and Reed 2013). Altered regulation of E3 ligases, which determine the identity and specificity of proteins for degradation (Ester Morreale and Walden 2016), could tip the balance of cells experiencing apoptotic stress toward or away from cell death.

Gene set enrichment analysis:

Thus far, we have focused our analysis on candidate modifiers surviving a specific statistical threshold in our GWAS. While this provides many new avenues for future analysis, it ignores the majority of the association data. We therefore performed gene set enrichment analysis (GSEA), using all GWAS variant data and their associated P-values. The gene nearest to each variant was assigned the variant’s P-value and used as GSEA input, using the method described (Subramanian ). Given a defined set of genes annotated with a certain GO function, GSEA determines whether the members of that set are randomly distributed throughout the ranked list or if they are found primarily at the top or bottom of that list. GO categories enriched at the top of the list functionally describe the phenotype of the gene set. While traditional GO analysis uses a set of genes based on a P-value cutoff, GSEA examines the entire gene set (Dyer ). GSEA identified 62 significantly associated gene sets (≥ 3 genes) at a corrected p-value of <0.05 (Table S7). The top gene set was synaptogenesis (GO:0007416, P = 3.7 × 10−3) and includes (SEMA6A), a conserved semaphorin-binding protein involved in axon guidance (Ayoob ; Gaudet ) and one of the top modifier candidates based on individual polymorphism analysis (Figure 4B, Table S7). Other genes in this category include those involved in synapse formation and organization, suggesting that regulating neuronal connectivity and synapse choice could play a role in the decision to apoptose or to survive.

The second most significantly enriched category was Wnt signaling (GO:0016055, P = 6.7 × 10−3), consisting of 51 enriched genes from our GMR-rpr analysis (Figure 4B, Table S7). One of these, , is also a candidate modifier gene (Table S4,S7). is a Drosophila ortholog of the genes encoding the co-receptors LRP5/6 in canonical Wnt signaling (Rives ). The second most significant single candidate gene in the GWA is swim (TINAGL1/TINAG), a secreted cysteine protease capable of binding the wingless (wg) ligand and enhancing its spread and signaling capabilities (Mulligan ). Also enriched for significant polymorphisms are four frizzled paralogs (Wnt receptors) and six paralogs of the Wnt ligand (Table S7). Other integral components of the canonical Wnt pathway, such as disheveled, axin, and CKIα, are enriched for associated polymorphisms, as are several peripheral and non-canonical regulators of Wnt signaling (Table S7). This striking association is reinforced by previous studies that have linked Wnt signaling with either the promotion or restraint of cell death (Pećina-Slaus 2010). Non-canonical Wnt signaling can activate JNK or calcium release from the ER, both of which can alter the decision to initiate apoptosis (Rasmussen ). It will be interesting to investigate Wnt signaling collectively as well as with individual candidates to determine how different branches of the pathway impact degenerative diseases.

GSEA also identified a number of genes and pathways involved in mitochondrial homeostasis and metabolism (Figure 4B), including malate metabolic processes (seven genes, GO:0006108, P = 0.011). These genes encode for malate dehydrogenase enzymes, six of which are localized to the mitochondrion (Figure 4B, Table S7). Malate dehydrogenase catalyzes the oxidation of malate to oxaloacetate in the last step of the TCA cycle prior to the entrance of acetyl-CoA (Minárik ). The presence of so many paralogs of this enzyme suggests that mitochondrial metabolism, and in particular the mitochondrial redox state, is a major regulator of apoptosis. Supporting this, one of the top candidates, Fum4 (FH), is also an essential enzyme in the TCA cycle (Table S4). The GSEA further supports this finding, as FAD binding is also enriched (48 genes, GO:0050660, P = 0.020) (Figure 4B). A primary function for these 48 enriched genes is the maintenance of redox homeostasis, 16 of which localize to the mitochondria. Another of these genes, the apoptosis-inducing factor , is activated independently from caspases by mitochondrial stress and is released into the cytoplasm, travels to the nucleus, and initiates the chromatin condensation and DNA fragmentation that immediately precedes cell death (Elmore 2007).

More generally, redox homeostasis in other cellular compartments is also implicated by GSEA (Table S7, Figure 4B). Three paralogs of aldehyde oxidase (Aox) and the NAD(P)H oxidoreductase (DUOX1) are enriched for associated polymorphisms; these oxidase enzymes are essential for maintaining an appropriate balance of reactive oxygen species in the cytoplasm. We identified four paralogs of acyl-coA oxidase (Acox), which is involved in the β-oxidation of very long chain fatty acids in the peroxisome, and an additional 4 genes involved in mitochondrial β-oxidation: (ETFA), Mcad (ACADM), (ACADS), and (ACADVL).

The involvement of enzymes regulating redox homeostasis, and more specifically redox homeostasis in the mitochondria, is consistent with -induced apoptosis. Both caspase-dependent and caspase-independent apoptotic pathways can be activated downstream of mitochondrial stress (Elmore 2007; Rasmussen ). Increasing the permeability of the mitochondrial membrane is sufficient to ensure activation of the apoptosome through the release of cytochrome-C (Elmore 2007). This, along with expression of the mitochondria-associated IAP inhibitors rpr/grim/hid, activates the caspase cascade (Sandu ). Damage to the mitochondria that increases permeability, such as through redox stress, is itself sufficient to activate apoptosis in a caspase-independent manner through the release of AIF (Elmore 2007). Importantly, nearly all the mitochondrial genes identified using GSEA analysis are expressed in the imaginal discs (Brown ). This suggests that mitochondrial function modifies apoptotic progression in a cell autonomous fashion, consistent with the known roles of the mitochondria in cell death.

Other metabolic processes such as sterol transport (GO:0015918, P = 0.013), leucine import (GO:0060356, P = 8.9 × 10−3), and fat body development (GO:0007503, P = 0.011) are enriched in the GSEA (Table S7, Figure 4B). Disruption of metabolic processes has long been known to induce oxidative and ER stress, both of which are capable of activating apoptosis through JNK/grm-rpr-hid signaling cascades or directly through mitochondrial stress (Kanda and Miura 2004). It will be interesting to explore how these metabolic processes alter apoptosis not only in this model of retinal degeneration, but in physiologically relevant cell types and tissues, such as the midgut, fat body, and insulin-producing cells.

The enrichment of multiple metabolic categories suggests that the impact of cellular and mitochondrial metabolism on redox homeostasis could play a major role in -induced degeneration. We hypothesize that these regulators of mitochondrial redox state and metabolism are directly and indirectly influencing the activation of mitochondrial proteins involved in the final decision to undergo apoptosis. Our GSEA emphasizes the importance of exploring not just individually associated genes but also their functional pathways and partners when identifying genetic modifiers of disease.

Functional analysis of candidate modifiers of apoptosis

To confirm the roles of our candidate genes in regulating apoptosis, we elected to test the impact of loss of modifier expression for the top nine candidate genes and the top seven candidate genes for which we were able to obtain transgenic RNAi lines. We crossed the RNAi strains targeting each of these modifiers into the GMR-rpr or GMR > line, and then measured the eye area in offspring carrying both the RNAi construct and the apoptosis model (Figure 5, Figure S3). Eye area was quantified and compared to a genetically matched control expressing only the apoptosis model (Figure S4A). We also crossed the RNAi lines to a strain expressing only GMR-GAL4 as a control for independent modifier effect. If loss of the candidate alone has a substantial impact on eye size or phenotype, it would suggest that this candidate is NOT specific to the model in question.

Figure 5

Knockdown of candidate modifiers significantly alters apoptosis-induced degeneration. RNAi against candidate modifiers was expressed under the control of GMR-GAL4 in the GMR-rpr model. The genetically matched attP2 line was crossed into the GMR-rpr line as a control (blue). Eye size in pixels was quantified for N = 11-15 flies per strain and plotted with the 25th-75th percentile of measurements in the central box. Measurements lying outside of 1.5 × interquartile range are indicated as points. Representative images of each line are found above the data for that line. Knockdown of or swim significantly reduces eye size in the GMR-rpr model of degeneration compared to controls. Loss of , , or results in a significant increase in eye size compared to controls. Loss of does not significantly alter eye size, but changes in pigmentation are similar in the presence or absence of GMR-rpr (Fig S4C). Loss of αManA1, , or do not produce a significant effect. RNAi lines with significant changes in eye size are indicated in red, while those that are not significantly changed are indicated in white. * P < 0.05, *** P < 0.0005.

For modifiers of -induced apoptosis, we found that knockdown of Ir40A (15952 ± 716 pixels, N = 15) expression resulted in enhancement of the apoptosis phenotype, showing a significant decrease in eye size compared to controls expressing only GMR > (17459 ± 492 pixels, N = 9) (Figure S3A). Knockdown of shab (KCNB1) (16598 ± 1143 pixels, N = 10), (CCNE1) (17462 ± 581 pixels, N = 15), or (CELG2) (16598 ± 1143 pixels, N = 10), resulted in a partial rescue, with a significant increase in eye size compared to controls expressing GMR > (18022 ± 884 pixels, N = 15) (Figure S3B). However, knockdown of either or in the GMR-GAL4 (no apoptosis control) background was also sufficient to increase eye size (Figure S3C), suggesting that these modifiers might be working independently from overexpression. No significant change in eye size was observed upon knockdown of cyt-c-P (CYCS) (16836 ± 1052 pixels, N = 15), (GRXCR1) (17102 ± 934 pixels, N = 15), or (TIAM1) (17218 ± 939 pixels, N = 15) as compared to controls expressing only GMR > (17459 ± 492 pixels, N = 9) (Figure S3A). While some of these genes may still be interesting modifiers that would benefit from follow-up studies, our results suggest that the modifiers identified in the GMR > screen are either of weak individual effect or, potentially, false positives. This also demonstrates the importance of validation studies such as these when analyzing GWA candidates.

We next focused our analysis on the modifiers. Knockdown of either () (16183 ± 875 pixels, N = 15) or swim expression (15518 ± 2418 pixels, N = 14) resulted in enhancement of the apoptosis phenotype, showing a significant decrease in eye size compared to controls expressing only GMR-rpr (17534 ± 1098 pixels, N = 11) (Figure 5). Knockdown of sema1a (18990 ± 746 pixels, N = 15), () (20323 ± 622 pixels, N = 15), or (20240 ± 617 pixels, N = 14) resulted in a partial rescue, with a significant increase in eye size compared to controls expressing GMR-rpr (Figure 5). No significant change in eye size was observed upon knockdown of (GZF1) (18525 ± 449 pixels, N = 12), (KARS) (17879 ± 1834 pixels, N = 12), αMan-1A (MAN1A2) (17842 ± 763 pixels, N = 15), or (SLC25A11) (18755 ± 787 pixels, N = 13) (Figure 5). In the absence of overexpression, RNAi revealed no significant change in phenotype from controls (Figure S4B,C). These results suggest that many of the top GWA candidate modifiers are capable of modifying the apoptotic phenotypes associated with the GMR-rpr model of degeneration.

It is of course theoretically possible that we may have observed a positive result for any five of nine randomly selected genes within the genome due to the fact that the genetic background is sensitized by the over-expression of rpr. Our results with the similar GMR > model suggest that this is not likely to be the case. We examined seven candidate modifiers of similar significance as the GMR-rpr candidates and observed only two that specifically affect eye size when apoptosis is activated. This suggests that while, as previously postulated, the GMR > candidates could be false positives, the GMR-rpr candidates are more likely to be true positives. Other groups using the DGRP have found this to be true as well (Vonesch ).

In the future we will also examine the impact of overexpression of candidate genes on the GMR-rpr model of apoptosis, as some candidate genes may exert a stronger influence under conditions of gain of function, rather than loss of function.

Conclusions

The primary goal of this study was to identify candidate genes and pathways that modify apoptosis and degenerative processes. Apoptosis is a primary cause of disease in a multitude of degenerative disorders (Mattson 2000). It is also a commonly targeted pathway for cancer therapies (Ouyang ). These and other diseases are subject to a large degree of phenotypic heterogeneity due to inter-individual differences in genetic background among patients (Queitsch ; Chow 2016). Our results point to stress-associated signals as important modifiers of apoptosis-induced degeneration. Overexpression of the apoptotic gene can be induced transcriptionally by either p53 or by stress-associated transcription factors such as Jun (Shlevkov and Morata 2012). Many of the candidate modifiers of the GMR-rpr model feed into the stress response through mitochondrial metabolism, Wnt signaling, and protein degradation. Their identification as modifier pathways indicates that activity is being modulated by feedback signals through these pathways and suggests that apoptotic diseases might be targeted through these stress-related pathways. Understanding how genetic diversity in the population impacts apoptosis could therefore lead to identification of prognostic predictors in the diagnosis of disease and of new therapeutic targets. The modifiers identified here inform our understanding of cell death regulation and could serve as therapeutic targets in a variety of apoptosis-related disorders.