The Gene scb-1 Underlies Variation in Caenorhabditis elegans Chemotherapeutic Responses.

Pleiotropy, the concept that a single gene controls multiple distinct traits, is prevalent in most organisms and has broad implications for medicine and agriculture. The identification of the molecular mechanisms underlying pleiotropy has the power to reveal previously unknown biological connections between seemingly unrelated traits. Additionally, the discovery of pleiotropic genes increases our understanding of both genetic and phenotypic complexity by characterizing novel gene functions. Quantitative trait locus (QTL) mapping has been used to identify several pleiotropic regions in many organisms. However, gene knockout studies are needed to eliminate the possibility of tightly linked, non-pleiotropic loci. Here, we use a panel of 296 recombinant inbred advanced intercross lines of Caenorhabditis elegans and a high-throughput fitness assay to identify a single large-effect QTL on the center of chromosome V associated with variation in responses to eight chemotherapeutics. We validate this QTL with near-isogenic lines and pair genome-wide gene expression data with drug response traits to perform mediation analysis, leading to the identification of a pleiotropic candidate gene, scb-1, for some of the eight chemotherapeutics. Using deletion strains created by genome editing, we show that scb-1, which was previously implicated in response to bleomycin, also underlies responses to other double-strand DNA break-inducing chemotherapeutics. This finding provides new evidence for the role of scb-1 in the nematode drug response and highlights the power of mediation analysis to identify causal genes.

Pleiotropy refers to the well established notion that a single gene or genetic variant affects multiple distinct traits (Paaby and Rockman 2013), and the discovery of pleiotropic genes can provide meaningful insights into the molecular mechanisms of these traits (Tyler ). It has become easier to identify pleiotropic genes with the advent of reverse-genetic screens and quantitative trait locus (QTL) mapping (Paaby and Rockman 2013). For example, pleiotropic QTL for diverse growth and fitness traits have been identified in organisms such as yeast (Cubillos ; Jerison ; Peltier ), Arabidopsis (McKay ; El‐Assal ; Fusari ), Drosophila (Brown ; McGuigan ), and mice (White ; Leamy ; Lin ). These studies have led to important questions in the field of evolutionary genetics regarding the ‘cost of complexity’ (Fisher; Orr 2000), as a single mutation might be beneficial for one trait and harmful for another (Wagner and Zhang 2011). Furthermore, human association studies have identified pleiotropic variants associated with different diseases (Sivakumaran ; Pavlides ; Chesmore ), highlighting both the ubiquity and importance of certain immune-related genes and oncogenes across unrelated diseases (Borrello ; Gratten and Visscher 2016). Perhaps the strongest evidence of pleiotropy exists for molecular phenotypes. Large-scale expression QTL (eQTL) mapping studies have identified single regulatory variants that control expression and likely the functions of hundreds of genes at once, opening a window into the mechanisms for how traits are controlled (Keurentjes ; Breitling ; Rockman ; Albert and Kruglyak 2015; Hasin-Brumshtein ; Albert ).

The nematode Caenorhabditis elegans provides a tractable metazoan model to identify and study pleiotropic QTL (Paaby and Rockman 2013). A large panel of recombinant inbred advanced intercross lines (RIAILs) derived from two divergent strains, N2 and CB4856 (Rockman and Kruglyak 2009; Andersen ), has been leveraged in several linkage mapping analyses (Li ; Gutteling , 2007a; Kammenga ; Seidel , 2011; Reddy ; McGrath ; Doroszuk ; Viñuela ; Rockman ; Bendesky , 2012; Bendesky and Bargmann 2011; Rodriguez ; Andersen ; Glater ; Snoek ; Balla ; Schmid ; Singh ; Zdraljevic , 2019; Lee ; Zamanian ; Evans ; Brady ). Quantitative genetic analysis using these panels and a high-throughput phenotyping assay (Andersen ) has facilitated the discovery of numerous QTL (Zamanian ), several quantitative trait genes (QTG) (Brady ) and quantitative trait nucleotides (QTN) (Zdraljevic , 2019) underlying fitness-related traits in the nematode. Additionally, three pleiotropic genomic regions were recently found to influence responses to a diverse group of toxins (Evans ). However, overlapping genomic regions might not represent true pleiotropy but could demonstrate the co-existence of tightly linked loci (Paaby and Rockman 2013).

Here, we use linkage mapping to identify a single overlapping QTL on chromosome V that influences the responses to eight chemotherapeutic compounds. We show that these drug-response QTL also overlap with an expression QTL hotspot that contains the gene , previously implicated in bleomycin response (Brady ). Although the exact mechanism of is yet unknown, it is hypothesized to act in response to stress (Riedel ) and has weak homology to a viral hydrolase (Kelley ; Zhang ). Together, these data suggest that the importance of expression might extend beyond bleomycin response. We validated the QTL using near-isogenic lines (NILs) and performed mediation analysis to predict that expression explains the observed QTL for four of the eight drugs. Finally, we directly tested the effect of loss of function on chemotherapeutic responses. We discovered that expression of underlies differential responses to several chemotherapeutics that cause double-strand DNA breaks, not just bleomycin. This discovery of pleiotropy helps to further define the role of by expanding its known functions and provides insights into the molecular mechanisms underlying the nematode drug response.

