Modeling of a negative feedback mechanism explains age-dependent genetic architecture in reproduction in domesticated C. elegans strains

Most biological traits and common diseases have a strong but complex genetic basis, controlled by large numbers of genetic variants with small contributions to a trait or disease risk. The effect-size of most genetic variants is not absolute, but can depend on a number of factors including the age and genetic background of an organism. In order to understand the mechanisms that cause these changes, we are studying heritable trait differences between two domesticated strains of C. elegans. We previously identified a major effect locus, caused by a mutation in a component of the NURF chromatin remodeling complex, that regulated reproductive output in an age-dependent manner. The effect-size of this locus changes from positive to negative over the course of an animal’s reproductive lifespan. Using a previously published macroscale model of egg-laying rate in C. elegans, we show how time-dependent effect-size can be explained by an unequal use of sperm combined with negative feedback between sperm and ovulation rate. We validate a number of key predictions of this model using controlled mating experiments and quantification of oogenesis and sperm use. By incorporating this model into QTL mapping, we identify and partition new QTLs into specific aspects of the egg-laying process. Finally, we show how epistasis between two genetic variants is predicted by this modeling as a consequence of unequal use of sperm. This work demonstrates how modeling of multicellular communication systems can improve our ability to predict and understand the role of genetic variation on a complex phenotype. Negative autoregulatory feedback loops, common in transcriptional regulation, could play an important role in modifying genetic architecture in other traits. AUTHOR SUMMARY Complex traits are influenced not only by the individual effects of genetic variants, but also how these variants interact with the environment, age, and each other. While complex genetic architectures seem to be ubiquitous in natural traits, little is known about the mechanisms that cause them. Here we identify an example of age-dependent genetic architecture controlling the rate and timing of reproduction in the hermaphroditic nematode C. elegans. Using computational modeling, we demonstrate how this age-dependent genetic architecture can arise as a consequence of two factors: hormonal feedback on oocytes mediated by major sperm protein (MSP) released by sperm stored in the spermatheca and life history differences in sperm use caused by genetic variants. Our work also suggests how age-dependent epistasis can emerge from multicellular feedback systems.


28
GWAS and QTL mapping, two common quantitative genetics techniques, are usually unable to 29 identify interacting genetic variants. GWAS, which can narrow down causative genetic variants 30 to small regions, are typically underpowered to identify statistically significant epistatic 31 interactions due to low natural allele frequencies and the large number of statistical tests that 32 must be performed (23). QTL mapping, on the other hand, has increased power to identify 33 interacting QTLs due to equal allele frequencies but identifies large regions in linkage 1 disequilibrium containing thousands of potential variants (24). Due to this inherent difficulty in 2 studying epistasis that occurs between genetic variants that segregate within a population, 3 studies of epistasis have typically focused on laboratory induced loss-of-function mutations 4 through mutagenesis or RNAi. It is unknown whether the mechanisms that cause epistasis in 5 these approaches will apply in natural populations.

6
As a more tractable model to understand how genetic variants impact a trait, we are focusing 7 our studies on two C. elegans strains, N2 and LSJ2, derived from an individual hermaphrodite 8 isolated in 1951. Due to a reproductive system that primarily relies on self-fertilization, the 9 population these two strains derived from was genetically identical, at which point they were

22
To address this deficiency, we studied the genetic basis of reproductive differences between the 23 N2 and LSJ2 strains at five different time points spanning their reproductive lifespan. Our goal 24 for this study was to identify examples of complex genetic architecture and understand their 25 molecular and cellular causes.

27
A major effect QTL surrounding nurf-1 has an age-dependent effect size

