A nematode-specific gene underlies bleomycin-response variation in Caenorhabditis elegans

Bleomycin is a powerful chemotherapeutic drug used to treat a variety of cancers. However, individual patients vary in their responses to bleomycin. The identification of genetic differences that underlie this response variation could improve treatment outcomes by tailoring bleomycin dosages to each patient. We used the model organism Caenorhabditis elegans to identify genetic determinants of bleomycin-response differences by performing linkage mapping on recombinants derived from a cross between the laboratory strain (N2) and a wild strain (CB4856). This approach identified a small genomic region on chromosome V that underlies bleomycin-response variation. Using near-isogenic lines and strains with CRISPR-Cas9 mediated deletions and allele replacements, we discovered that a novel nematode-specific gene (scb-1) is required for bleomycin resistance. Although the mechanism by which this gene causes variation in bleomycin responses is unknown, we suggest that a rare variant present in the CB4856 strain might cause differences in the potential stress-response function of scb-1 between the N2 and CB4856 strains, thereby leading to differences in bleomycin resistance. Article summary We performed linkage mapping on a panel of recombinant lines generated between two genetically divergent strains of Caenorhabditis elegans and identified a bleomycin-response QTL. We generated CRISPR-Cas9 deletions and reciprocal allele-replacement strains for all six candidate genes across the QTL confidence interval. Deletions of one gene, H19N07.3, caused increased bleomycin sensitivity in both divergent genetic backgrounds. This gene might act in stress responses and detoxification in nematodes. We further compared our linkage mapping to a genome-wide association mapping and showed that a rare expression variant in the CB4856 strain likely underlies bleomycin-response differences.

individual patients vary in their responses to bleomycin. The identification of genetic differences 48 that underlie this response variation could improve treatment outcomes by tailoring bleomycin 49 dosages to each patient. We used the model organism Caenorhabditis elegans to identify genetic 50 determinants of bleomycin-response differences by performing linkage mapping on 51 recombinants derived from a cross between the laboratory strain (N2) and a wild strain 52 (CB4856). This approach identified a small genomic region on chromosome V that underlies 53 bleomycin-response variation. Using near-isogenic lines and strains with CRISPR-Cas9 54 mediated deletions and allele replacements, we discovered that a novel nematode-specific gene 55 (scb-1) is required for bleomycin resistance. Although the mechanism by which this gene causes 56 variation in bleomycin responses is unknown, we suggest that a rare variant present in the 57 CB4856 strain might cause differences in the potential stress-response function of scb-1 between 58 the N2 and CB4856 strains, thereby leading to differences in bleomycin resistance. 59

INTRODUCTION
was extracted from a previously established variant dataset (Cook et al. 2016  2018). A 5% genome-wide error rate was calculated by permuting phenotype and genotype data 213 1,000 times. Mappings were repeated iteratively, each time using the marker with the highest 214 LOD score as a cofactor until no significant QTL were detected. The annotate_lods function was 215 used to calculate the percent of variance in RIAIL phenotypes explained by each QTL and a 95% 216 confidence interval around each peak marker, defined by any marker on the same chromosome 217 within a 1.5-LOD drop from the peak marker. 218 11 of 48 219