Materials and Methods

Strains

Animals were grown at 20° on modified nematode growth media (NGMA) containing 1% agar and 0.7% agarose to prevent burrowing and fed OP50 (Ghosh ). The two parental strains, the canonical laboratory strain, N2, and the wild isolate from Hawaii, CB4856, were used to generate all recombinant lines. 208 recombinant inbred advanced intercross lines (RIAILs) generated previously by Rockman et al. (Rockman and Kruglyak 2009) (set 1 RIAILs) were phenotyped for expression QTL mapping (detailed below). A second set of 296 RIAILs generated previously by Andersen et al. (Andersen ) (set 2 RIAILs) was used more extensively for drug phenotyping and linkage mapping. The set 2 RIAILs were used for linkage mapping because they addressed the three main disadvantages of the set 1 RIAILs detailed previously (Andersen ), namely a structured population, the laboratory-derived variant in (Sterken ), and the incompatibility (Seidel , 2011). Because of these limitations, the set 2 RIAILs were generated using QX1430 and CB4856. QX1430 is from the N2 strain background but contains a transposon insertion in and the CB4856 allele introgressed on chromosome X (Andersen ). Near-isogenic lines (NILs) were generated by backcrossing a selected RIAIL for several generations to the parent strain (N2 or CB4856) (Brady ) using PCR amplicons for insertion-deletion (indels) variants to track the introgressed region. NILs were whole-genome sequenced to verify introgressions were only in the targeted genomic intervals. CRISPR-Cas9-mediated deletions of were described previously (Brady ). All strains are available upon request or from the C. elegans Natural Diversity Resource (Cook et al. 2017). Primers used to generate ECA1114 can be found in the Supplemental Information.

High-throughput fitness assays for linkage mapping

For dose responses and RIAIL phenotyping, we used a high-throughput fitness assay described previously (Andersen ). In summary, populations of each strain were passaged and amplified on NGMA plates for four generations. In the fifth generation, gravid adults were bleach-synchronized and 25-50 embryos from each strain were aliquoted into 96-well microtiter plates at a final volume of 50 µL K medium (Boyd ). The following day, arrested L1s were fed HB101 bacterial lysate (Pennsylvania State University Shared Fermentation Facility, State College, PA; (García-González )) at a final concentration of 5 mg/mL in K medium and were grown to the L4 larval stage for 48 hr at 20° with constant shaking. Three L4 larvae were sorted into new 96-well microtiter plates containing 10 mg/mL HB101 bacterial lysate, 50 µM kanamycin, and either diluent (1% water or 1% DMSO) or drug dissolved in the diluent using a large-particle flow cytometer (COPAS BIOSORT, Union Biometrica; Holliston, MA). Sorted animals were grown for 96 hr at 20° with constant shaking. The next generation of animals and the parents were treated with sodium azide (50 mM in 1X M9) to straighten their bodies for more accurate length measurements. Animal length (median.TOF), optical density integrated over animal length (median.EXT), and brood size (norm.n) were quantified for each well using the COPAS BIOSORT. Nematodes get longer (animal length) and become thicker and more complex (optical density) over developmental time. Because animal length and optical density are highly correlated, we calculated a fourth trait (median.norm.EXT) that normalizes optical density by animal length (median.EXT / median.TOF). Phenotypic measurements collected by the BIOSORT were processed and analyzed using the R package easysorter (Shimko and Andersen 2014) as described previously (Brady ). Differences among strains within the control conditions were controlled by subtracting the mean control-condition value from each drug-condition replicate for each strain using a linear model (drug_phenotype ∼ mean_control_phenotype). In this way, we are addressing only the differences among strains that were caused by the drug condition and the variance in the control condition does not affect the variance in the drug condition. For plotting purposes, these residual values (negative and positive residuals) were normalized from 0 to 1 where 0 refers to the smallest residual phenotypic value in that condition and 1 refers to the largest.

Dose-response assays