28
We previously performed QTL mapping on reproductive rate using 94 recombinant inbred 29 strains (RILs) generated between LSJ2 and CX12311 (29). CX12311 is a strain that derives the 30 majority (>99%) of its DNA from N2 with the exception of a small amount of DNA backcrossed 31 from the CB4856 wild strain that surrounds the npr-1 and glb-5 genes ( Figure 1B) (26). Novel 32 mutations in these two genes became fixed in the N2 lineage and result in pleiotropic effects on 1 a large number of phenotypes. Use of the CX12311 strain allows us to avoid studying their 2 effects. To study the role age plays on reproduction, we chose five time points that spanned the 3 reproductive lifespan of the CX12311 animals. Egg-laying rate was quantified by counting the 4 number of eggs laid by six animals for six hours on agar plates seeded with E. coli bacteria 5 ( Figure 1C). We previously identified a major effect QTL centered over the nurf-1 gene that 6 accounted for ~50% of the observed phenotypic variation (26). To study how the animal's age 7 affected the effect-size of this locus, we first segregated the 94 RIL strains based upon their 8 genotype at nurf-1 (Figure 1D -top panel). At the first three time points, RIL strains with the 9 N2 genotype laid more eggs than RIL strains with the LSJ2 genotype, however, at the fourth 10 and fifth time point, this relationship flipped -animals with the LSJ2 allele of nurf-1 laid more 11 eggs than animals with the N2 allele. To visualize this effect more clearly, we also plotted the

17
To verify these observations, we next assayed a near isogenic line (NIL) constructed by 18 backcrossing the region surrounding nurf-1 from LSJ2 into the CX12311 strain ( Figure 1B) 19 along with the CX12311 and LSJ2 parental strains. The CX12311 strain (containing the N2 20 allele of the nurf-1 locus) laid more eggs than the NIL nurf-1 strain for the first three time points but 21 fewer eggs at the fourth and fifth time points, again resulting in a time-dependent effect size that 22 flips from positive to negative ( Figure 1E -top and bottom panel). The effect size of nurf-1 23 calculated using both the RIL and NIL strains and showed good qualitative agreement ( Figure   24 1D and 1E -bottom panel). Due to additional segregating variants in the RIL strains, we do 25 not expect these two calculations to be identical. In both experimental designs, the crossing of 26 the two lines was due to a decrease in the egg-laying rate in the strain containing N2 nurf-1 as 27 opposed to an increase in egg-laying of the LSJ2 nurf-1 strains. The LSJ2 strain was statistically 28 indistinguishable at the first three time points from the NIL strain. However, at the last two time 29 points, LSJ2 laid additional eggs resulting in a larger effect-size at these two points. These 30 results demonstrate that the nurf-1 locus can have both positive and negative effects on egg-31 laying in a manner correlated with the animal's age.

32
Difference in early egg-laying rate caused by differences in germline stem cell 1 production, oocyte maturation and/or rate of fertilization 2 The rate that eggs are laid on an agar plate is dependent on a large number of factors: size and 3 rate of mitosis of a germline progenitor pool, speed of meiosis/differentiation of these cells, 4 maturation and growth of oocytes, ovulation and fertilization of the primary oocyte to produce an 5 egg, and finally the rate of active expulsion of an egg through a controlled motor program that 6 results in the opening of the vulva (Figure 2A). To characterize which factor might be affected 7 by nurf-1, we first characterized the number of fertilized eggs in each strain using DIC 8 microscopy. If the rate fertilized eggs were laid was affected in LSJ2 and NIL nurf-1 strains but the 9 rate of production of fertilized eggs remained the same, we would expect LSJ2 and NIL nurf-1 to 10 contain more eggs in their uterus than the CX12311 strain. We measured the number of 11 fertilized eggs in each of these three strains at 24 and 48 hours after the L4 stage ( Figure 2B).

12
In contrast to this prediction, we found that both LSJ2 and NIL nurf-1 had significantly less fertilized 13 eggs than the CX12311 strain. We conclude that production of fertilized eggs must be affected 14 in the LSJ2 and NIL nurf-1 strain. The lower number of unlaid fertilized eggs in these strains is 15 potentially a consequence of the reduced rate of production of fertilized eggs.

