Ras/MAPK Modifier Loci Revealed by eQTL in Caenorhabditis elegans

The oncogenic Ras/MAPK pathway is evolutionarily conserved across metazoans. Yet, almost all our knowledge on this pathway comes from studies using single genetic backgrounds, whereas mutational effects can be highly background dependent. Therefore, we lack insight in the interplay between genetic backgrounds and the Ras/MAPK-signaling pathway. Here, we used a Caenorhabditis elegans RIL population containing a gain-of-function mutation in the Ras/MAPK-pathway gene let-60 and measured how gene expression regulation is affected by this mutation. We mapped eQTL and found that the majority (∼73%) of the 1516 detected cis-eQTL were not specific for the let-60 mutation, whereas most (∼76%) of the 898 detected trans-eQTL were associated with the let-60 mutation. We detected six eQTL trans-bands specific for the interaction between the genetic background and the mutation, one of which colocalized with the polymorphic Ras/MAPK modifier amx-2. Comparison between transgenic lines expressing allelic variants of amx-2 showed the involvement of amx-2 in 79% of the trans-eQTL for genes mapping to this trans-band. Together, our results have revealed hidden loci affecting Ras/MAPK signaling using sensitized backgrounds in C. elegans. These loci harbor putative polymorphic modifier genes that would not have been detected using mutant screens in single genetic backgrounds.

Here we go beyond mutant screens of let-60(gf) in a single genetic background in C. elegans by incorporating the effects of multiple genetic backgrounds. This provides the opportunity to explore the genetic variation for identifying novel modifiers of the Ras/MAPK pathway. We recently mapped modifiers affecting Ras/MAPK signaling associated with vulval development in a population of recombinant inbred lines derived from a cross between wild-type Hawaiian CB4856 and Bristol N2 (Schmid et al. 2015). The lines were sensitized by introgression of the G13E gain-of-function mutation (n1046) in the Ras gene let-60 [mutation introgression recombinant inbred lines (miRILs)]. This mutation is analogous to mutations causing excess cell division in human tumors (Kyriakakis et al. 2015). Hawaiian CB4856 males were crossed with Bristol N2 let-60 mutants. Random segregation of the two parental genomes was allowed, except for the let-60 mutation which was kept homozygous from the F2 generation onwards. After 10 generations of self-fertilization, to drive all regions to homozygosity, independent miRILs were successfully obtained, each carrying the mutation. Quantitative trait loci (QTL) analysis on the vulval index (VI) of these miRILs in combination with allele-swapping experiments revealed the polymorphic monoamine oxidase A (MAOA) gene amx-2 as a negative regulator of Ras/MAPK signaling (Schmid et al. 2015).
Here, we extended the study by Schmid et al. (2015) by studying the transcriptional architecture of the same let-60(gf)-sensitized miRILs. First, we measured the transcriptome of 33 miRILs using microarrays. Second, to gain further insight into the underlying molecular mechanisms, pathways, and new modifiers, we mapped gene expression QTL (eQTL). eQTL are polymorphic loci underlying variation in gene transcript abundances and can be used to detect loci which affect many transcript levels and transcriptional pathways (Jansen and Nap 2001). We found that the let-60(gf) mutation mainly reveals novel trans-eQTL, which are eQTL distant from the regulated gene. In contrast, cis-eQTL, genes that colocalize with the eQTL, were mostly independent of the let-60(gf) mutation. Of the trans-eQTL, 77 genes have an eQTL in a trans-band (eQTL hotspot) colocalizing with amx-2. Comparing the transcriptional profiles of transgenic lines with the N2 amx-2 allele with the CB4856 amx-2 allele showed the involvement of amx-2 in 79% (61/77) of the genes mapping to this locus. Through networkassisted gene expression analysis, we found evidence that amx-2 indirectly affects gene expression mapping to this trans-band.
Strains were maintained on NGM agar seeded with OP50 bacteria at 20°, strains were age synchronized by bleaching and were harvested 72 hr postsynchronization. One sample consisted of multiple pooled, nonstarved populations from approximately four separate 9-cm NGM plates.
RNA isolation, cDNA synthesis, and cRNA synthesis The RNA of the samples was isolated using the RNEasy Micro Kit from Qiagen (Hilden, Germany), following the "Purification of Total RNA from Animal and Human Tissues" protocol with a modified lysing procedure. As prescribed, 150 ml RLT buffer and 295 ml RNAse-free water were used to lyse the samples, but with an addition of 800 mg/ml proteinase K and 1% b-mercaptoethanol. This suspension was incubated for 30 min at 55°and 1000 rpm in a Thermomixer (Eppendorf, Hamburg, Germany). Thereafter the protocol as supplied by the manufacturer was followed.
For gene expression analysis, 200 ng of RNA (as quantified by NanoDrop) was used in the "Two-Color Microarray-Based Gene Expression Analysis; Low Input Quick Amp Labeling" protocol, version 6.0, from Agilent (Agilent Technologies, Santa Clara, CA).

