The genetic basis of natural variation in a phoretic behavior

Phoresy is a widespread form of commensalism that facilitates dispersal of one species through an association with a more mobile second species. Dauer larvae of the nematode Caenorhabditis elegans exhibit a phoretic behavior called nictation, which could enable interactions with animals such as isopods or snails. Here, we show that natural C. elegans isolates differ in nictation. We use quantitative behavioral assays and linkage mapping to identify a genetic locus (nict-1) that mediates the phoretic interaction with terrestrial isopods. The nict-1 locus contains a Piwi-interacting small RNA (piRNA) cluster; we observe that the Piwi Argonaute PRG-1 is involved in the regulation of nictation. Additionally, this locus underlies a trade-off between offspring production and dispersal. Variation in the nict-1 locus contributes directly to differences in association between nematodes and terrestrial isopods in a laboratory assay. In summary, the piRNA-rich nict-1 locus could define a novel mechanism underlying phoretic interactions.

S pecies dispersal has been an important topic in evolutionary biology since Charles Darwin's era. He was fascinated by how one species facilitates the migration of another species. For example, he described a phoretic interaction between ducks and freshwater snails after observing just-hatched snails attached to a duck's foot, suggesting a dispersal mechanism for the wide-range distribution of snails 1 . Since that time, numerous reports of phoresy have accumulated [2][3][4] , but the genetic bases of these interactions remain elusive. To address the genetic underpinnings of a natural phoretic interaction, we investigated the association between the stress-resistant, long-lived dauer larvae 5 of the nematode Caenorhabditis elegans and terrestrial isopods 6 . Interspecific association and phoretic dispersal of dauer larvae are conserved among Caenorhabditis species 7 (Supplementary Movie 1), and this interaction is facilitated by a dauer-specific phoretic behavior called nictation where dauer larvae lift and wave their bodies, presumably to increase interactions with more mobile species 8,9 (Supplementary Movie 2).
Because natural variants are known to contribute to phenotypic differences in wild populations, quantitative genetic approaches have been used to successfully elucidate the molecular bases of the natural variation underlying a variety of C. elegans traits [10][11][12][13] . We investigated natural variation in the nictation behavior of C. elegans wild isolates, utilizing a previously established quantitative assay to measure the fraction of nictating dauers among a population of moving dauers 8,14 . We find significant differences among the strains, including two strains (N2 from Bristol, England and CB4856 from Hawaii, USA) for which recombinant inbred lines have been constructed to enable linkage-mapping approaches 15 . Using this mapping approach, we identified a quantitative trait locus (QTL) on chromosome IV that contributes to variation in nictation behavior. This locus contains a large cluster of Piwi-interacting small RNAs (piRNAs), which are all regulated by the Piwi Argonaute PRG-1 16 . Because we could not perturb single piRNAs out of many hundreds, we generated prg-1 mutant strains to test the causal connection of piRNAs to nictation behavior. These results show that the N2 strain contains piRNAs that likely inhibit genes that mediate nictation. We go on to demonstrate that this QTL underlies a phoretic behavior with terrestrial isopods and contributes to a potential fitness trade-off.