16
We next measured the number of large oocytes undergoing oogenesis in these three strains 17 (region marked oogenesis in Figure 2A). We hoped to distinguish between two possible 18 mechanisms modifying the rate of production of fertilized eggs: 1) the rate of production of 19 mature oocytes was decreased in LSJ2/NIL nurf-1 strains or 2) the rate of fertilization of a mature 20 oocyte was decreased in the LSJ2/NIL nurf-1 strains. We reasoned that in the former case, the 21 number of large oocytes undergoing oogenesis would be lower in LSJ2 and NIL nurf-1 due to the 22 decrease in production. In the latter case, if oocyte production was unaffected but the rate they 23 were fertilized was lowered, we reasoned that the number of mature, large oocytes would 24 increase over time. However, when we measured the number of large oocytes in each strain at 25 24 and 48 hours, we found no statistically significant difference between any of the three strains 26 ( Figure 2C). We believe this likely indicates the presence of a homeostatic mechanism that 27 keeps the number of oocytes undergoing maturation constant, despite a difference in how often 28 they are created or fertilized.

29
Finally, we measured the rate that progenitor germline stem produce new cells through mitosis, 30 as these cells provide the source of meiotic cells that later differentiate into oocytes (marked as 31 mitotic zone in Figure 2A). There is not a one-to-one relationship between the rate of mitosis 32 and oocyte production due to cannibalization of a subset of these cells, but it is thought that 33 there is a relationship between the size of the progenitor pool and subsequent rate of 1 reproduction (31). Progenitor cells can be distinguished from cells in the transition zone based 2 upon nuclear morphology. We observed a significant difference between CX12311 and NIL nurf-1 3 at the 24 hour time point in the number of progenitor cells, suggesting that the difference in egg-4 laying rate could be caused by this difference (Figure 2D). However, the NIL nurf-1 strain was also    also will decrease (35, 36). This suggests a possible mechanism tying an animal's age to its 33 egg-laying rate. As an animal ages, each egg it creates through fertilization reduces the number 1 of sperm by one. Animals with N2 nurf-1 lay more eggs during the first three time points 2 compared to LSJ2 nurf-1 animals and will consequently have fewer sperm at the fourth and fifth 3 time point. This reduction in sperm will lower sperm hormone present at the primary oocyte, 4 which counterbalances the positive effect of the N2 nurf-1 allele on egg-laying during the 5 transition from between the third and fourth timepoint.