Microarray hybridization and scanning
The Agilent C. elegans (V2) Gene Expression Microarray 4X44K slides were used to measure gene expression. As recommended by the manufacturer, the samples were hybridized for 17 hr and scanned by an Agilent High Resolution C Scanner. The intensities were extracted using Agilent Feature Extraction Software (version 10.7.1.1). The data were normalized using the Limma package in "R" (version 3.3.0 x64). For within-array normalization, the "Loess" method was used and the "Quantile" method for between-array normalization (Zahurak et al. 2007;Smyth and Speed 2003). Unless mentioned otherwise, the log 2 transformed intensities were used in subsequent analysis.
Expanding the genetic map Previously, the strains were genotyped with 73 fragment length polymorphism (FLP) markers (as described in Schmid et al. 2015 andZipperlen et al. 2005) (Table S2).
To increase the resolution of the genetic map, cis-eQTL of a large C. elegans eQTL experiment were used to further pinpoint the crossovers (Rockman et al. 2010). The gene expression of the miRILs was transformed to the mean of the two parental lines (R) by where Y stands for the untransformed intensities of strain x and spot i (1, 2, 3, ..., 45220). The obtained values were correlated to the cis-eQTL effects from Rockman et al. (2010). Since the current study uses a newer version of the C. elegans microarray, only the spots that occurred in both designs were compared. The expression markers were generated following a procedure that is a variation of the method described in West et al. (2006). Per 20 consecutive cis-eQTL, the correlation between the N2-eQTL effect and the transformed intensities (R) were calculated. In this way, the correlation would be negative if the genotype stemmed from CB4856, and positive if the genotype stemmed from N2.
In this way, 424 gene expression markers were generated per strain. To control for quality, the markers were filtered for calling the correct genotype in the four MT2124 replicas, the four CB4856 replicas, and in .50% of the samples an absolute correlation .0.5 was required. This pruning resulted in 204 reliable gene expression markers. The quality of these markers was controlled by predicting the genotype by expression in RIL AH2244, in which the gene expression was measured twice. There, we found that 10/204 selected gene expression markers were in disagreement with the genotype. This could be reduced to 0 by comparing only markers with an absolute correlation .0.6. Therefore, all the correlations with an absolute value .0.6 were assigned the predicted marker, which corresponds to an error rate of ,0.01 per strain. The genotypes of the remaining markers were manually inferred from the surrounding markers. Furthermore, the genotypes at the ends of the chromosomes were inferred from the distal most-assigned marker. This brings the size of the resulting map to 289 markers. The correlations and assigned genotypes can be found in Table S3.

Evaluating the genetic map
The 289 marker set was analyzed by correlation analysis for markers describing unique crossover events and to see if there are any strong linkages between chromosomes ( Figure S1). This led to the selection of 247 markers indicating the border of a crossover event. To reduce the chances of false positives, we only used markers where the frequency of the least-occurring genotype in the population was .15%. It has to be noted that this excluded most of chromosome IV [as the genotype was predominantly N2 due to selection for strains including the let-60(gf) mutation (Schmid et al. 2015)]. Furthermore, also the peel-1/zeel-1 region on chromosome I was excluded for the same reason (Seidel et al. 2008).
eQTL mapping and threshold determination Mapping of eQTL mapping was done in "R" (version 3.3.0 x64) and the gene expression data were fitted to the linear model where y is the log 2 -normalized intensity as measured by microarray of spot i (i = 1, 2, ..., 45220) of miRIL j. This is explained by the genotype (either CB4856 or N2) on marker location x (x = 1, 2, ..., 247) of miRIL j.
A permutation approach was used to determine an empirical false discovery rate (Snoek et al. 2017;Vinuela et al. 2010). First, the log 2normalized intensities were randomly distributed per gene over the genotypes. This randomized data set was tested using the eQTL mapping model. This procedure was repeated 10 times. The threshold was determined using where, at a specific significance level, the false discoveries (FDS) were the averaged outcome of the permutations and read discoveries (RDS) were the outcome of the eQTL mapping. The value of m 0, the number of true null hypotheses tested, was 45220 RDS; and for the value of m, the number of hypotheses tested, the number of spots (45220) was taken. Because this study only used a limited set of strains, a more lenient threshold of 0.1 was taken as the q-value (Benjamini and Yekutieli 2001), which resulted in a threshold of 2log 10 (p) . 3.2.