Generation of near-isogenic lines (NILs): A near-isogenic line (NIL) is genetically identical to 220
another strain aside from a small genomic region that is derived from an alternate strain. NILs 221 are used to test the effect of modifications to particular genomic regions in a consistent genetic 222 background. To make each NIL, males and hermaphrodites of the desired RIAIL (QX131 for 223 ECA230 and QX450 for ECA232) and parental background (CB4856 for ECA230 and N2 for 224 ECA232) were crossed in bulk, then male progeny were crossed to the parental strain in bulk for 225 another generation. For each NIL, eight single-parent crosses were performed followed by six 226 generations of propagating isogenic lines to ensure homozygosity of the genome. For each cross, 227 PCR was used to select non-recombinant progeny genotypes within the introgressed region by 228 amplifying insertion-deletion (indel) variants between the N2 and CB4856 genotypes on the left 229 and right side of the introgressed region. NIL strains ECA411 and ECA528 were generated by 230 crossing ECA230 and CB4856 in bulk. Heterozygous F1 hermaphrodites were crossed to 231 CB4856 males and the F2 L4 hermaphrodites were placed into individual wells of a 96-well 232 plate with K medium and 5 mg/mL bacterial lysate and grown to starvation. After starvation, 233 each well of the 96-well plate was genotyped to identify recombinants in the desired region. 234 Recombinant strains were plated onto 6 cm plates and individual hermaphrodites were 235 propagated for several generations to ensure homozygosity across the genome. NIL strains were 236 whole-genome sequenced as described above to confirm their genotypes (File S5, File S6). 237 Reagents used to generate all NIL strains and a summary of each introgressed region are detailed 238 in the Supplementary Information. 239 240 Two-factor genome scan: We performed a two-factor genome scan to identify potentially 241 epistatic loci that might contribute to traits of interest (either bleomycin responses or gene-242 expression levels). We used the scantwo function in the R qtl package to perform this analysis. 243 Each pairwise combination of loci was tested for correlations with trait variation in the RIAILs. 244 The summary of each scantwo output includes the maximum interactive LOD score for each pair 245 of chromosomes. To determine a threshold for significant interactions, we performed 1000 246 permutations of the scantwo analysis and selected the 95 th percentile from these permutations. 247 For the bleomycin-response variation scantwo, the significant interaction threshold was a LOD 248 score of 4.08. The significant interaction threshold for the scb-1 expression variation scantwo 249 was 4.05. The bleomycin-response variation scantwo summary is available as File S8, and the 250 scantwo summary for scb-1 expression variation is available as File S17. 251 252 Identification of genes with non-synonymous variants: The get_db function within the cegwas 253 R package was used to query WormBase build WS245 for genes in the QTL confidence interval 254 (V:11,042,745-11,189,364). Our initial linkage mappings used the 1,454 marker set (Rockman 255 and Kruglyak 2009) and had a QTL confidence interval larger than the interval found using the 256 whole-genome marker set (described above). Because this confidence interval was larger and 257 more conservative, we kept it for subsequent testing of candidate genes. This expanded interval 258 contained an additional 20 kb on either side of the whole-genome marker set confidence interval 259 (File S9). The snpeff function within the cegwas R package was used to identify variants within 260 the region of interest (File S11  . We performed linkage mapping as described above on these expression data and 315 identified significant peaks for 3,298 probes (File S13). For the H19N07.3 probe 316 (A_12_P104350), the full annotated-LOD data frame is provided (File S14). 317 318 Hemizygosity high-throughput assay: Hemizygosity assays were used to test the effect of zero 319 (two deletion alleles), one (one deletion allele and one wild-type allele), or two (two wild-type 320 alleles) functional copies of a gene product on bleomycin responses. If the phenotype is affected 321 by gene function, one would expect to see bleomycin responses scale with the number of 322 functional alleles of the gene present in each strain. For each heterozygous/hemizygous 323 genotype, 30 hermaphrodites and 60 males were placed onto each of four 6 cm plates and 324 allowed to mate for 48 hours. The same process was completed for homozygous strains to 325 remove biases introduced by the presence of male progeny in the assay. Mated hermaphrodites 326 were transferred to a clean 6 cm plate and allowed to lay embryos for eight hours. After the egg 327 lay period, adults were manually removed from egg lay plates, and embryos were washed into 15 328 mL conicals using M9 and a combination of pasteur pipette rinsing and scraping with plastic 329 inoculation loops. Embryos were rinsed and centrifuged four times with M9 before being 330 resuspended in K medium and 50 µM kanamycin. Embryos hatched and arrested in the L1 larval 331 stage in conicals overnight at 20º with shaking. The next morning, 50 L1 larvae were sorted into 332 each well of a 96-well plate containing K medium, 10 mg/mL bacterial lysate, 50 µM 333 kanamycin, and either 1% distilled water or 1% distilled water plus bleomycin using the 334 BIOSORT (dose-response assay, Figure S10, File S1) or by titering (hemizygosity assays, 335 stop sites for all genes in the region. File S11 reports predicted non-synonymous variants 407 between the N2 and CB4856 strains in the region. File S12 is derived from the Rockman et al. 408 2010 RIAIL microarray expression data, and reports the expression measurements for each of the 409 13,107 microarray probes across 209 RIAILs. File S13 contains all significant QTL identified by 410 linkage mapping of File S12 data. File S14 contains the annotated linkage mapping of the 411 H19N07.3 expression data. File S15 reports the H19N07.3 expression and residual median 412 optical density for strains of the RIAIL panel that were assayed for both of those traits. File S16 413 contains H19N07.3 RNA-seq expression data for populations of young adults of N2 and 414 CB4856. File S17 is a summary of the scantwo analysis for H19N07.3 expression in the RIAILs 415 and reports the maximum interaction LOD score for each chromosome pair. File S18 contains 416 control-regressed phenotypic data for all wild isolates assayed in response to bleomycin. File 417 S19 contains genome-wide association mapping for the phenotypes in File S18. File S20 418 contains genotype information for each strain measured in File S18 across all variants within the 419

Genetic differences underlie bleomycin-response variation 427
Bleomycin causes double-stranded DNA breaks, which ultimately lead to cytotoxicity of rapidly 428 dividing cell populations. Therefore, exposure to bleomycin can affect the development of C. 429 elegans larvae as well as germ-cell production of adult animals. We used a high-throughput 430 assay (HTA) to measure the effects of bleomycin on development and brood size (Figure S1, 431 Materials and Methods). To determine the concentration of bleomycin that would maximize 432 among-strain while minimizing within-strain phenotypic variation, we used the HTA to perform 433 a dose-response assay. We assessed bleomycin responses for four divergent strains (N2, 434 CB4856, JU258, and DL238) across each of 26 HTA traits (File S1, Figure S2). Lee et al. 2017). We used these RIAILs to identify genetic variants that contribute to differential 449 bleomycin responses between the N2 and CB4856 strains. Using our HTA, we measured each of 450 the 26 fitness parameters for 249 RIAILs (Materials and Methods, File S3). Correlations 451 between each pairwise combination of the 26 HTA measurements revealed several clusters of 452 highly correlated traits ( Figure S3). Therefore, the summary statistics measured by the 453 BIOSORT should not be considered independent traits for linkage mapping. We selected median 454 optical density (median.EXT) for future analyses, which is related to both animal length and 455 optical extinction, because this trait was highly correlated with many of the 26 HTA traits and 456 was highly heritable (H 2 = 0.73, File S2). 457 458

The QTL on the center of chromosome V strongly impacts bleomycin response 459
We performed linkage mapping on the residual median optical density measurements in 460 bleomycin ( Figure S4, Figure S5, File S4, Materials and Methods) and identified four 461 significant quantitative trait loci (QTL, Figure 1A). The QTL on the center of chromosome V 462 was highly significant (explained 43.58% of the total variation and 55.60% of the genetic 463 variation) with a LOD score of 32.57, and it was detected for 25 of the 26 HTA traits (Figure S5,  464 File S4). The QTL 95% confidence interval was approximately 147 kb. Strains with the N2 465 allele at the peak marker had a lower median optical density in bleomycin and were interpreted 466 to be more sensitive than those RIAILs with the CB4856 allele ( Figure 1B). 467

468
We isolated this QTL in a controlled genetic background by generating near-isogenic lines 469 (NILs, File S5, File S6) that each contain a genetic background derived from either the N2 or 470 22 of 48 CB4856 strain and a region of chromosome V from the opposite parental genotype. We used the 471 HTA to test these strains in response to bleomycin (File S3). The NIL with the N2 genotype 472 across this QTL introgressed into the CB4856 background (ECA230) was statistically more 473 sensitive to bleomycin than CB4856 (Figure S6, File S7, p = 1.3e-14, Tukey HSD). This 474 phenotype indicated that the N2 genotype within the introgressed region (which includes the 475 QTL confidence interval) confers sensitivity to bleomycin. However, the reciprocal NIL with the 476 CB4856 locus introgressed into the N2 background (ECA232) had a bleomycin-response 477 phenotype that was not significantly different from the N2 strain ( Figure S6, File S7, p = 0.053, 478 Tukey HSD), suggesting that interacting loci could underlie bleomycin responses in a 479 background-dependent manner. We performed a two-factor genome scan to map potential 480 epistatic loci but did not identify a significant interaction between the QTL on chromosome V 481 and other loci (Figure S7, File S8). However, the failure to detect significant interacting QTL 482 could be because we have too few recombinant strains or because too few replicates of each 483 RIAIL were phenotyped. Alternatively, more than two loci might underlie the transgressive 484 phenotype of ECA230 and a two-factor genome scan might not be able to capture this 485 complexity. 486 487 Nonetheless, because ECA230 recapitulated the expected QTL phenotype, we generated two 488 NILs (ECA411 and ECA528) that narrowed this introgressed region to more precisely locate the 489 causal variant (File S5, File S6). In addition, the N2 region on the left of chromosome V was 490 removed from both NIL strains to ensure that this region of introgression did not underlie the 491 phenotypic difference between ECA230 and CB4856. The genotypes of ECA411 and ECA528 492 differ in a small region of chromosome V that includes the QTL confidence interval (Figure 2, 493 File S5, File S6). Both of these strains were more sensitive to bleomycin than the background 494 parental strain, CB4856. This result could suggest that the introgressed region shared between 495 these strains, which does not include the QTL, conferred some bleomycin-response variation 496 between the N2 and CB4856 strains (Figure 2). Alternatively, the hypersensitivity of these NILs 497 could suggest the presence of Dobzhansky-Muller incompatibilities between the N2 and CB4856 498 genotypes (Snoek et al. 2014) that might affect stress responses of the NILs. However, ECA528 499 was much more sensitive to bleomycin than ECA411 (Figure 2, File S7). Because ECA528 has 500 the N2 genotype across the QTL region and ECA411 has the CB4856 genotype, these results 501 suggest that the QTL genotype strongly affects bleomycin sensitivity (Figure 2 Because the recombination rate in the centers of C. elegans chromosomes is lower than 509 chromosome arms (Rockman and Kruglyak 2009), it was difficult to generate additional NILs to 510 narrow the QTL region further. Therefore, we took a targeted approach and created CRISPR-511 Cas9 directed modifications of candidate genes in the region. The 147 kb confidence interval on 512 chromosome V contains 93 genes, including pseudogenes, piRNA, miRNA, ncRNA, and 513 protein-coding genes (File S9). Given the narrow confidence interval, we expanded our search to 514 include an additional 20 kb on each side of the 147 kb interval (Materials and Methods). Of the 515 118 genes included in the wider region, five genes, C45B11.8, C45B11.6, jmjd-5, srg-42, and 516 cnc-10, contain predicted non-synonymous variants between the N2 and CB4856 strains ( Figure  517 3A, File S10, File S11). These variants could cause differential bleomycin sensitivity between 518 the N2 and CB4856 strains. 519 520 To test these genes in bleomycin-response variation, we systematically deleted each of the 521 candidate genes in both the N2 and CB4856 backgrounds. We used CRISPR-Cas9 mediated 522 genome editing to generate two independent deletion alleles of each gene in each genetic 523 background to reduce the possibility that off-target mutations could cause phenotypic differences 524 (Materials and Methods, Supplemental Information). We tested the bleomycin response of each 525 deletion allele in comparison to the N2 and CB4856 parental strains ( Figure 3B, File S3). The 526 deletion alleles of C45B11.8, C45B11.6, srg-42, and cnc-10 each had a bleomycin response 527 similar to the respective parent genetic background, which suggested that the functions of each 528 of these genes did not affect bleomycin responses ( Figure 3B, File S7, p > 0.05, Tukey HSD). 529 By contrast, the jmjd-5 deletion alleles in the N2 and the CB4856 backgrounds were each more 530 resistant to bleomycin than their respective parental strains ( Figure 3B, File S7, ECA1047 vs. 531 CB4856 p = 3.8e-10, ECA1048 vs. CB4856 p = 0.026, ECA1051 vs. N2 p = 7.4e-4, ECA1052 532 vs. N2 p = 2.9e-6, Tukey HSD). However, we also noted that these strains were more sensitive in 533 the control condition than other deletion strains ( Figure S8, File S1, File S7). Therefore, the 534 relative increased bleomycin resistance observed in the jmjd-5 deletion strains could be caused 535 by their increased sensitivity in the control condition. 536

537
We tested if the non-synonymous variant in jmjd-5 between the N2 and CB4856 strains caused 538 bleomycin-response differences. At residue 338 of JMJD-5, the N2 strain has a proline, whereas 539 25 of 48 the CB4856 strain has a serine (S338P, Figure S9A). We used CRISPR-Cas9 to generate 540 reciprocal allele replacements of the jmjd-5 single-nucleotide polymorphism that encodes the 541 putative amino-acid change in the N2 background jmjd-5(N2 to CB4856) and in the CB4856 542 background jmjd-5(CB4856 to N2) (Materials and Methods, Supplemental Information). We 543 created two independent allele replacements in each genetic background and measured each 544 strain for bleomycin-response differences as compared to the parental strains ( Figure S9B, File  545   S3). Although the allele-replacement strains with the CB4856 allele in the N2 genetic 546 background jmjd-5(N2 to CB4856) were significantly different from the N2 parental strain, these 547 strains were more sensitive to bleomycin than N2 ( Figure S9B, File S7, ECA576 vs. N2 p = 548 0.006, ECA577 vs. N2 p = 1.6e-6, Tukey HSD). This increased sensitivity was unexpected, 549 because the CB4856 allele at the jmjd-5 locus should confer resistance. However, the NIL with 550 the CB4856 genotype across the QTL was not different from the N2 parental strain (ECA232 in 551 Figure S6, File S7), suggesting that the QTL might only confer increased sensitivity when the 552 N2 allele is in the CB4856 background. Therefore, it remained unclear whether an allele 553 replacement of jmjd-5 in the N2 parental background could confer resistance. Neither of the two 554 strains with the N2 allele in the CB4856 background, jmjd-5(CB4856 to N2), conferred a 555 significantly more sensitive phenotype than the CB4856 parental strain (Figure S9B, File S7). 556 Given that the QTL explained 43.58% of phenotypic variation among the RIAILs, the causal 557 variant should have a clear impact on bleomycin response. Additionally, the NILs with the N2 558 allele at the QTL introgressed into the CB4856 background displayed a significant increase in 559 bleomycin sensitivity compared to the parental CB4856 strain (Figure S6, File S7). Taken 560 together, the phenotypes of the reciprocal allele-replacement strains showed that the amino-acid 561 change in JMJD-5 likely does not underlie bleomycin-response variation between the N2 and 562 26 of 48 CB4856 strains, although deletion of this gene did cause resistance to bleomycin regardless of 563 the genetic background. We performed a reciprocal hemizygosity assay to test if natural variation 564 in jmjd-5 function affected bleomycin responses ( Figure S10, Figure S11, Figure S12). The 565 results of this assay supported the previously identified increase in bleomycin resistance of 566 homozygous jmjd-5 deletions in both parental backgrounds ( Figure S11, File S7, p < 0.05, 567 Tukey HSD), which again might be caused by an increased sensitivity in the control condition 568 ( Figure S12). However, the increases in bleomycin resistance between each jmjd-5 deletion 569 strain and the strain with the same genetic background were similar, and the reciprocal 570 hemizygous strains show equivalent bleomycin responses ( Figure S11). Taken together, these 571 results suggest that natural variation in jmjd-5 function does not underlie this QTL. 572

573
The nematode-specific gene H19N07.3 impacts bleomycin variation 574 Because none of the genes with a non-synonymous variant between the N2 and CB4856 strains 575 explained the QTL, we explored other ways in which natural variation could impact bleomycin 576 responses. We used the 13,001 SNPs to perform linkage mapping on the gene expression data of 577 the RIAILs and identified 4,326 expression QTL (eQTL) across the genome (Rockman et al. 578 2010, File S12, File S13, Materials and Methods). Of the 118 genes in the 187 kb surrounding 579 the bleomycin-response QTL, expression for 50 genes were measured in the previous microarray 580 study. We identified a significant eQTL for eight of these 50 genes, four of which mapped to 581 chromosome V (File S13). eQTL for two of those four genes, H19N07.3 and cnc-10, mapped to 582 the center of chromosome V and overlapped with the bleomycin-response QTL. Because cnc-10 583 did not underlie bleomycin response variation (Figure 3B) (Figure 4A and B, File S14). The length 586 of animals and expression of scb-1 was correlated in the RIAIL strains ( Figure 4C, File S15, 587 r 2 = 0.61, p < 9.5e-13 Spearman's correlation). Although this gene does not have a non-588 synonymous variant between the N2 and CB4856 strains, natural variation in gene expression 589 could impact bleomycin responses. Each H19N07.3 deletion strain was more sensitive to bleomycin than the respective parental 595 strain ( Figure 5, File S3, File S7, ECA1133 vs. CB4856 p < 1.4e-14, ECA1134 vs. CB4856, p < 596 1.4e-14, ECA1132 vs. N2, p = 6.9e-5, ECA1135 vs. N2, p = 0.006, Tukey HSD). These results 597 suggest that H19N07.3 function is required for resistance to bleomycin. Therefore, we renamed 598 this gene scb-1 for sensitivity to the chemotherapeutic bleomycin. Unlike with the jmjd-5 599 deletion strains, the scb-1 deletion strains had no significant differences in the control condition 600 (File S7). Therefore, the bleomycin sensitivity of the scb-1 deletion strains were not caused by 601 control-condition phenotypes. 602 603 Because an scb-1 non-synonymous variant has not been identified between the N2 and CB4856 604 strains, changes to protein function likely do not cause bleomycin response differences. RIAILs 605 with the CB4856 allele at the QTL peak marker have increased expression of scb-1 and 606 increased bleomycin resistance compared to RIAILs with the N2 allele (Figure 1, Figure 4). 607 Therefore, scb-1 expression differences might cause the bleomycin-response variation between 608 28 of 48 the parental strains. We performed RNA-seq of the N2 and CB4856 strains to assess scb-1 609 expression differences between the parental strains and did not identify a significant increase in 610 expression in the CB4856 strain ( Figure S13, File S16, p = 0.20, Wald test; p = 0.17, likelihood 611 ratio test). This result could be caused by the low sample size (n = 3) in the RNA-seq 612 experiment, or the RIAIL strains could have a novel variant that arose during strain construction 613 that causes scb-1 expression variation. Alternatively, the expression difference observed in the 614 RIAIL strains could be attributed to epistatic loci. We performed a two-factor genome scan to 615 identify epistatic loci that underlie scb-1 expression variation in the RIAILs and identified two 616 significant interactions: one between loci on chromosomes IV and X and another between loci on 617 chromosomes II and V ( Figure S14, File S17). This result might suggest that epistatic loci 618 underlie scb-1 expression variation in the RIAILs and could explain why scb-1 expression is not 619 variable in the parental strains. 620

621
To test the role of natural variation in scb-1 function, we performed a reciprocal hemizygosity 622 test in control and bleomycin conditions (Figure 6, File S3). These results matched the increase 623 in sensitivity of homozygous deletions in both parental backgrounds observed previously. The 624 hemizygous strain with the CB4856 allele of scb-1 had a bleomycin phenotype similar to the 625 heterozygous strain, whereas the hemizygous strain with the N2 allele of scb-1 was more 626 sensitive to bleomycin than the heterozygous strain. Taken together, these results suggest that 627 natural variation in scb-1 function underlies the bleomycin-response difference between the N2 628 and CB4856 strains. 629 630 Differences in scb-1 function might be regulated by a rare variant 631