6
We tested this hypothesis more formally using a previously published macroscale model (37), 7 which stipulates the egg-laying rate is proportional to the product of the number of oocytes with 8 the number of sperm ( Figure 3B). This model is defined by four time-independent parameters: 9 k o , which specifies how rapidly oocytes are created, k c , which specifies a carrying capacity of 10 the gonad, k f , which defines the fertilization rate, and S 0 , which specifies the number of self-11 sperm created by the animal. In this report, we excluded the k c parameter, as we found that the . This is in good qualitative agreement with our previous demonstration 23 that LSJ2 animals lay more eggs over the course of their lifetime than CX12311 (29). However, 24 the quantitative agreement was rather poor for the LSJ2 strain. While the predicted S 0 value of 25 CX12311 was in good agreement with our previously measured fecundity (263 vs. 256), the 26 predicted value of S 0 was significantly higher than the measured fecundity (434 vs. 319). This is 27 most likely due to the inability of the model to fully account for the genetic changes in the LSJ2 28 strain. The average residuals for the LSJ2 strain were significantly higher than CX12311 (2.0 vs. This analysis showed that a constant change to oocyte generation rate could have a time-1 dependent effect on egg-laying resulting in sign-switching at later time points. Our model 2 explicitly predicts that the reduction in CX12311 animals in egg-laying rate at later time points is 3 not due to changes in reproductive capacity but rather decreases in sperm number. To test this 4 prediction, we first measured the number of fertilized eggs and large oocytes in CX12311 5 animals at 66 hours ( Figure 4A). In line with our expectations, we observed a reduction in the 6 number of fertilized eggs (~10) contained in the CX12311 animals compared to earlier time 7 points ( Figure 2B). This indicates that the decrease in egg-laying rate cannot be explained by 8 reduction of in the rate of laying fertilized eggs in the uterus. We also observed a large number 9 of large oocytes (~22) in CX12311 animals, an increase over the previous two time points 10 ( Figure 4A and 2C). This suggests that the change in egg-laying rate is not due to changes in 11 oocyte maturation rate, as mature oocytes are available for fertilization. These observations are 12 consistent with our hypothesis that decreased sperm and MSP at later timepoints are 13 responsible for the decreased rates of egg-laying in CX12311.
14 Our modeling predicts that the CX12311 strain will have fewer sperm later on in life due to their

20
As a final method to test this model, we took advantage of C. elegans' androdioecious mating 21 system to increase the sperm available to the hermaphrodites at later time points by mating 22 CX12311 and NIL nurf-1 hermaphroties to CX12311 males, which result in the transfer of ~1000 23 sperm to the hermaphrotite. After mating young adult animals, we separated the 24 hermaphrodites from males and measured the egg-laying rate at two time points late in life.

25
Consistent with our predictions, increasing sperm number prevented the sign-switching from 26 occurring at both of these times in unmated animals ( Figure 4D). This indicates that decreases 27 in sperm number are the primary reason for decreases in egg-laying rate at late time points.

28
Identification of modifier QTLs that also effect the egg-laying process

29
We next investigated whether additional QTLs affected egg-laying rate differences between 30 LSJ2 and CX12311. The presence of a major effect QTL (such as the nurf-1 locus) can mask 31 the effects of smaller QTLs so we performed additional scans using the genotype of nurf-1 as 32 an additive or interacting covariate ( Figure 5A) (38). We identified five QTLs that were significant at the genome-wide level: one QTL on the center of chromosome I, one QTL on the 1 center of II, one large QTL on the center and right arm of chromosome IV, one QTL on the right 2 arm of V, and one QTL in the center of X (Figure S1 -S5). The identification of the QTL on 3 chromosome I was expected, as it contains a previously described missense mutation in nath-4 10 that was shown to affect egg-laying (30). However, the other four QTLs do not contain any 5 genetic variants that are known to affect egg-laying.

6
In order to analyze how these additional five modifier QTLs regulated the egg-laying rate at the  slope is positive and the other slope is negative). This analysis indicated that the effect of these 21 modifier QTLs were also complex. Each QTL was age dependent and showed non-linear 22 interactions with the nurf-1 locus. This mapping suggests that the egg-laying differences 23 between N2 and LSJ2 are multigenic, involve extensive epistatic interactions, and are highly 24 age-dependent.

25
Modeling also predicts age-dependent epistasis at latter timepoints 26 Further inspection of the age-dependence and epistasis of the modifier QTLs with nurf-1 27 revealed a number of interesting trends (Figure 5B). Epistasis was most likely to occur at the 28 first and last two time points (7 out of 8 significant effects) and less likely to be observed at the Finally, we decided to test whether we could use the macroscale modelling of the egg-laying 19 process to improve our QTL mapping. Instead of performing QTL mapping on each of the five 20 time points, we used this data to estimate a k o and S 0 for each RIL strain (assuming as before 21 that k c = 0 and k f was the same for all RIL strains). The distribution of these fitted parameters 22 are plotted in Figure 6A and 6B. We next performed QTL mapping to identify genetic regions 23 that influence these rates. As expected, we identified a number of loci that we found from 24 previous analysis. The QTL surrounding nurf-1 was identified as a regulator of both the oocyte 25 generation rate (k o ) and the number of self-sperm (S 0 ). The modifier QTLs on II and X regulated 26 the oocyte generation rate but not the number of sperm number, while the QTL on V had an 27 effect on sperm number but not the oocyte generation rate. Interestingly, this analysis also 28 identified a new QTL on chromosome III as a regulator of sperm number (Figure 6B),