Results
Natural variation in nictation behavior. To establish whether natural strains of C. elegans vary in a potential dispersal behavior, we focused on nictation-a behavior where the long-lived dauer larvae lift and wave their bodies to increase interactions with larger, more mobile species. We used a previously established quantitative assay to measure the fraction of nictating dauers among a population of moving dauers 8,14 for a collection of genetically diverse wild strains ( Fig. 1a and Supplementary Data 1). We found that 12 strains differed significantly (one-way analysis of variance (ANOVA), p < 0.001) in nictation behaviors over a fourfold range, indicating that natural strains likely have diverse phoretic interactions mediated by nictation. We found that a smaller fraction of CB4856 dauers nictate than dauers of the N2 strain (Fig. 1a), and recombinant inbred lines have been constructed to enable linkage-mapping approaches 15 using these two strains. Additionally, N2 dauers show the fifth highest mean nictation fraction out of the 12 strains measured, indicating that differences between these two strains fall within the range of natural behavioral variation. Identification and fine mapping of the nict-1 QTL. We sought the genetic causes underlying this difference in nictation by measuring the behaviors of 186 recombinant inbred advanced intercross lines (RIAILs) derived from a cross between N2 and CB4856 ( Fig. 1b and Supplementary Data 2). Linkage mapping revealed a single significant QTL, which we named nict-1, on the right arm of chromosome IV (Fig. 1c). This locus does not overlap any of laboratory-derived loci, including npr-1, glp-5, and nath-10, that were discovered previously [17][18][19][20] . Recombinant strains with the N2 nict-1 genotype exhibited a higher nictation ratio compared to those strains with the CB4856 genotype ( Fig. 1d), consistent with the parental difference. c Genomic features of NIL-narrowed 73 kb nict-1 QTL interval are shown. Protein-coding genes (top) are colored orange (higher expression in the N2 nict-1 strain), blue (higher expression in CB4856 nict-1 strain), or gray (no expression difference between the nict-1 genotypes). The blue dot denotes a predicted functional effect in Y105C5A.1272. The locations of piRNAs and other genomics features (tRNAs, snoRNAs, etc.) are shown in gray in the middle and bottom rows. d Tukey box plots of normalized nictation fractions of N2 (orange), LJ1203 (nict-1 CB4856 > N2, blue), ECA586 (prg-1(ean28) in the N2 genetic background, gray), CB4856 (blue), LJ1213 (nict-1 N2 > CB4856, orange), ECA584 (prg-1(ean30) in the CB4856 genetic background, gray) are shown with data points plotted behind. The horizontal line is the median, and the box denotes the 25th to 75th quantiles of the data. The vertical line represents the 1.5 interquartile range. All strains were scored in 21 biological replicates, except for strain ECA584 where 19 biological replicates were scored. Each biological replicate was composed of three technical replicates To validate and to narrow the nict-1 QTL, we generated nearisogenic lines (NILs) by crossing the nict-1 QTL genomic region from N2 into the CB4856 genetic background, as well as performing the reciprocal cross. The phenotypes of these NIL strains confirmed the QTL effect in both genetic backgrounds, with the N2 nict-1 interval promoting higher nictation than the CB4856 nict-1 interval (Fig. 2a, b and Supplementary Data 3). By generating recombinants across the confirmed nict-1 region, we created nine additional NILs containing smaller N2 genomic regions in the CB4856 background. The differences in nictation behaviors of these NILs narrowed the nict-1 QTL to a 73 kb region (Fig. 2a, b and Supplementary Data 3).
piRNAs contribute to differences in nictation behavior. The nict-1 QTL genomic region contains four protein-coding genes, three pseudogenes, two non-coding RNA genes, and 289 annotated 21U piRNAs genes 21 (Fig. 2c). Using characterized sequence variation between the two parental strains 22-24 , we found that only one (Y105C5A.1272) of the four protein-coding genes contained variation between the N2 and CB4856 strains. The variant in Y105C5A.1272 encodes a putative serine-to-glycine change that could alter gene function (Supplementary Data 4). Additionally, we investigated gene expression of the four protein-coding genes and the three putative psuedogenes by quantitative RT-PCR. We found that the N2 strain expresses Y105C5A.1272 gene, whereas the CB4856 strain does not (Supplementary Data 5). We tested two independently generated deletion alleles of this gene in the N2 nict-1 genetic background ( Supplementary Fig. 2). Although each deletion removes a large proportion of the Y105C5A.1272 coding sequence, the two mutant strains did not show the same effect on nictation (Supplementary Fig. 1 and Supplementary Data 6), suggesting that Y105C5A.1272 does not underlie differences in nictation between the two strains. Of the other three protein-coding genes in the interval (Fig. 2c), we found that Y105C5A.15 had a detectable expression difference between the parental strains and the nict-1 NILs (Supplementary Data 5).
To test whether this gene expression difference could cause the nictation effect, we assayed a deletion allele of Y105C5A.15 in the N2 background and found that it does not have an effect on nictation (Supplementary Figs. 1 and 2 and Supplementary Data 6), suggesting that none of the protein-coding genes in the nict-1 interval play a role in nictation differences between the N2 and CB4856 strains.
Next, we investigated variation in pseudogene sequences and expression in both the N2 and CB4856 strains. We found that the three pseudogenes have variants that would eliminate gene function and remain psuedogenes in both genetic backgrounds. Additionally, expression of these three psuedogenes did not correlate with the differences in the nict-1 QTL genotype (Supplementary Data 5). These results led us to consider the 289 21U piRNA genes. Because specific perturbations of individual piRNA genes often do not cause loss of gene regulation and the 289 piRNA genes are distributed throughout the 73-kb nict-1 QTL region, single gene perturbations are not feasible. However, the more than 12,000 21 U piRNA genes, including those genes in this region, require the Piwi Argonaute encoding gene, prg-1 17 . Therefore, we globally perturbed piRNA functions by deleting exons one through seven of the prg-1 gene in both the N2 and CB4856 genetic backgrounds using the CRISPR/Cas9 genome-editing system [25][26][27] . Loss of prg-1 in the N2 genetic background resulted in a nictation fraction similar to the CB4856 strain (Fig. 2d). Furthermore, loss of prg-1 in the CB4856 genetic background had a similar nictation fraction as the CB4856 strain, implicating N2-specific piRNAs in the control of this phoretic behavior. Additionally, an independent prg-1 lossof-function allele in the N2 genetic background had a similar effect on nictation ( Supplementary Fig. 3 and Supplementary Data 6). These small RNA genes regulate the expression of transposons and also endogenous genes 28 . Any one or multiple of these targeted genes could play a role in the regulation of nictation behavior. Because expression of prg-1 is restricted to the germline 16 , it remains elusive how differences in piRNAs can underlie differences in nictation behaviors. The mean N2 nict-1 QTL allele frequencies before (origin) and after (destination) transmission competition assays are plotted. Error bars are standard deviations. c The experimental scheme of the transmission competition assay is shown. Transmission competition assay chambers contain two nematode culture plates. The "origin" plate contained a mixture of dauer animals from the CB4856 and LJ1213 (nict-1 N2 > CB4856) strains and was covered with medical gauze to facilitate nictation. The "destination" plate was placed 1-2 cm away with no nematodes on this plate. Terrestrial isopods were added to the chamber and allowed to roam freely for 24 h. Nematodes depend on the isopods for transfer to the "destination" plate from the "origin" plate because no transfer was observed without isopods presence. Eleven different biological replicates were scored nict-1 controls a trade-off between nictation and reproduction.
The nict-1 QTL likely contributes to different fitness consequences for wild C. elegans strains. The NIL strains had brood sizes that were significantly different from their respective parental genetic backgrounds ( Fig. 3a and Supplementary Data 7). The offspring production of nict-1 from N2 crossed into the CB4856 background (LJ1213) was lower than that of the CB4856 parent, and the brood size of nict-1 from CB4856 crossed into the N2 background (LJ1203) was higher than N2. Therefore, the nict-1 QTL likely confers a trade-off between nictation and reproduction with the N2 version, promoting nictation but inhibiting offspring production. This result implies that the N2 strain has a functional piRNA that acts to promote nictation and the same or an independent piRNA that acts to inhibit reproduction. Alternatively, a linked N2 variant in a proteincoding gene could exert an opposing effect on reproduction in the N2 nict-1 genotype. Next, we analyzed a correlation between the nictation fraction and previously determined 22 offspring production of the wild strains ( Supplementary Fig. 4). We did not find a strong correlation between these two traits (r = −0.1666), suggesting that the trade-off between nictation and reproduction might be specific to the nict-1 QTL between the N2 and CB4856 strains.
Natural populations of C. elegans exhibit local genetic diversity 29 , suggesting that intraspecific competition occurs among wild strains in the natural environment. To examine how quantitative variation in nictation behavior relates to the dispersal of C. elegans in a competitive environment, we developed a transmission competition assay to quantitatively analyze isopod phoresy (Fig. 3b, c). To mimic this natural phoretic interaction, we used the terrestrial isopod, Porcellio scaber, so that we can test whether C. elegans can hitchhike to a more favorable environment. Isopods were previously reported as a natural carrier of C. elegans 6, 7 , and dauers can readily attach to them. We induced dauers from two different nict-1 genotypes in the same culture to determine whether the N2 nict-1 genotype, which confers a higher nictation fraction, could confer a higher transmission frequency than the CB4856 nict-1 genotype. Because of differences in offspring production caused by the trade-off discussed above, the culture contained four times as many dauers of the CB4856 nict-1 genotype than the N2 nict-1 genotype ( Fig. 3b and Supplementary Data 8). Despite this disparity on the origin plate, we found that dauers with the N2 nict-1 genotype were transported to the destination plate at a rate four times higher than the rate of the CB4856 nict-1 genotype. This transport was dependent on the presence of terrestrial isopods. These results suggest that variation in nict-1 underlies this natural phoretic interaction. Furthermore, given the potential intraspecific competition in the natural environment, our results demonstrate that the nict-1 QTL could control the ability of strains to colonize new bacteria-rich environments.