Four genetically divergent strains (N2, CB4856, JU258, and DL238) were treated with increasing concentrations of each of the eight drugs using the high-throughput fitness assay described above. The dose of each drug that provided a reproducible drug-specific effect that maximizes between-strain variation while minimizing within-strain variation across the four traits was selected for the linkage mapping experiments. The chosen concentrations are as follows: 100 µM amsacrine hydrochloride (Fisher Scientific, #A277720MG) in DMSO, 50 µM bleomycin sulfate (Fisher, #50-148-546) in water, 2 µM bortezomib (VWR, #AAJ60378-MA) in DMSO, 250 µM carmustine (Sigma, #1096724-75MG) in DMSO, 500 µM cisplatin (Sigma, #479306-1G) in K media, 500 µM etoposide (Sigma, #E1383) in DMSO, 500 µM puromycin dihydrochloride (VWR, #62111-170) in water, and 150 µM silver nitrate (Sigma-Aldrich, #S6506-5G) in water.

Linkage mapping

Set 1 and set 2 RIAILs were phenotyped in each of the eight drugs and controls using the high-throughput fitness assay described above. Linkage mapping was performed on each of the drug and gene expression traits using the R package linkagemapping (https://github.com/AndersenLab/linkagemapping) as described previously (Brady ). The cross object derived from the whole-genome sequencing of the RIAILs containing 13,003 SNPs was loaded using the function load_cross_obj(“N2xCB4856cross_full”). The RIAIL phenotypes were merged into the cross object using the merge_pheno function with the argument set = 1 for expression QTL mapping and set = 2 for drug phenotype mapping. A forward search (fsearch function) adapted from the R/qtl package (Broman ) was used to calculate the logarithm of the odds (LOD) scores for each genetic marker and each trait as -n(ln(1-R where R is the Pearson correlation coefficient between the RIAIL genotypes at the marker and trait phenotypes (Bloom ). A 5% genome-wide error rate was calculated by permuting the RIAIL phenotypes 1000 times. The marker with the highest LOD score above the significance threshold was selected as the QTL then integrated into the model as a cofactor and mapping was repeated iteratively until no further significant QTL were identified. Finally, the annotate_lods function was used to calculate the effect size of each QTL and determine 95% confidence intervals defined by a 1.5 LOD drop from the peak marker using the argument cutoff = proximal.

Modified high-throughput fitness assay for NIL validation

NILs and deletion strains were tested using a modified version of the high-throughput fitness assay detailed above. Strains were propagated for two generations, bleach-synchronized in three independent replicates, and titered at a concentration of 25-50 embryos per well of a 96-well microtiter plate. The following day, arrested L1s were fed HB101 bacterial lysate at a final concentration of 5 mg/mL with either diluent or drug. After 48 hr of growth at 20° with constant shaking, nematodes were treated with sodium azide (5 mM in water) prior to analysis of animal length and optical density using the COPAS BIOSORT. As only one generation of growth is observed, brood size was not calculated. A single trait (median.EXT) was chosen to represent animal growth generally, as the trait is defined by integrating optical density over length. Because of the modified timing of the drug delivery, lower drug concentrations were needed to recapitulate the previously observed phenotypic effect. The selected doses are as follows: 12.5 µM amsacrine in DMSO, 12.5 µM bleomycin in water, 2 µM bortezomib in DMSO, 250 µM carmustine in DMSO, 125 µM cisplatin in K media, 62.5 µM etoposide in DMSO, 300 µM puromycin in water, and 100 µM silver in water.

Expression QTL analysis

Microarray data for gene expression using 15,888 probes were previously collected from synchronized young adult populations of 208 set 1 RIAILs (Rockman ). Expression data were corrected for dye effects and probes with variants were removed (Andersen ). Linkage mapping was performed as described above for the remaining 14,107 probes, and a significance threshold was determined using a permutation-based False Discovery Rate (FDR). FDR was calculated as the ratio of the average number of genes across 10 permutations expected by chance to show a maximum LOD score greater than a particular threshold vs. the number of genes observed in the real data with a maximum LOD score greater than that threshold. We calculated the FDR for a range of thresholds from 2 to 10, with increasing steps of 0.01, and set the threshold so that the calculated FDR was less than 5%.

Local eQTL were defined as linkages whose peak LOD scores were within 1 Mb of the starting position of the probe (Rockman ). eQTL hotspots were identified by dividing the genome into 5 cM bins and counting the number of distant eQTL that mapped to each bin. Significance was determined as bins with more eQTL than the Bonferroni-corrected 99th percentile of a Poisson distribution with a mean of 3.91 QTL (total QTL / total bins) (Brem ; Rockman ; Evans ). We identified nine eQTL hotspots (II, IVL, IVC, IVR, VL, VC, VR, XL, and XC). To avoid false positives, we increased the LOD threshold for QTL to be counted in the hotspot analysis to a LOD > 5 or LOD > 6. At a LOD > 5, six of the nine eQTL hotspots persist (IVL, IVR, VC, VR, XL, and XC), and at a LOD > 6, three persist (IVL, IVR, and XL). We further looked for spurious eQTL hotspots in ten permuted datasets. At a LOD > 5, we identified four hotspots, and at a LOD > 6, we identified one hotspot.

Mediation analysis

A total of 159 set 1 RIAILs were phenotyped in each of the eight drugs and controls using the standard high-throughput fitness assay described above. Mediation scores were calculated with bootstrapping using the mediate function from the mediation R package (version 4.4.7) (Tingley ) for each QTL identified from the set 1 RIAILs and all 49 probes (including , A_12_P104350) that mapped to the chromosome V eQTL hotspot using the following models:The output of the mediate function can be summarized as follows: the total effect of genotype on phenotype, ignoring expression (tau.coef); the direct effect of genotype on phenotype, while holding expression constant (z0); the estimated effect of expression on phenotype (d0); the proportion of the total effect that can be explained by expression data (n0). This mediation proportion (n0) can be a useful way to identify the impact of gene expression on the overall phenotype. However, cases of inconsistent mediation (where the direct effect is either smaller than or in the opposite direction of the indirect mediation effect) render this measurement uninterpretable with values greater than one or less than zero (MacKinnon ). We used the estimated effect of expression on phenotype (z0) as the final mediation score for this reason. Because the effect size can be positive or negative, mediation scores range from -1 to 1, and we evaluated the absolute value of mediation estimates to compare across traits. Each mediation estimate generated a p-value, indicating confidence in the estimate, derived from bootstrapping with 1000 simulations. The likelihood of mediating a given QTL effect was calculated relative to the other 48 probes with an eQTL in the region (Table S1). Traits in which was at or above the 90th percentile of this distribution were prioritized over other traits.

Statistical analysis

Broad-sense heritability was calculated from the dose response phenotypes using the lmer function in the lme4 R package (Bates ) with the formula phenotype ∼1 + (1|strain) for each dose. For the NIL and deletion high-throughput assays, statistical significance of phenotypic differences between each strain pair was tested using the TukeyHSD function (R Core Team 2017) on an ANOVA model with the formula phenotype ∼ strain to assess differences between strains in the control-regressed phenotype data.

Data availability

File S1 contains the results of the original dose response high-throughput fitness assay. File S2 contains the residual phenotypic values for all 159 set 1 RIAILs, 296 set 2 RIAILs, and parent strains (N2 and CB4856) in response to all eight chemotherapeutics. File S3 contains the linkage mapping results for the set 2 RIAILs for all 32 drug-response traits tested in the high-throughput fitness assay. File S4 is a VCF that reports the genotype of ECA1114. File S5 contains the simplified genotype of all the NILs in the study. File S6 contains the raw pruned phenotypes for the NIL dose response with the modified high-throughput fitness assay. File S7 contains the pairwise statistical significance for all strains and high-throughput assays. File S8 contains the microarray expression data for 14,107 probes from Rockman et al. 2010. File S9 contains the linkage mapping results for the expression data obtained with the set 1 RIAILs. File S10 contains the location of each eQTL hotspot and a list of genes with an eQTL in each hotspot. File S11 contains the linkage mapping results from the set 1 RIAILs for all 32 drug-response traits tested in the high-throughput fitness assay. File S12 contains the pairwise mediation estimates for all 32 drug-response traits and all 49 probes. File S13 contains the raw pruned phenotypes for the deletion modified high-throughput fitness assay. The datasets and code for generating figures can be found at https://github.com/AndersenLab/scb1_mediation_manuscript. Supplemental material available at figshare: https://doi.org/10.25387/g3.12250091.

Results

Natural variation on chromosome V underlies differences in responses to several chemotherapeutics

We measured C. elegans development and chemotherapeutic sensitivity as a function of animal length (median.TOF), optical density (median.EXT), and brood size (norm.n) with a high-throughput assay developed using the COPAS BIOSORT (see Methods) (Andersen ; Zdraljevic , 2019; Evans ; Brady ). Animal length and optical density (animal thickness and composition) are both measures of nematode development, and brood size is a measure of nematode reproduction (Andersen ). Because optical density is calculated as a function of length and these traits are related, a fourth trait that captures the optical density normalized by length (median.norm.EXT) was also included. We exposed four genetically divergent strains (N2, CB4856, JU258, and DL238) to increasing doses of eight chemotherapeutic compounds. Five of these compounds (bleomycin, carmustine, etoposide, amsacrine, and cisplatin) are known to cause double-strand DNA breaks and/or inhibit DNA synthesis (Dorr 1992; Ketron ; Dasari and Tchounwou 2014; Montecucco ; Nikolova ). The remaining three compounds either inhibit protein synthesis (puromycin) (Azzam and Algranati 1973), inhibit the proteosome and subsequent protein degradation (bortezomib) (Piperdi ), or cause cellular toxicity in a poorly defined way (silver nitrate) (Kaplan ) (Table 1). In the presence of each drug, nematodes were generally shorter, less optically dense, and produced smaller broods compared to non-treated nematodes (Figure S1, File S1). We observed significant phenotypic variation among strains and identified a substantial heritable genetic component for most traits (average H2 = 0.52 +/− 0.53).

Table 1

Main mechanism of action for eight chemotherapeutic drugs

Drug	Drug class	Mechanism of action
Amsacrine	Topoisomerase inhibitors	DNA intercalation and inhibition of topoisomerase II, causing DNA double-strand breaks, cell cycle arrest, and cell death
Bleomycin	Antitumor antibiotic	Forms complexes with iron that reduce molecular oxygen to form free radicals which in turn cause DNA single- and double-strand breaks
Bortezomib	Proteosome inhibitors	Reversibly inhibits the 26S proteosome and inhibits nuclear factor (NF)-kappaB causing disruption of various cell signaling pathways, cell cycle arrest, and cell death.
Carmustine	Alkylating agents	Alkylates and cross-links DNA causing cell cycle arrest and cell death
Cisplatin	Alkylating agents	Alkylates and cross-links DNA causing cell cycle arrest and cell death
Etoposide	Topoisomerase inhibitors	Binds to and inhibits topoisomerase II causing an increase of DNA single- and double-strand breaks, cell cycle arrest, and cell death
Puromycin	Aminonucleoside antibiotic	Acts as analog of 3′ terminal end of aminoacyl-tRNA and incorporates itself into growing polypeptide chain causing premature termination and inhibition of protein synthesis
Silver	NA	Multi-faceted induction of apoptosis

We exposed a panel of 296 RIAILs (set 2 RIAILs, see Methods) to all eight chemotherapeutics at a selected concentration that both maximizes among-strain and minimizes within-strain phenotypic variation (File S2). Linkage mapping for all four traits for each of the eight drugs (total of 32 traits) identified 79 QTL from 31 traits (one trait had no significant QTL), several of which have been identified previously (Zdraljevic ; Evans ; Brady ) (File S3, Figure S2). Strikingly, a QTL on the center of chromosome V was linked to variation in responses to all eight compounds (Figure 1). In all cases, the CB4856 allele on chromosome V is associated with greater resistance to the drug than the N2 allele (Figure S2, File S2, File S3). We previously identified this genomic interval as a QTL hotspot, defined as a region heavily enriched for toxin-response QTL (Evans ). Because several of the chemotherapeutics share a similar mechanism of action, a single pleiotropic gene might underlie the observed QTL for multiple drugs.

Figure 1

A large-effect QTL on the center of chromosome V underlies responses to several chemotherapeutics. Linkage mapping results with the set 2 RIAILs for a representative trait for each drug are shown (amascrine: median.norm.EXT, bleomycin median.TOF:, bortezomib: median.TOF, carmustine: norm.n, cisplatin: median.TOF, etoposide: norm.n, puromycin: median.TOF, silver: median.TOF). Genomic position (x-axis) is plotted against the logarithm of the odds (LOD) score (y-axis) for 13,003 genomic markers. Each significant QTL is indicated by a red triangle at the peak marker, and a blue rectangle shows the 95% confidence interval around the peak marker. The percentage of the total variance in the RIAIL population that can be explained by each QTL is shown above the QTL. The dotted vertical line represents the genomic position of .

In order to isolate and validate the effect of this QTL, we constructed reciprocal near-isogenic lines (NILs) by introgressing a genomic region on chromosome V from the resistant CB4856 strain into the sensitive N2 background and vice versa (File S4, File S5). We used a modified high-throughput assay (see Methods) to measure length and optical density of a population of animals grown in the presence of the drug for 48 hr (from larval stages L1 to L4). In this modified assay, less drug was required to observe the same phenotypic effect as before (Figure S3, File S6). Statistical significance was calculated in a pairwise manner for each strain (see Methods; File S7). For all eight chemotherapeutics tested, the strain with the N2 introgression was significantly more sensitive than its CB4856 parent and/or the strain with the CB4856 introgression was significantly more resistant than its N2 parent (Figure 2, File S6, File S7). These data confirm that one or more genetic variant(s) within this region on chromosome V cause increased drug sensitivities in N2.

Figure 2

Near-isogenic lines validate the chromosome V QTL. (A) NIL genotypes on chromosome V are shown, colored orange (N2) and blue (CB4856). From top to bottom, strains are N2, ECA232, ECA1114, and CB4856. The dotted vertical line represents the location of . (B) NIL phenotypes in eight chemotherapeutics (12.5 µM amsacrine, 12.5 µM bleomycin, 2 µM bortezomib, 250 µM carmustine, 125 µM cisplatin, 62.5 µM etoposide, 300 µM puromycin, and 100 µM silver) are plotted as Tukey box plots with strain (y-axis) by relative median optical density (median.EXT, x-axis). Statistical significance was calculated for each strain pair (File S7). Significance of each strain compared to its parental strain (ECA232 to N2 and ECA1114 to CB4856) is shown above each strain pair and colored by the parent strain against which it was tested (ns = non-significant (p-value > 0.05); *, **, ***, and **** = significant (p-value < 0.05, 0.01, 0.001, or 0.0001, respectively).

Expression QTL mapping identifies a hotspot on the center of chromosome V

Genetic variation can affect a phenotype most commonly through either modifications of the amino acid sequence that lead to altered protein function (or even loss of function) or changes in the expression level of the protein. In the latter case, measuring the intermediate phenotype (gene expression) can be useful in elucidating the mechanism by which genetic variation causes phenotypic variation. More specifically, cases with overlap between expression QTL (eQTL) and drug-response QTL suggest that a common variant could underlie both traits and provide evidence in support of causality for the candidate gene in question (Huang ; Sasaki ).

To identify such cases of overlap between expression QTL and the drug-response QTL on chromosome V, we need genome-wide expression data for the RIAILs. In a previous study, expression of 15,888 probes were measured using microarrays for a panel of 208 RIAILs (set 1 RIAILs, see Methods) between N2 and CB4856 (Rockman and Kruglyak 2009) (File S8). This study used the variation in gene expression as a phenotypic trait to identify eQTL using linkage mapping with 1,455 variants (Rockman ). They identified 2,309 eQTL and three regions with significantly clustered distant eQTL (eQTL hotspots), suggesting that these regions are pleiotropic, wherein one or more variant(s) are affecting expression of multiple genes. We recently performed whole-genome sequencing for these strains and identified 13,003 informative variants (Brady ). Using this new set of variants, we re-analyzed the eQTL mapping by performing linkage mapping analysis for a selected 14,107 of the 15,888 probes without genetic variation in CB4856 (Andersen ). We identified 2,540 eQTL associated with variation in expression of 2,196 genes (Figure 3A, File S9). These eQTL have relatively large effect sizes compared to the drug-response QTL. On average, each eQTL explains 23% of the phenotypic variance in gene expression among the RIAIL population. Half of the eQTL (50.2%; 1,276) mapped within 1 Mb of the gene whose expression was measured and were classified as local (see Methods) (Albert and Kruglyak 2015). The other half (49.7%; 1,264) were found distant from their respective gene, and over a third (37%; 940) were found on different chromosomes entirely. In general, eQTL effect sizes increased, max LOD scores decreased, and confidence intervals became smaller compared to the original mapping results (File S9). These differences and the additional eQTL observed between this analysis and the original are possibly caused by the integration of new genetic markers. Additionally, we found several differences in methodology between the current approach and the previous one. These differences include ignoring the population structure of the set 1 RIAILs, adding the forward-search marker-regression linkage mapping, and altering the linkage mapping method itself (see Methods, (Rockman )).

Figure 3

Expression QTL mapping identifies several hotspots. (A) The genomic locations of the eQTL peaks derived from linkage mapping using the set 1 RIAILs (x-axis) are plotted against the genomic locations of the probe (y-axis). The size of the point corresponds to the effect size of the QTL. eQTL are colored by the LOD score, increasing from purple to pink to yellow. The diagonal band represents local eQTL, and vertical bands represent eQTL hotspots. (B) Quantification of eQTL hotspots identified by overlapping distant eQTL. The number of distant eQTL (y-axis) in each 5 cM bin across the genome (x-axis) is shown. Bins above the red line are significant and marked with an asterisk. The bins with the blue asterisks are most significant and have been identified in a previous analysis. The dotted vertical line represents the genomic position of . Gray rectangles below the plot represent locations of the drug-response QTL hotspots previously identified.

We noticed regions of the genome that appeared to be enriched for distant eQTL. We identified eQTL hotspots in a similar manner to the previous study (see Methods) and found a total of nine eQTL hotspots (Figure 3B, File S10). Six of the nine eQTL hotspots withstood more stringent filtering methods (see Methods), and three (left of chromosome IV, right of chromosome IV, and left of chromosome X) were the most significant. These three hotspots also overlap with the most significant eQTL hotspots in the previous study (Rockman ). Notably, three of the eQTL hotspots (center of chromosome IV, right of chromosome IV, and center of chromosome V) overlap with the previously identified drug-response QTL hotspots on chromosomes IV and V (Figure 3B) (Evans ). The overlap of these eQTL and drug-response QTL hotspots could provide strong candidate genes whose expression underlies the differences in nematode drug responses generally. Expression of one gene of interest, , has been previously implicated in response to bleomycin (Brady ) and resides within the eQTL hotspot region on the center of chromosome V (File S10, Table S1). Although the exact mechanism of how responds to bleomycin is unknown, its putative hydrolase activity (Kelley ; Zhang ; Brady ) suggests that it might act to break down chemotherapeutic compounds. These data suggest that variation in expression of and responses to these eight chemotherapeutics (including bleomycin) could be mechanistically linked through the metabolic breakdown of chemotherapeutic drugs.

Mediation analysis suggests that scb-1 expression plays a role in responses to several chemotherapeutics

Mediation analysis seeks to explain the relationship between an independent and a dependent variable by including a third intermediary variable. We can use mediation analysis to understand how certain genetic variants on chromosome V (independent variable) affect drug responses (dependent variable) through differential gene expression of genes within the eQTL hotspot (mediator variable) (Figure S4). We measured brood size, animal length, and optical density in response to all eight chemotherapeutics in the set 1 RIAILs and performed linkage mapping for these traits (File S2, File S11, Figure S5). Although the power to detect QTL with these strains is lower than in our original mapping set (set 2 RIAILs; see Methods) (Andersen ), we still identified overlapping QTL on chromosome V for half of the drugs tested (bleomycin, cisplatin, silver, and amsacrine) (Figure S5, File S11).

We calculated the effect that variation in expression of had on drug-response traits compared to the other 48 genes with an eQTL in the chromosome V eQTL hotspot using mediation analysis (see Methods). We estimated that the effect of expression variation of on bleomycin response is 0.65 (set 1 RIAILs, Figure 4, Figure S6, File S12). Moreover, out of all 49 genes with an eQTL in the region (Table S1), was a clear mediation score outlier. All of the remaining three chemotherapeutics with a QTL on the center of chromosome V in the set 1 RIAIL mapping showed moderate evidence of mediation, with falling well above the 90th percentile of mediation estimates for all genes with an eQTL in this region (Figure 4, Figure S6, File S12). We further performed this mediation analysis on all 32 drug-response traits, regardless of the presence of a QTL in the set 1 RIAIL panel (Figure S6, File S12). Etoposide and puromycin also showed evidence of mediation. This in silico approach indicated that expression of might be an intermediate link between genetic variation on chromosome V and responses to several of the chemotherapeutic drugs tested.

Figure 4

Mediation analysis for the eQTL hotspot on the center of chromosome V. Mediation estimates calculated as the indirect effect that differences in expression of each gene plays in the overall phenotype (y-axis) are plotted against the genomic position of the eQTL (x-axis) on chromosome V for 49 probes (including (red diamond)) that map to the chromosome V eQTL hotspot (set 1 RIAILs). A representative trait for each drug from the set 1 linkage mapping analysis are shown: amascrine (median.EXT), bleomycin (median.EXT), cisplatin (median.TOF), and silver (median.norm.EXT). The 90th percentile of the distribution of mediation estimates for each trait are represented by the horizontal gray lines. The confidence intervals for the QTL (set 1 RIAILs) are shown with the vertical blue dotted lines. The confidence of the estimate increases (p-value decreases) as points become less transparent.

Expression of scb-1 affects responses to several chemotherapeutics that cause double-strand DNA breaks

To empirically test whether expression modulates the chromosome V QTL effect for each drug, we used the modified high-throughput assay (see Methods) to expose two independently derived strains with deletions (Brady ) to each chemotherapeutic (Figure 5, Figure S7, File S13). Statistical significance was calculated in a pairwise manner for each strain (see Methods; File S7). Because RIAILs with the CB4856 allele on chromosome V express higher levels of than RIAILs with the N2 allele (File S8, File S9), we expect that loss of will cause increased drug sensitivity in the CB4856 background but might not have an effect in the N2 background. We validated the results of Brady et al. and confirmed that ablated expression causes hyper-sensitivity to bleomycin in both N2 and CB4856 (Figure 5, Figure S7, Figure S8 File S7, File S13). We also observed similarly increased sensitivity to cisplatin with deletions in both backgrounds. Furthermore, removing shows moderately increased sensitivity in the CB4856 background for amsacrine and in the N2 background for carmustine. The remaining four drugs did not show a significantly different phenotype between the parental N2 and CB4856 strains, suggesting these traits might be less reproducible or that variation does not underlie these drug differences. Overall, these results provide evidence for the pleiotropic effect of , which appears to mediate responses to at least four of the eight chemotherapeutic drugs.

Figure 5

Testing the role of in drug responses. (A) Strain genotypes on chromosome V are shown, colored orange (N2) and blue (CB4856). From top to bottom, strains are N2, ECA1132, ECA1134, and CB4856. Deletion of is indicated by a gray triangle. The dotted vertical line represents the location of . (B) Phenotypes of strains in eight chemotherapeutics (12.5 µM amsacrine, 12.5 µM bleomycin, 2 µM bortezomib, 250 µM carmustine, 125 µM cisplatin, 62.5 µM etoposide, 300 µM puromycin, and 100 µM silver) are plotted as Tukey box plots with strain (y-axis) by relative median optical density (median.EXT, x-axis). Statistical significance was calculated for each strain pair (File S7). Significance of each deletion strain compared to its parental strain (ECA1132 to N2 and ECA1134 to CB4856) is shown above each strain pair and colored by the parent strain against which it was tested (ns = non-significant (p-value > 0.05); *, **, ***, and **** = significant (p-value < 0.05, 0.01, 0.001, or 0.0001, respectively).

Discussion

In this study, we identified overlapping QTL on the center of chromosome V that influence sensitivities to eight chemotherapeutic drugs. Because five of these drugs are known to cause double-strand DNA breaks, we hypothesized that this genomic region might be pleiotropic – a single shared genetic variant affects the responses to each drug. Because this variant might affect drug responses by regulating gene expression levels, we looked for the co-existence of drug-response QTL and expression QTL on chromosome V. We identified 2,540 eQTL and nine eQTL hotspots, including a region on the center of chromosome V. We calculated the mediation effect of all 49 genes with an eQTL that maps to this hotspot region and identified as a candidate gene whose expression influences the responses to several chemotherapeutics. We used CRISPR-Cas9-mediated deletion strains to empirically validate the role of in the chemotherapeutic response. In addition to bleomycin (Brady ), we discovered that responses to cisplatin, amsacrine, and carmustine are affected by expression. In this study, we found evidence that several overlapping QTL are representative of pleiotropy at the gene level and further elucidated the function of as a potential response to double-strand DNA break stress.

Mediation of drug-response QTL using gene expression to identify causal genes

Mediation analysis often suggests potential candidate genes that underlie different traits (Huang ; Sasaki ) and could be applied to drug responses. Using C. elegans strains and high-throughput assays, we can rapidly validate hypotheses generated by mediation analysis. Three of the eight chemotherapeutics that map to an overlapping drug-response QTL and were potentially mediated by were validated using targeted deletion strains.

Although mediation analysis provided moderate evidence that expression of could also play a role in sensitivity to etoposide and puromycin, we observed no experimental evidence of this relationship. Additionally, we have evidence that expression of might mediate response to carmustine. However, mediation analysis disagrees. The discrepancy between the mediation analysis and validation of causality using targeted deletion strains could be partially explained by one of several possibilities. First, different traits were measured in each experiment. The mediation analysis used traits measured over 96 hr of growth in drug conditions spanning two generations, but the causality test used traits measured over 48 hr of growth in drug conditions within one generation. Second, the precision of our mediation estimates was likely reduced by the poor quality drug traits for the set 1 RIAIL panel (Andersen ). Indeed, bortezomib, carmustine, etoposide, and puromycin did not map to the center of chromosome V using the set 1 RIAILs (Figure S5). Expression data for the set 2 RIAIL panel would likely generate more accurate mediation estimates, especially if data were collected using RNA sequencing to avoid the inherent reference bias of microarray data (Zhao ). Third, our mediation analysis was performed using expression data collected in control conditions and phenotype data collected in drug conditions. This analysis will only provide evidence of mediation if the baseline expression differences affect an individual’s response to the drug. Collecting expression data from drug-treated nematodes could help us learn more about how gene expression varies in response to treatment with the chemotherapeutic. Finally, as we only directly assessed the complete loss of in drug sensitivity, it is still possible that reduction of function (or change in function) caused by a single nucleotide variant or other structural variation in CB4856 could validate the role of in responses to these drugs.

This study demonstrates the power of pairing genome-wide linkage mapping of gene expression and drug response data using simple colocalization as well as more complex mediation analysis techniques. In addition to providing a resource for candidate gene prioritization within a QTL interval, mediation analysis can help to identify the mechanism by which genetic variation causes phenotypic differences. This type of approach could be even more powerful using genome-wide association (GWA) where the lower linkage disequilibrium between variants also has smaller confidence intervals in some genomic regions. Smaller intervals have fewer spurious overlapping eQTL, which could help to narrow the list of candidate genes. Although mediation analysis is only effective if a change in expression is observed and might not be useful for identifying effects from protein-coding variation, many current studies show that the majority of genetic variants associated with complex traits lie in regulatory regions (Hindorff ). Whole-genome expression analysis could provide the missing link to the identification of causal genes underlying complex traits.

New evidence for the pleiotropic function of scb-1

We identified eight chemotherapeutics with a QTL that mapped to a genomic region defined as a QTL hotspot on the center of chromosome V (Evans ). Multiple genes in close proximity, each regulating an aspect of cellular growth and fitness, might underlie each QTL independently. Alternatively, genetic variation within a single gene might regulate responses to multiple (or all) of the eight drugs tested, particularly if the gene is involved in drug transport or metabolism or if the drug mechanisms of action were shared (e.g., repair of double-strand DNA breaks). Expression of , a gene previously implicated in modulating responses to bleomycin, was found to reduce sensitivity to half of the drugs tested. This pleiotropic effect of provides new evidence for the function of the gene and possible molecular mechanisms underlying nematode drug responses. It is hypothesized that SCB-1 might function as a hydrolase that metabolizes compounds like bleomycin (Brady ) or somehow plays a role in the nematode stress response (Riedel ). Both hypotheses are consistent with our data, explaining why nematodes with low expression of are highly sensitive to the compound. Furthermore, all four of these chemotherapeutics, whose responses are mediated by expression of , are known to cause double-strand DNA breaks (Dorr 1992; Ketron ; Dasari and Tchounwou 2014; Nikolova ). Although the results for bortezomib, puromycin, and silver were inconclusive, we found no clear evidence that expression of dictates their responses. Together, these data suggest a potential role for specifically in response to stress induced by double-strand DNA breaks. However, the lack of sensitivity in etoposide, which also causes double-strand DNA breaks (Montecucco ), indicates that this response might be more complex.

The exact variant that causes the differential expression of is still unknown. Importantly, lies within an eQTL hotspot region where it is hypothesized that genetic variation at a single locus might regulate expression of the 49 genes with an eQTL in this region. It is possible that the same causal variant that regulates expression of could also underlie the QTL for the remaining four chemotherapeutics through differential expression of other genes. For example, mediation analysis for both bortezomib and etoposide indicated that expression variation of a dehydrogenase () may underlie their differential responses (File S12). Alternatively, the causal variants underlying these drug-response QTL might be distinct but physically linked in the genome. This result would suggest a cluster of genes essential for the nematode drug response. Overall, our study highlights the power of using mediation analysis to connect gene expression to organismal traits and describes a novel function for the pleiotropic gene .