Statistical power calculations
The statistical power of the mapping was determined at the significance threshold of 2log 10 (p) . 3.2. Using the genetic map (33 strains with 247 markers), 10 QTL were simulated per marker location, explaining 20-80% of the variation (in increments of 5%). Next to a peak, random variation was also introduced, based on a standard normal distribution (m = 0 and s = 1). Peaks with corresponding explanatory power were simulated in this random variation (e.g., a peak size of 1 corresponds to 20% explained variation). Furthermore, random data without peaks was also generated. These simulated data were mapped as described above and for each simulated peak it was determined: (i) if it was detected correctly, (ii) how precise the effect size estimation was, and (iii) how precise the location determined was. An overview of the results can be found in Table S4.

eQTL analysis
The distinction between cis-and trans-eQTL was made on the distance between the gene and the eQTL peak. Given the relatively low number of unique recombinations (due to the limited set of strains), the cis-eQTL window was set at 2 Mb. This means that if a gene lies within 2 Mb of the QTL peak or the confidence interval, it is called a cis-eQTL. The confidence interval was based on a 2log 10 (p) drop of 1.5 compared to the peak. Trans-bands were identified based on a Poisson distribution of the mapped trans-eQTL (as in Rockman et al. 2010). The number of trans-eQTL was counted per 1-Mb bin, leading to the identification of 52 bins with trans-eQTL. Since we mapped 1149 trans-eQTL peaks (spots) to 52 bins, we expected 22.10 trans-eQTL per bin with trans-eQTL assigned to it. Based on a Poisson distribution, any bin linking .30 spots with a trans-eQTL would have a significance of P , 0.05, and .38 trans-eQTL would yield a significance of P , 0.001. In this way, six trans-band loci were identified.
Comparative analysis against previous eQTL experiments was done on remapped experiments downloaded from WormQTL van der Velde et al. 2014) and EleQTL (http:// www.bioinformatics.nl/EleQTL). This data set consists of four different experiments representing nine different conditions. The first set contains eQTL in two temperature conditions, 16°and 24°, measured in the L3 stage (Li et al. 2006). The second set contains eQTL over three life stages: L4 juvenile animals grown at 24°, reproducing adult animals (96 hr) grown at 24°, and aging animals (214 hr) grown at 24° ( Vinuela et al. 2010). The third set contains eQTL from a single experimental condition (young adults grown at 20°) measured on a large RIL panel (Rockman et al. 2010). The fourth set contains eQTL from three experimental conditions over the course of a heat-shock treatment: a control condition (L4 animals grown for 48 hr at 20°), a heat-shock condition (L4 animals grown for 46 hr at 20°and exposed to 2 hr of 35°), and a recovery condition (similar to heat shock, only followed by 2 hr at 20°) (Snoek et al. 2017). Each of these experiments was compared at FDR = 0.05 to the eQTL mapped at FDR = 0.10 in the let-60(gf)-sensitized miRILs.
Since the data set of Rockman et al. (2010) was used for expanding the genetic map, analysis was also conducted excluding this study. The reason for excluding this set from analysis was that a bias could be introduced in overlap with the cis-eQTL. The main conclusion that cis-eQTL are less unique than trans-eQTL also stands with this analysis. Excluding the Rockman et al. data, resulted in detection of 999/1516 (65.9%) cis-eQTL and 171/898 (19.0%) trans-eQTL present in previous experiments.

Enrichment analysis
Gene group enrichment analysis was done using a hypergeometric test on the unique genes (not on the spots) with the following criteria: Bonferroni-corrected P , 0.05; size of the category, n . 3; and size of the overlap, n . 2.
amx-2 allelic comparison Four independent transgenic lines, two of amx-2(lf);Si[amx-2(Bristol)]; let-60(gf) and two of amx-2(lf);Si[amx-2(Hawaii)];let-60(gf), were used to investigate the effect of the amx-2 CB4856 and N2 allele on gene expression. These lines contain single-copy insertions on LGII of the N2 or CB4856 alleles of amx-2 in an N2 amx-2(lf); let-60(gf) background. The effect of different alleles on gene expression was measured by microarray for each independent transgenic line and compared to the effects of the eQTL mapping to the trans-band closest to the position of amx-2.

Connectivity network analysis
To investigate the connectivity and function of the affected genes we used WormNet (version 3) (Cho et al. 2014) and GeneMANIA plug-in for Cytoscape (version 3.4.0) (Montojo et al. 2010;Shannon et al. 2003). WormNet was used to investigate enrichment in connectivity within groups of genes. For example, groups of coexpressed genes mapping to the same trans-band were assumed to share the same regulator. Further evidence for this coexpression and regulation can be found in WormNet if these genes are significantly more connected than by chance. GeneMANIA was used to find the closest neighbors, those genes with the most connections, of amx-2 and let-60 to locate the genes by which amx-2 modifies RAS signaling.

Data availability
All strains used can be requested from the authors. The transcriptome data sets generated and the mapped eQTL profiles can be interactively accessed via http://www.bioinformatics.nl/EleQTL. Moreover, the raw transcriptomic data are also available at ArrayExpress (E-MTAB-5856).

RESULTS eQTL mapping in miRIL population
Using microarrays, we measured gene expression levels in 33 miRILs sensitized for RAS/MAPK signaling by introgression of a let-60 (gf) mutation in a segregated N2/CB4856 genetic background ( Figure S2) (Schmid et al. 2015). Analysis of the statistical power of this population with a genetic map of 247 informative markers showed that we can detect 77% of the QTL explaining 40% of the variation (Table S4). We detected 2303 genes (represented by 3226 array spots) with at least one eQTL (FDR = 0.1, 2log 10 (p) . 3.2; Figure 1, Table 1, and Table S5). For 1516 of these genes, a cis-eQTL was found, indicating regulation from within a window of 2 Mb around the affected gene. For 898 genes a trans-eQTL was found, showing regulation more distal from the affected transcript. Most cis-eQTL (1074 out of 1516; $71%) show a positive effect for the N2 allele; whereas for the trans-eQTL, the positive and negative effects were equally distributed over the N2 and CB4856 allele ( Table 1).
The eQTL were not equally distributed across the genome. We detected different clustered trans-eQTL as "hotspots" or trans-bands. Trans-bands are frequently found in eQTL experiments and indicate a locus with one or multiple alleles affecting the expression of many distant genes. To identify the trans-bands in our experiment, we calculated the chance of overrepresentation of trans-eQTL per marker, using a Poisson distribution (as in Rockman et al. 2010). This method was applied per 1-Mb window on the genome. The trans-bands that were identified in adjacent windows were merged. This way, six transbands were identified (Poisson distribution, P , 0.001; Table 2).
Specificity of let-60 eQTL By comparative analysis of different environments and/or populations, it has been shown that eQTL trans-bands can be population, environment, or age specific (for example, see Li et al. 2006). To investigate if and which eQTL are specifically found in the sensitized miRILs, we compared the detected eQTL to the eQTL found in nonsensitized RILs (Li et al. 2006;Rockman et al. 2010;Vinuela et al. 2010;Snoek et al. 2017). These data sets contain eQTL mapped over nine different conditions (see Materials and Methods), which were compared at an FDR = 0.05 with the eQTL mapped in the sensitized miRILs. We found that $73% (1112 out of 1516) of the genes with a cis-eQTL in the miRIL population were also found in one or more of the previous studies (Figure 2 and Figure S3B). When taking the allelic effects into account, of the 1074 genes with a cis-eQTL that gave a higher expression linked to the N2 allele, 816 ($76%) were found in other studies. For the 456 genes with a cis-eQTL that gave a higher expression linked to the CB4856 allele, 310 ($68%) were found in other studies. This shows that the majority of the cis-eQTL detected in the miRIL population can be also be detected in other experiments.
In contrast, for the trans-eQTL detected in our experiment, not taking location into account, we only found $24% (214 out of 898) of the genes had a trans-eQTL in at least one of the previous experiments (Figure 2   (2196) 898 (1149) The numbers in brackets indicate the number of spots. a The discrepancy between the sum of N2 higher, CB4856 higher, and the total is due to genes being represented by multiple spots, which often represent different splice variants.
and Figure S3C). To further compare the trans-eQTL overlap, it was investigated whether the trans-eQTL colocalized on the same chromosome. Using that criterion, $9% (80 out of 898) of the genes with a trans-eQTL were found in a previous study. This showed that the majority of the genes with a trans-eQTL were specifically detected in the sensitized miRIL population. This observation also extended to the trans-bands. To test the specificity of the trans-bands, we counted the number of times a gene within the miRIL trans-bands had a colocating trans-eQTL in other studies. All trans-bands were found to be specific for the let-60(gf) miRILs (Poisson distribution, P , 0.001; Table 2). This provides additional evidence that these trans-bands are specific for interaction between the genetic background and the let-60(gf) mutation. Altogether, we identified 404 genes with a novel cis-eQTL and 684 genes with a novel trans-eQTL ( Figure 3A). As our method can be biased by the way the genetic map was constructed (using expression differences, see Materials and Methods), we conducted the same analyses using the original FLP map. Our conclusions on the detection of novel cis-and trans-eQTL were not affected by the different genetic map (for details, see File S1). This implies that a substantial majority of trans-eQTL showed a mutation-dependent induction. This is also illustrated by showing the eQTL that were consistently found ( Figure  3B), which mainly shows the cis-diagonal-and a few trans-eQTL.
Functional allelic differences Enrichment analysis on groups of genes with an eQTL mapping to the same trans-band can be used to uncover the functional consequences of the allelic variation at the selected locus. Enrichment analysis was done on GO-terms, KEGG, anatomy terms, disease phenotypes, development expression, phenotypes, transcription-factor binding sites, and RNAi phenotypes (Table S6).
The genes with cis-eQTL were enriched for genes in the gene classes math, bath, btb, and fbx. These classes contain many genes that function as intra-and extracellular sensors involved in gene-environment interactions, and were found highly polymorphic between N2 and CB4856, but also between other wild isolates (Thompson et al. 2015;Volkers et al. 2013). The genes with trans-eQTL with higher expression linked to CB4856 loci were enriched for genes related to development in general and gonad development specifically. These enrichments were likely to stem from the trans-bands on chromosome I: 12-15 Mb and II: 3-6 Mb specifically. These two transbands mainly consisted of genes with a higher expression linked to the CB4856 locus and were enriched for genes associated with reproduction and/or development.

Allelic effect of amx-2 on gene expression
The trans-band on chromosome I: 12-15 Mb colocalizes with a previously identified QTL for vulval induction, for which amx-2 was identified as polymorphic regulator (Schmid et al. 2015). To determine if the trans-eQTL mapping to the amx-2 locus were affected by the allelic difference between N2 and CB4856, we measured gene expression in transgenic lines containing single-copy insertions of either allele of amx-2 in a amx-2(lf); let-60(gf) background. For each allele, two independent strains were created and the transcriptome was measured using microarrays (Schmid et al. 2015). From these measurements, the allelic effects of amx-2 on gene expression were determined and compared with the trans-eQTL effects measured in the miRIL population. These were found to correlate significantly ( Figure S4; R 2 $0.34; P , 2 · 10 28 ). Specifically, 61 of the 77 genes with an eQTL colocalizing with amx-2 showed the same allelic effect in the transgenic strains carrying the two amx-2 variants. Furthermore, analysis of the effect directions showed that the trans-band on chromosome I (Table 2) consists of two separate trans-bands (see Figure S4); the first colocalizing with amx-2 and the second lying more distal on chromosome I.
To further investigate the mechanism by which amx-2 modifies let-60 Ras/MAPK signaling, we used the well-established C. elegans connectivity networks WormNet (Cho et al. 2014) and GeneMANIA (Montojo et al. 2010;Shannon et al. 2003). WormNet showed that the genes mapping to the amx-2 locus were more connected than expected by chance (P , 3 · 10 213 ). Visualizing the connectivity of these genes together with amx-2 and let-60 using GeneMANIA identified one major gene cluster ( Figure 4). As expected, genes affected by the transband were found to be highly connected and at the core of the cluster, indicating a shared biological function. In contrast, 25% of the genes with a cis-eQTL at the trans-band location were not connected, and those that were connected had only one to three connections. Interestingly, let-60 was at the core of the cluster of genes with a trans-eQTL, whereas amx-2 was at the periphery. We found six genes, lin-40 (also called egr-1), egl-27, trr-1, pbrm-1, ceh-26 (also called pros-1), and emb-5, that were previously found to have a genetic interactions n    (Li et al. 2006;Rockman et al. 2010;Vinuela et al. 2010;Snoek et al. 2017). The histogram shows how many of the eQTL mapped in this study were found back in these independent conditions. Whereas the majority of cis-eQTL was mapped in previous studies, there was only a minority of trans-eQTL that had been mapped.

DISCUSSION
Ras/MAPK signaling is strongly affected by gain-of-function mutations in let-60, as shown by the strong effect on vulva development in C. elegans (Beitel et al. 1990). The allele let-60(n1046) hyper-activates Ras/MAPK signaling, leading to a multivulva phenotype and the differentiation of more than three vulval precursor cells. In the miRIL population, which carries a let-60(n1046) mutation in a segregating N2/CB4856 background, a wide range of VI was found. While in the full N2 background, VI of the let-60(gf) mutation was 3.7 induced vulval cells on average, the VI in the miRIL population varied from 3.0 to 5.7 induced vulval cells, illustrating the strong modulatory effect of the genetic background on the mutant phenotype (Schmid et al. 2015). By measuring gene expression across the different miRILs, we gained insight into the genetic background effects on the transcriptional architecture of let-60(gf)-sensitized miRILs.
The let-60(gf) mutation affects gene expression in trans To our knowledge, this is the first study where a mutation, introgressed into a panel of RILs, is used to investigate the genetic architecture of transcript variation by eQTL analysis. Unfortunately, the direct effect of the let-60(gf) allele on the genetic background is impossible to determine in this population, since there is no population with the same genetic background lacking the let-60(gf) mutation. Therefore, we used previously published eQTL studies in C. elegans as an estimation for the miRIL-specific eQTL. We found that introgression of let-60(gf) leads to an overrepresentation of specific trans-eQTL, while the cis-eQTL showed a high overlap with eQTL observed in the absence of let-60(gf) allele. Although our study is based on a relatively small number of strains, the verification of identified eQTL against other eQTL studies in C. elegans shows that at least 73% of the identified cis-eQTL were previously discovered. This supports our power analysis, showing that our study has enough statistical power to detect eQTL with larger effects. Given that cis-eQTL are most likely genes with polymorphisms in or near the gene itself, it is expected that most have been previously discovered. In support of this, most of the cis-eQTL genes show higher expression in N2, making it likely that a relatively large proportion of these cis-eQTL indeed stem from polymorphisms causing mishybridization on the microarrays instead of true expression differences (Rockman et al. 2010;Alberts et al. 2007). Therefore, most of the cis-eQTL will not represent gene expression changes that can be linked to the let-60(gf) mutation in the miRIL population. That is also exactly what we found when constructing the regulatory network by connecting cis-and trans-eQTL to let-60; cis-eQTL are only loosely, or not at all, connected to the network.
Compared to cis-eQTL, the trans-eQTL identified in the miRIL population were hardly found in previous eQTL studies in C. elegans (only 24% were previously detected). However, it should be noted that trans-eQTL are themselves are also environment and age dependent (Vinuela et al. 2010;Francesconi and Lehner 2014;Li et al. 2006;Snoek et al. 2017), therefore it is possible that not all the trans-eQTL are indeed let-60 induced. For example, they might be age related (Vinuela et al. 2010). Another aspect that influences detection of trans-eQTL is the power of the study, as trans-eQTL explain less of the variation than cis-eQTL. However, by being stringent in the overlap (requiring only one match over nine different environments), it is likely that most of these eQTL would be detected. This makes it likely that many of the identified trans-eQTL are specifically linked to an effect of the let-60(gf) mutation in the genetic background. This extends the idea that environmental Figure 3 eQTL mapped in the let-60(gf) miRIL population compared to other studies (Li et al. 2006;Rockman et al. 2010;Vinuela et al. 2010;Snoek et al. 2017). As in Figure  perturbations can reveal additional genetic variation (Li et al. 2008) by perturbing the genetic environment (Kammenga et al. 2008;Duveau and Felix 2012;Schmid et al. 2015). The advantage of mutational perturbation and perturbation via induced responses in eQTL studies (for example, see Snoek et al. 2012 andSnoek et al. 2017) is that it contextualizes the transcriptional response. Yet, as the phenotypic influence of let-60(gf) in the miRIL population was investigated before, we know it results in a phenotype. In the study presented here, the context is the let-60(gf) mutation, which results in a phenotype where VI is affected. It is therefore interesting that the trans-bands colocalize with QTL mapped for VI in the same miRIL population: trans-band I: 12-15 Mb overlaps with QTL1b, trans-band II: 7-8 Mb with QTL2, and trans-band V: 13-14 Mb with QTL3 (Schmid et al. 2015). This colocalization adds to the plausibility that the novel trans-eQTL detected in the miRIL population are indeed due to the let-60(gf) mutation and its interaction with the genetic background.
A trans-band on chromosome I links to variation in the amx-2 gene Using the same let-60(gf) miRIL population, we have previously identified the polymorphic MAOA amx-2 as a background modifier that negatively regulates the Ras/MAPK pathway and the VI phenotype (Schmid et al. 2015). The trans-band on chromosome I colocalizes with the amx-2 QTL for VI [QTL1b (Schmid et al. 2015)], it is therefore possible that the modifier amx-2 also affects gene expression. More specifically, the trans-band on chromosome I originally identified in the miRILs consists of two separate regulatory loci. The N2 allelic effect at the amx-2 locus is mostly negative; whereas for the more distal part of the original trans-band, the N2 allelic effect is mostly positive. Yet, it seems unlikely that amx-2 is a direct gene expression effector. As amx-2 is a mitochondrial MAOA, a catalyzer of neuropeptide oxidative deamination, it will probably not influence gene expression directly (Tipton et al. 2004). Placing the trans-eQTL of the amx-2 trans-band in a gene interaction network supports this line of reasoning: although let-60 is in the center of the trans-eQTL, amx-2 is only peripheral.
By measuring gene expression in transgenic lines expressing the N2 or CB4856 allelic variants of amx-2 in an amx-2(lf);let-60(gf) genetic background, we attempted to link the allelic effect of amx-2 to the trans-eQTL mapping to the amx-2 locus. As a significant correlation between the expression differences in the transgenic lines and the trans-eQTL effects mapping to the amx-2 locus is found, these trans-eQTL can be confidently linked to allelic variation in amx-2. As discussed Figure 4 Interaction network of genes colocalizing with amx-2. Genes with an eQTL colocalizing with amx-2 were selected together with amx-2 (yellow) and let-60 (pink). Interactions were obtained from GeneMANIA (Montojo et al. 2010;Shannon et al. 2003). Node color indicates cis in black and trans in green. Node shape indicates eQTL effect, N2 . CB4856 as a s and CB4856 . N2 as a h. Node border color gradient indicates expression in amx-2 transgenic lines (blue N2-amx-2-allele . CB4856-amx-2 allele, red CB4856-amx-2 allele . N2 amx-2 allele). Node size indicates number of edges connected. Edges in blue show coexpression links and edges in pink show genetic interactions as found by Lee et al. (2010), Lehner et al. (2006), and Byrne et al. (2007). Genes connected by genetic interactions have boldface gene names (as well as amx-2). before, it is unlikely that amx-2 is the direct cause of the transcriptional variation, but rather it acts through an indirect mechanism. One route of effect might be through amx-2-mediated degradation of serotonin (5-HT) to 5-HIAA, which both affect VI (as discussed in Schmid et al. 2015). Subsequently, the affected VI might result in different gene expression levels. It is interesting to remark that the amx-2 trans-band is characterized by downregulation of expression related to the N2 genotype as well as a decreased VI for that genotype.
As this places amx-2 in the causal chain of events, we think it is most likely that amx-2 does not affect gene expression directly but via another gene involved in Ras/MAPK signaling, possibly directly linked to let-60. Therefore, we hypothesize that amx-2 affects Ras/MAPK signaling via one or more of the six previously found genes with a genetic interaction with let-60: lin-40 (also called egr-1), egl-27, trr-1, pbrm-1, ceh-26 (also called pros-1), and emb-5. For example, egl-27 and lin-40 both have a trans-eQTL mapping to the amx-2 locus, and are the two MTA (metastasis-associated protein) homologs found in C. elegans Herman et al. 1999;Solari et al. 1999). The proteins act in the NURD chromatin remodeling complex, which has previously been shown to antagonize Ras-induced vulval development (Solari and Ahringer 2000). Chromatin remodeling provides a more likely mechanism of action than the molecular role of amx-2 itself.
Implications for understanding the Ras/MAPK pathway How can our results help form a better understanding of the Ras/MAPK pathway? Our previous investigation on VI in the miRIL population resulted in the identification of three QTL harboring polymorphic Ras-signaling modifiers. Expanding our research to genetic variation affecting gene expression in a let-60(gf)-sensitized RIL population uncovered six trans-bands (or "eQTL hotspots"). As mentioned before, these trans-bands overlap with the QTL mapped for VI in Schmid et al. (2015). There are three more trans-bands that do not overlap with QTL for VI, but are also likely to represent modifier loci of the Ras/MAPK pathway, possibly underlying other Ras/MARK-associated phenotypic differences.
The main advantage of studying the genetic architecture of gene expression in the miRIL population is that it creates more insight in the genes and pathways affected by allelic variation acting on the Ras/MAPK pathway. It identifies hidden genetic variation; genetic variants that are unlocked under altered environmental conditions or when the genetic background is modified (for a review, see Rockman 2014 andKammenga 2017). Identification of background modifiers of disease pathways is important for gaining insight into individual-based differences of disease contraction. Mutant phenotypes can be strongly affected by the genetic background. For example, large variation in traits between different backgrounds in many different mutated genes has been observed in C. elegans (Paaby et al. 2015;Vu et al. 2015;Duveau and Felix 2012). A major discovery was that this variation in mutant phenotypes could be predicted from gene expression variation (Vu et al. 2015). These results are in line with the discovery of trans-bands at the location of each VI QTL. Furthermore, results presented here (and previously in Elvin et al. 2011) link variations in individual gene expression levels to enhancing or diminishing the severity of a Mendelian disorder caused by a so-called "major gene." More specifically, these differences seem to stem from a couple of modifier loci harboring polymorphic regulators.
So far, we only explored the amx-2 trans-band, and it is likely that the other loci also contain polymorphic modifiers of the Ras/MAPK pathway. Genetic perturbation by let-60(gf) leads to a strong increase in specific trans-acting eQTL organized over six trans-bands, thus supporting the involvement of multiple genes as modifiers. Given the fact that we detected amx-2 both via an indirect transcriptional response, but also mechanistically, we feel confident that the detected loci indeed harbor other modifiers. A single gene mutation apparently has the capacity to unlock a huge number of novel interactions controlled by many genes across different genetic backgrounds. Further studies should be conducted to identify and characterize these underlying polymorphic regulators.

Conclusion
Introgression of a mutation in a segregating genetic background allows for identification of polymorphic modifier loci. By introducing a let-60(gf) mutation in a N2/CB4856 genetic background, in the form of the miRIL population, and measuring the transcriptome, we identified six trans-bands specific for the let-60(gf) miRIL population. The majority of trans-eQTL are specific for the miRIL population, showing that genetic variation in gene expression can be specifically modified by a background mutation. Therefore, genetic perturbation can be viewed as analogous to environmental perturbation, which also results in specific trans-eQTL. We demonstrated the involvement of amx-2 and the allelic variation between N2 and CB4856 in gene expression variation originating from chromosome I. Yet, we think it is unlikely that amx-2 directly affects gene expression variation. Instead, we prefer the hypothesis that allelic variation in amx-2 indirectly affects gene expression, possibly through the NURD complex.