Discussion
We propose that the nict-1 QTL controls a hitchhiking behavior and phoretic dispersal, as well as a trade-off between dispersal and reproduction, two traits necessary for the survival and evolution of the C. elegans species. Behavioral diversity arises from the gene-environment interface where a biological system confronts an ever-changing natural niche, and multiple genes are involved in such dynamic interactions 30 . Although this study identifies a novel difference between the laboratory strain and a wild strain, the observed difference in nictation does not appear to be laboratory-derived because the N2 strain is not exceptional in this behavior, as it is for other traits optimized in the laboratory 31 . Our study provides strong evidence that the nictation behavior is controlled by differences in regulatory small RNAs, which could act via a large number of distributed effects to create robust environmental adaptations 32 . Involvement of the piRNA pathway in neuronal plasticity is suggested from studies of the Aplysia central nervous system and the mammalian brain, implying conserved regulatory roles among metazoans [33][34][35][36] . Studies of a large collection of diverse wild strains are required to understand the broad applicability of the piRNA-rich nict-1 QTL effect across the C. elegans population. This proposed molecular mechanism of behavioral control enables future studies of this phoretic interaction and its precise genetic causes across C. elegans and related nematode species.
Dauer induction and nictation assays. Ten to twenty L4 larvae or young adults were transferred to synthetic pheromone plates containing agar (10 g/L), agarose (7 g/L), NaCl (2 g/L), KH 2 PO 4 (3 g/L), K 2 HPO 4 (0.5 g/L), cholesterol (8 mg/L), and synthetic pheromone-ascaroside 1, 2, 3 (2 mg/L each) 38, 39 seeded with Escherichia coli OP50 at 25°C for dauer induction 8,14 . After 4 days, dauers were morphologically identified by their dark intestines and radially constricted bodies. Micro-dirt chips were made by pouring 3.5% agar solution onto polydimethylsiloxane (PDMS) mold 14 . Solidified agar micro-dirt chip was detached from the PDMS mold and dried for 90 min at 37°C. More than 30 dauers were collected by glass capillary (Kimble Chase Life Science and Research Products LLC.) using M9 buffer and mounted on a micro-dirt chip. After 10-30 min, when dauers were actively moving, a fraction of nictating dauers among moving dauers on a micro-dirt chip was measured three times consecutively, and the mean value of three technical replicates was represented as the nictation fraction for that biological replicate. The mean nictation fraction was calculated from the independent biological replicates of multiple nictation assays. We calculated the normalized nictation value by fitting a linear model, nictation fraction~assay date. Number of independent biological replicates for each assay is described in figure legends.
Quantitative genetic analyses. The nictation fractions of 186 RIAIL strains from an advanced intercross between N2 and CB4856 were scored as described above. The phenotype data and genotype data were entered into R and scaled to have a mean of zero and a variance of one for linkage analysis. QTLs were detected by calculating logarithm of odds (LOD) scores for each marker and each trait as −n(ln(1−r 2 )/2ln(10)), where r is the Pearson correlation coefficient between RIAIL genotypes at the marker and phenotype trait values 40,41 . We randomly permuted the phenotype values of each RIAIL while maintaining correlation structure among phenotypes 1000 times to estimate significance empirically. Confidence intervals were defined as the regions contained within a 1.5 LOD drop from the maximum LOD score.
NIL construction. LJ1203 snuIR3 (nict-1 QTL, CB4856 > N2) was made from QX4, a RIAIL with the CB4856 nict-1 QTL interval on an otherwise N2 chromosome IV. QX4 was backcrossed to the N2 parent six times while selecting for the CB4856 nict-1 QTL by single-nucleotide polymorphism (SNP) genotyping. LJ1204 snuIR4 (nict-1 QTL, N2 > CB4856) was made from QX162, a RIAIL with the N2 nict-1 QTL interval in an otherwise CB4856 chromosome IV. QX162 was backcrossed to Analysis of genes in the QTL interval. C. elegans protein-coding and non-coding RNA genes were extracted using a custom Python script (Wormbase GFF version WS248). Variants within the QTL interval that distinguish N2 from CB4856 and that are predicted to have "moderate" or "high" functional effects were identified using the R cegwas package 37,43 and were plotted with respect to gene models.
Brood size assays. Twenty young adults of each strain were transferred on 10 cm NGM plates seeded with OP50 and allowed to lay embryos. After 1-2 h, adult animals were removed and fewer than 100 embryos were grown at 20°C. After 2-3 days, 10 L4 larvae or young adults were singled and transferred to 55 mm NGM lite agar plates 44 seeded with OP50 at 25°C. Animals were transferred to new plates every 24 h for 3 days. The number of offspring on each plate was counted manually after 2 days from transfer. Normalized fecundity was calculated using a linear model, brood size~assay date.
Transmission competition assays. Ten synchronized young adults from CB4856 and LJ1213 (nict-1 N2 > CB4856) were placed together on a single pheromone plate, where dauers were induced by the same method as described above. After 4 days, when hundreds of dauers were observed, more than 20 dauers were singly collected by glass capillary and genotyped at the haw64888 SNP (IV, 15,696,457) in the nict-1 QTL region using single-animal lysis and a PCR-based DraI RFLP assay 42 . After that step, medical gauze was mounted on the pheromone plate as a platform to enable nictation of the remaining dauer larvae. The plate with nictating dauers and a new plate with OP50 food were placed 1-2 cm apart within a small plastic box. Four to six Porcellio scaber terrestrial isopods, collected from Seoul National University campus, were also placed in the box. The isopods moved freely for 24 h. We genotyped the transferred animals at the haw64888 SNP (IV, 15,696,457) in the nict-1 QTL using the same method as described above.