29
suggesting that modeling of the egg-laying process not only helps us understand the modifier 30 QTLs effect but also can help identify new QTL loci.

31
In order to test these results, we utilized a NIL strain surrounding the QTL on V. We crossed the 32 NIL V strain to the NIL nurf-1 strain to create the double NIL nurf-1;V strain ( Figure 6D). This combination of strains allowed us to test the effect of the modifier QTL in both nurf-1 genotypes.

1
For each of these strains, we measured the total amount of progeny produced. These NIL 2 strains showed that the QTL V increased the brood size in both nurf-1 backgrounds, but had a 3 larger effect size in the LSJ2 nurf-1 genotype (60 vs. 30) indicating epistasis does exist between 4 these two loci. This qualitatively agrees with the QTL mapping data (Figure 6B -right panel) 5 but discrepancies were observed between the magnitude of the total broodsizes, suggesting 6 additional improvements to the model might be beneficial.     not segregated appropriately by age or genetic background. Our results provide a framework to understand how age-dependence can also arise through emergent properties of a cellular 1 network, which we believe to be the major scientific contribution of this work. Our work indicates 2 that the complex seesaw of effects that nurf-1 has on reproductive output can be explained 3 using two major considerations: (1) a hormonally-mediated negative feedback loop linking 4 sperm with oocyte maturation, and (2) the effect nurf-1 has on the rate of oogenesis.

5
Consequently, the molecular details of how nurf-1 modifies protein and cellular function are 6 unnecessary to explain its age-dependence. Our modeling experiments also demonstrated how 7 sign epistasis could arise in an age-dependent manner strictly through age-independent 8 changes in the oocyte production rate. The origin of genetic epistasis is often thought of in terms 9 of biochemical properties of proteins, through physical interactions between two proteins, or 10 through multiple changes to a parallel or linear signaling pathway (40-42). However, these 11 mechanisms were typically identified through analysis of mutations or genetic perturbations 12 generated in the laboratory that have a strong negative effect on fitness. It is interesting to us 13 that such mechanisms don't seem to be at play here. This is the second example we have

23
However, as the amount of protein product increases, the transcription factor will turn off the 24 expression of its own mRNA. Since the amount of protein product will be higher in the strain 25 with the higher initial expression, the transcription will turn off sooner in that strain, and the 26 genetic variant will soon appear to have a negative effect on expression.

27
Embryo/oocytes/sperm/germline stem cells were quantified by fixing whole animals in 95% 28 ethanol and staining nuclei with 1.5 mg/mL DAPI contained in Vectashield mounting medium.

29
Embryos were scored using DIC on a 10x/0.30 UPlanFL N objective (Olympus). DAPI stained 30 oocytes were imaged and scored between the spermatheca and posterior gonadal arm with a 31 40x/1.3 PlanApo objective (Olympus). Sperm and germline stem cells were identified using z-1 stack to capture multiple planes within each spermatheca or progenitor zone respectively. The 2 transition zone was identified using morphological criteria based upon the crescent shape of the 3 nucleus. All images and scoring were processed in ImageJ.

2
We solved these ODEs numerically using a Dormand-Prince explicit solver or using an 3 analytical solution (below). To calculate best fits, we assumed that the LSJ2, CX12311, RILs 4 and NIL nurf-1 strains had unique values of k o , but shared the k f parameters. These parameters 5 were estimated using a Levenberg-Marquardt non-linear least squares algorithm. To calculate 6 effect-size, we subtracted the two fits.

7
We used Mathematica to analytically solve the ODE equations: We used these equations to estimate the parameters using a Levenberg-Marquardt non-linear 10 least squares algorithm.