Anthropogenic selection along directions of most evolutionary resistance

Selection from harvesting, habitat fragmentation or climate warming is often bivariate against both a late maturity and large body size in parallel. This bivariate nature of anthropogenic selection remains largely overlooked, and the ability of populations to respond to it is poorly evaluated. Here, we use a laboratory experiment over 6 generations to measure the ability of wild-caught medaka fish (Oryzias latipes) to evolve in response to opposite anthropogenic selection regimes. Specifically, we selected medaka for maturity at 75 days-post-hatch and imposed random-(control), down- (anthropogenic, fishing-like) or up- (inverse anthropogenic) size selection. That is, we imposed a common positive selection differential on maturity but antagonistic selection differentials on body size, which generated antagonistic selection gradients on maturity as well, such that inverse anthropogenic-selected medaka experienced a negative selection gradient on maturity. As predicted from selection gradients, inverse anthropogenic-selected medaka evolved faster somatic growth and delayed maturation, but with no correlated change in survival rate, fertility, egg size or pituitary gene expressions of luteinizing, follicle-stimulating or growth hormone. In contrast, anthropogenic-selected medaka were unable to respond to selection gradients for a smaller body size and earlier maturity, indicating a functional limitation on this particular combination of traits. Our results demonstrate that humans may unintentionally select along directions of most evolutionary resistance in the phenotypic space, and that a consideration of the multivariate nature of selection is paramount to our ability to predict potential for evolutionary rescue in the face of global change.


INTRODUCTION
Human-induced perturbations such as habitat fragmentation, harvesting or climate warming often converge towards both increasing background mortality and selecting against a large body size in wild populations (Edeline, 2016). Such anthropogenic mortality regimes generate bivariate selective pressures acting in parallel on maturity and body size. Increased mortality, either size-independent or directed against large-sized individuals, directly selects for an earlier maturity by reducing the fitness of late-maturing individuals relative to early-maturing individuals . In parallel, Experimental harvesting experiments in the laboratory are potentially free of such problems because they make it possible to accurately target the traits under selection, to fully control the pattern and intensity of artificial selection, as well as to standardize environmental variation so that the effects of phenotypic plasticity are alleviated. Despite these obvious strengths of experimental approaches in the laboratory only a handful of such harvesting experiments have been performed in invertebrates  and fish van Wijk et al., 2013;. These experiments have shown that size-at-age or size at maturity in populations subject to small-vs. large-sized harvesting may van Wijk et al., 2013) or may not  evolve in the predicted direction. We provide a more detailed review of these experiments as Supplementary material.
Such harvest-simulating experiments remain few, and it is thus unclear whether exploited populations are always able to respond to size-selective harvesting. Importantly, none of these previous experimental designs have examined potential for bivariate selection acting both on maturity and body size in parallel. More generally, the many size-selective experiments performed so far on model organisms such as Drosophila (Partridge et al., 1999), chicken (Dunnington et al., 2013) or mice (Falconer & Mackay, 1996) have not selected on maturity in parallel. We are aware of one single study having examined response to parallel selection on body size and development time in the hawkmoth Manduca sexta (Davidowitz et al., 2012(Davidowitz et al., , 2016. Results show that traits responded roughly as predicted but the response of development time was weak and asymmetric (stronger response to selection for a long than for a short development time), making the bivariate response to selection essentially erratic.
Whether and how other populations or species can adapt to a bivariate anthropogenic selective regime remains unknown. Additionally, our understanding of the molecular and physiological regulations that underlay adaptive response to anthropogenic selection remains limited (Duffy et al., 2013;van Wijk et al., 2013;Uusi-Heikkila et al., 2017). This lack of knowledge not only hampers our ability to predict life-history evolution (Davidowitz et al., 2012(Davidowitz et al., , 2016, but also hinders the development of molecular diagnosis tools to evaluate potential for (and signature of) evolutionary rescue in wild populations. To contribute filling these gaps in our knowledge, we examined the ability of a wild population of medaka fish (Oryzias latipes) to respond to anthropogenic selection in the laboratory. To keep selected lines synchronised all fish were similarly selected for maturity at 75 days-post-hatch (dph), an age at which 86% of the fish were mature on average. We contrasted selection by applying down, up and random size-selection. Selection for maturity at 75 dph prevented late-maturing fish to reproduce, and thus mimicked a positive anthropogenic selection on maturity. Onto this high-mortality background common to all lines, the added contrast in selection differentials on body size generated antagonistic selection gradients both on body size and on maturity as well (see Results). While down size selection generated selection gradients for increased maturity probability and smaller body size, up size selection generated selection gradients for decreased maturity and larger body size. Selection gradients are indeed not necessarily of the same sign as selection differentials when traits are phenotypically correlated, as it is the case for size and maturity. Our size-selection protocol thus generated three selected lines through (i) bivariate AS for both earlier maturation and slower somatic growth in parallel, (ii) bivariate inverse-anthropogenic selection for both later maturation and faster somatic growth in parallel, and (iii) bivariate intermediate selection that served as a control (intermediate) selective (CS) regime by imposing mild selection for earlier maturation and as little as possible selection on somatic growth in parallel.
We selected medaka during 6 generations, measuring each generation a total of 14 phenotypic and neuroendocrine traits (Table 1). We isolated selected pairs and raised their offspring in individual tanks, keeping track of the pedigrees along the experiment. For the first time in a harvest-induced evolution experiment, this made it possible to standardize the number of offspring per individuals, maximize effective population sizes, and to control for the level of inbreeding throughout the selection procedure.
We made three specific predictions for medaka response to selection: (1) compared to CS medaka, AS medaka should evolve earlier maturation and slower somatic growth rates. We predicted and opposite pattern in IS medaka. (2) Life-history traits are often correlated, and selection on size-at-age only has been shown to induce correlated responses of reproductive and larval viability traits (Walsh et al., 2006). Therefore, we predicted that evolution in AS medaka should be paralleled by evolution towards a higher size-specific fecundity and/or larger egg sizes (Roff, 1992), as well as towards reductions in 5 115 120 125 130 135 140 size at hatch and larval survival (Walsh et al., 2006). We predicted an opposite response in IS medaka.
(3) The neuroendocrine control of vertebrate development involves pituitary hormones. Reproductive investment is stimulated by production of the luteinizing (LH), follicle-stimulating (FSH) and growth (GH) hormones, while somatic growth is stimulated by production of GH (Rousseau & Dufour, 2007;Zohar et al., 2010). Hence, compared to CS medaka we predicted altered GH, LH and FSH pituitary expressions in selected medaka, with potentially opposite alteration patterns in AS and IS lines.
Our results validate prediction (1) but in IS medaka only, because AS medaka did not show any response to selection. Prediction (2) was partially validated since egg size increased in AS medaka, but size-at-hatch decreased in both AS and IS medaka, and larval survival was not affected by selection.

Fish origin and maintenance
Our starting medaka population descended from 100 wild-caught individuals sampled in Kiyosu (Toyohashi, Aichi Prefecture, Japan) in June 2011. The genome of the Kiyosu population is free of any significant structure and shows a high degree of polymorphism, indicating no recent population bottleneck (Spivakov et al., 2014). These 100 breeders were maintained in five 20L aquariums and eggs were collected daily from July to September 2011. Hatched larvae were stocked in six 10 m 3 outdoor ponds.
In 2013, around 100 adult fish were transferred from outdoor ponds to the laboratory where all the 9 subsequent generations (dubbed F -1 to F 7 ) were maintained under a constant cycle of 14h of light -10h of darkness in 3L aquariums connected to a continuous flow-through system ensuring good water quality. Temperature was maintained between 26 and 27.5°C. Fish were fed ad libitum with a mixed diet of dry food (Marin start, Le Gouessant Aquaculture) delivered 4 times per day using automatic microfeeders (Eheim 3581), and live food (Artemia salina nauplii and/or Turbatrix aceti) delivered once a day, 5 days per week. These light, temperature and food conditions are expected to provide optimal growth and maturation conditions to medaka (Kinoshita et al., 2009).
Each generation, larvae were introduced in their aquariums at a density of 18.6 ± 3.9, 18.1 ± 3.5, 19.2 ± 3.0 (mean ± SD) larvae per aquarium in the control-selected (CS), anthropogenic-selected (AS) and inverse anthropogenic-selected (IS) lines, respectively. Selection was performed on aquariums of at least 10 individuals at 60 dph (see selection protocol below), in which initial larval densities per aquarium were 19.6 ± 1.6, 19.2 ± 1.9, 19.8 ± 1.0 (mean ± SD) in the CS, AS and IS lines, respectively.

Breeding design, pedigree and fish numbers
Fish initially transferred from outdoor ponds to the laboratory were allowed to mate randomly in groups of 3-6 fish per aquarium to produce the F -1 generation. In F -1 and F 0 , we randomly bred 54 (F -1 ) and 56 (F 0 ) pairs, respectively, to alleviate maternal and grand maternal effects, as well as to break any 7 160 165 170 175 180 genetic structure or linkage disequilibrium that could remain from the wild-caught population. In subsequent generations (F 1 to F 7 ) we proceeded with selection (see below). Eggs from each breeding pair were pooled for incubation in the same jar in a common recirculation system, and larvae from the same clutch were transferred to the same growth aquarium so as to form sibling families. This way, we were able to keep track of individual pedigrees so as to estimate individual inbreeding rate as 2k-1, where k is one's kinship coefficient with oneself (as calculated from the pedigree data using the kinship2 R package (Sinnwell et al., 2014)).
In each generation we performed selection (see below) on 636 fish on average (212 fish/line).
Specifically, we kept 4 breeders (2 males and 2 females) in each of 10 families per line to form the subsequent generation (20 breeding pairs/line/generation). By controlling pairing and normalizing the number of offspring per pair, we were able to maximize the effective population size. Additionally, we formed breeding pairs so as to minimize inbreeding. Assuming no inbreeding in F 1 , mean inbreeding rate was 9.6% (± 1.9 SD) by F 7 (11.5 ± 2.4% in the IS line, 8.5 ± 0.9% in the AS line, 9.1 ± 1.0% in the CS line). This corresponds to an average effective population size ("inbreeding effective numbers" sensu Crow & Kimura, 1970) of N e = 30.2 (31.7 in CS, 34.0 in AS and 24.8 in IS lines), i.e., corresponds to an unusually large N e for pedigree-recording artificial selection in a vertebrate species.

Phenotyping and hormonal measurements
Each generation, eggs from each breeding pair were collected during a period corresponding to mother's 88 to 92 dph. Eggs were counted and each clutch was photographed and ImageJ was then used to measure individual egg perimeters (9795 eggs measured from F 1 to F 7 ). Hatched larvae were collected during a 5-day time window so as to synchronize birth dates as much as possible. Individual age was calculated as the median hatching date of each sibling family.
At 0 (hatching), 15, 60 and 75 dph each single individual was photographed and ImageJ was then used to measure standard body length (from the tip of the snout to the base of the caudal fin, 16808 individual measurements from F 1 to F 7 ). Additionally, each individual at each phenotyping was sexed as immature (I), female (F) or male (M) according to their secondary sexual characters (Yamamoto, 1975). This method detected the onset of secondary sexual characteristics, which was a nondestructive proxy for the onset of maturity. All fish manipulations were performed after anaesthesia in tricaine In addition to phenotyping, a subsample of fish were individually measured for pituitary mRNA levels of β-subunits of gonadotropin hormones (LHβ and FSHβ) and GH. At about 40 dph in each generation from F 1 to F 7 , 10 to 15 fish per line were randomly sampled and dissected for endocrine measurements (233 fish measured for all three hormones from F 1 to F 7 ). F 0 preliminary data indicated that the onset of secondary sexual characteristics occurred roughly between 40 and 60 dph, and we chose to dissect fish at 40 dph so as to sample fish at the initiation of puberty. Fish were phenotyped as described above, sacrificed and dissected under a binocular microscope for the pituitary which was immediately immersed in 250 µL Trizol (Ambion) and stored at -20°C. Pituitary mRNA levels of LHβ, FSHβ and GH were measured using reverse transcription quantitative real-time polymerase chain reaction (RT-qPCR). Further details on the RT-qPCR procedure are provided in Supplementary Methods.

Selection differentials and gradients
We applied a common, positive selection differential on maturity and a contrasted, line-dependent selection differential on body length that, combined, produced a control-selected (CS), an anthropogenic-selected (AS), and an inverse-anthropogenic (IS) lines. Specifically, the size-dependent selection differential was applied both on families at 60 dph and on mature individuals at 75 dph, an age at which on average 86% of the fish were mature. At 60 dph, we discarded families of less than 10 individuals (viability selection) and, among the remaining families, we kept 10 families at random (CS) or that had the smallest (AS) or largest (IS) average standard body length. At 75 dph, we individuallyselected breeders (see above) among mature fish only based on their individual standard body length, generating selection differentials of +12.5% maturity probability and +0.71 mm, +10.4% maturity probability and -0.89 mm, and +12.8% maturity probability and +2.12 mm in the CS, AS and IS lines, respectively. This selection procedure resulted in keeping on average 12% of individuals per line at each generation (number of breeders / total number of fish before selection at 75 dph).
Selection differentials are not a proper measure of selection when selection acts on several correlated characters (here body size and maturity showed a 53% correlation at 75 dph), in which case the measurement of selection should rely on the estimation of selection gradients (Lande & Arnold, 1983;9 220 225 230 235 240 245 Phillips & Arnold, 1989). We estimated directional selection gradients as the coefficients of the linear regression of standard body length (mm) and maturity coded as 0 or 1 (both strandardized to zero mean) on relative fitness ω (%) defined as ω i =ν i /ν G [i] , L [i ] where ν i is the number of progeny produced by individual i surviving to age 75 dph in the next generation, and ν G [i ], L[ i] is mean absolute fitness in the generation and line (population) to which individual i belongs.

Line replication
In the context of selection experiments where the number of individuals is limited, line replication trades off with increasing effective population size N e . Maximizing N e should prime, because a large N e decreases genetic drift, limits the effect of linkage disequilibrium on selection limits, and delays the unavoidable increase in inbreeding (Robertson, 1960;Hill & Robertson, 1966), see e.g. Weber & Diggins (1990) for experimental evidence. In particular, avoiding genetic drift and inbreeding is crucial when studying the evolution of correlated characters (Phillips et al., 2001). Given this trade off, and because the logistical constraints linked to controlled pairing, pedigree tracking and intensive phenotyping limited the total number of fish to be processed at each generation to an average 636/generation, we chose to derive three large-population lines (3 x N e = 30) rather than replicating small-population treatments. By doing this, we lost the ability to measure experimentally the respective contributions of drift vs. selection to the observed phenotypic (non)divergence among the lines -note, however, that we could still estimate it theoretically from the effective population size, which is precisely known from the pedigree. In any case, an accurate measurement of the phenotypic divergence due to genetic drift requires a very high number of line replicates (Lynch, 1988).
Previous laboratory harvesting experiments in fish have used line duplicates in simple designs with no control for effective population size, pedigrees or inbreeding, but however report adult population size N. Based on a median N e /N = 0.23 in random-mating populations (Palstra & Fraser, 2012), N e was expected to range from 5.5 to 17.9 in   Table 1 in Results below).

Data analyses
The aim of our statistical analyses was to estimate and test significance for the overall effect of anthropogenic selection on medaka life-history and neuroendocrine traits. A visual appreciation of time series response to selection on body length (Fig. 2B) shows that the divergence between the three selected lines somehow stabilized from generation F 3 . Hence, in our analyses we pooled data from generations F 3 to F 7 and treated generation as a random effect. A detailed description of the statistical models used is provided as Supplementary Methods.
We visualized the effect of anthropogenic selection on the maturation process using probabilistic maturation reaction norms (PMRNs). Maturation reaction norms have been developed to account for the plastic effects of the growth trajectory on the maturation process, such that a shift in the maturation reaction norm may be interpreted as an evolutionary shift in maturation (Stearns & Koella, 1986).
Technically, the maturation reaction norm approach is most robustly implemented by modelling the maturity data using logistic regressions, in which case maturation reaction norms become so-called PMRNs (Heino et al., 2002). The PMRN approach has often been applied to assessing potential for evolutionary shifts in the maturation process of exploited fish populations (Heino & Dieckmann, 2008).
We computed line-specific PMRNs, defined as the length-and age-dependent 50% probability for an immature medaka to initiate maturity (as informed by the onset of secondary sexual characteristics), using the methods of  and Van Dooren et al. (2005). The methods consisted in (1) computing maturity "ogives", (2) computing maturation probabilities and (3) computing line-specific PMRNs. More details are provided as Supplementary Methods.

MCMC parameter estimation
All statistical models used in analysing data (Table S2) were fitted using Markov chain Monte Carlo (MCMC) in JAGS 4.2.0 (Plummer, 2003). For each model, we ran three independent MCMC chains thinned at a period of 5 iterations until parameter convergence was reached, as assessed using the Gelman-Rubin statistic (Gelman & Rubin, 1992). Weakly informative priors for regression parameters were defined as normal distributions with zero mean and 1000 standard deviation and for variance 11 285 290 295 300 305 310 parameters as a uniform distribution between 0 and 100.
We tested the significance of effects from posterior parameter distributions in a test equivalent to a twoway t test. In these tests, the MCMC P-value was twice the proportion of the posterior for which the sign was opposite to that of the mean posterior value. We further assessed goodness of fit of our models by using a Bayesian P-value (Gelman et al., 1996). Briefly, we computed residuals for the actual data as well as for synthetic data simulated from estimated model parameters (i.e., residuals from fitting the model to ''ideal'' data). The Bayesian P-value is the proportion of simulations in which ideal residuals are larger than true residuals. If the model fits the data well, the Bayesian P-value is close to 0.5.
Bayesian P values for our models ranged from 0.49 to 0.66, indicating excellent model fit. All models were fitted using an "effect" parametrization (Table S2), i.e., by setting one level of each factor as a reference levels as is done by default in R. All data and codes used in this study are provided as Supplementary Material. 12 315 320

Selection differentials and gradients
As expected from our selection protocol, selection differentials were different among the selected lines for standard body length but not for maturity, i.e., all lines similarly experienced a selection differential for increased maturity probability but contrasted selection differentials on body size (Fig. 2, grey, Tshaped arrows). In contrast to selection differentials, selection gradients show that our selection procedure produced also contrasted selection gradients on maturity among the selected lines ( Fig. 2, black, curved arrows). Selection gradients, not differentials are needed to measure selection on correlated traits such as maturity and body size (Lande & Arnold, 1983;Phillips & Arnold, 1989).
Therefore, while control-selected (CS) fish experienced a mild positive gradient for maturity, anthropogenic-selected (AS) medaka were strongly selected for maturity because small-sized mature fish were few. In contrast, anthropogenic-selected (IS) fish, despite experiencing a positive selection differential on maturity, were in fact negatively selected on maturity because largest-sized medaka were on average larger than the bulk of mature medaka.
Specifically, from F 1 to F 7 , selection gradients on maturity expressing % change in relative fitness with a change from immaturity to maturity were on average +81%, +215% and -65% in the CS, AS and IS lines, respectively. Selection gradients on standard body length expressing % change in relative fitness with a 1 mm increase in standard body length were on average +4%.mm -1 , -27%.mm -1 , +43%.mm -1 in the CS, AS and IS lines, respectively. Relative to CS medaka, selection gradients were +134% and -146% on maturity and -31%.mm -1 and +40%.mm -1 on standard body length in the AS and IS lines, respectively. Therefore, relative to the CS line our artificial selection design accurately reproduced the typical bivariate anthropogenic selection for both slower somatic growth and earlier maturation in the AS line and an almost exact opposite pattern in the IS line.
We used simulations to more deeply explore the conditions that may give rise to a negative selection gradient on maturity despite a positive selection differential (R code provided as Supplementary material). We simulated individual growth and maturation trajectories in which individuals varied randomly for somatic growth rate, asymptotic length, and the length-and age-dependent probability to mature. We then applied bivariate selection at 75 dph by removing all immature individuals (positive differential on maturity) and by keeping only the largest individuals among mature (positive differential 13 330 335 340 345 350 on body size). Computation of selection gradients show that selection gradients on maturity become negative when the body size-maturation correlation and the strength of selection on body size both increase (Fig. S2).

Body size and maturation response to selection
The average body length in all lines was substantially affected by a generation environmental effect, with no obvious trend (Fig S1). In contrast, these was a substantial trend towards decreased maturity probability at 75 day-post-hatch (dph) in all the lines (not shown). In order to account for these environmental effects, all subsequent analyses were performed by comparing selected lines with the CS line. Qualitative (directions of) responses to selection for the 14 measured traits are summarized in Table 1, while quantitative statistical results are provided in Supplementary Table S2.
As predicted from selection gradients, IS medaka evolved towards a larger standard body length at 75 dph in both mature (Figs. 3A and S1A, Table S2) and immature fish (Figs. 3B and S1B, Table S2).
Additionally, compared to CS, IS medaka had a lower maturity probability at an average body length and age (Table S2), indicating evolution towards delayed maturity. We visualized this effect using probabilistic maturation reaction norms (PMRNs), which show the estimated combination of age and lengths at which maturation probability is 50%. The PMRN itself and its slope were larger in IS compared to CS medaka (Fig. 4), confirming that AS medaka had evolved delayed maturation compared with CS medaka.
In contrast to predictions, neither standard body length at 75 dph (Figs. 3 and S1, Table S2) nor maturity probability at an average body length and age (Table S2) responded to selection in AS medaka. Accordingly, the PMRNs of CS and AS medaka largely overlapped (Fig. 4), indicating that overall maturation response to selection in AS medaka was negligible. Therefore, our prediction (1) was validated only in IS and not in AS medaka, indicating an asymmetric life-history response to anthropogenic selection.

Response to selection in non directly-selected traits
Our prediction (2) was that evolution in selected life-history traits should be paralleled by evolution in correlated traits, and in particular in size-specific fecundity, egg sizes, size at hatch and larval survival 14 360 365 370 375 380 (Roff, 1992;Walsh et al., 2006). However, only one out of 12 non directly-selected traits did respond to selection. Egg size was significantly increased in AS compared to CS medaka (Tables 1 and S2).
In contrast with prediction (2), we found that medaka body length at hatch was significantly decreased in both AS and IS medaka compared to CS medaka, and we found no effect of selection on survival (Table 1, Table S2). Noticeably, body length at hatch was also the only of the 14 monitored traits that was significantly influenced by inbreeding, more inbred individuals having a larger size at hatch (Table   S2). None of the other traits investigated showed any significant effect of inbreeding (Table S2).
Our prediction (3) was also invalidated since pituitary expression levels of growth hormone (GH), of the β subunits of luteinizing hormone (LH) and follicle-stimulating hormone (FSH) were among the traits that were not significantly influenced by selection (Table 1, Table S2). Additionally, pituitary gene expressions for the three hormones were highly positively correlated (Table S2).

DISCUSSION
Medaka response to bivariate selection on maturity and body size was asymmetric. simply be impossible in some populations or species. Such an asymmetric response to selection was not found in previous experiments on Atlantic silverside , zebra fish by  or guppy (van Wijk et al., 2013), but compares with the results of  in zebra fish showing that the magnitude of response to size-dependent selection was traitspecific and contingent upon the direction of selection. More generally, artificial selection experiments in insects suggest that weak, asymmetric or erratic responses to selection are commonplace when development time is selected simultaneously with a morphological trait (Scheiner & Istock, 1991;Zijlstra et al., 2003;Fischer et al., 2006;Davidowitz et al., 2016). Asymmetric and reversed responses might even be a general (and not well understood) feature of responses to multivariate selection (Lynch & Walsh, 2018 p.607).
In our experiment, lack of response to selection in AS medaka could not be ascribed to an absence of selection (cumulated selection gradients in the AS relative to CS line were +802% on maturity and -183% on body length at 75 days-post-hatch), nor due to differences in survival, fecundity or population size which were all identical among the selected lines, nor due to inbreeding which was by F 7 identical among the CS and AS lines. Instead, the absence of evolution in AS relative to CS medaka lines suggests that physiological, developmental or functional constraints limit the range of pathways that life-history evolution may follow in medaka (Lande, 1979;Arnold, 1992;Schluter, 1996).
Here in particular, it is apparent that there is little (no) scope for the evolution of a combined smaller body size and earlier maturity in medaka. This particular functional constraint (sensu Arnold, 1992) is likely to be commonplace in short-lived, r-selected species because maturation is expected to occur shortly after reaching a functional size limit in these species. In contrast, in long-lived, K-selected species, delayed maturation allows for an evolutionary response to r-like selection, and thus increases probability for evolutionary rescue in response to anthropogenic selection. In Atlantic cod (Gadus morhua) for instance, several decades of rapid evolution towards earlier maturation and slower somatic growth preceded stock collapse (Olsen et al., 2004;Swain, 2011), suggesting that evolution might have rescued the populations until a functional size limit was reached.
So far, most theoretical, empirical or experimental studies of harvest-induced evolution have focused on univariate selection acting either on maturity or on size-at-age. This assumption of univariate 16 420 425 430 435 440 445 selection simplifies analyses but seems unrealistic since human-induced selection is often (at least) bivariate by acting both on maturity and body size in parallel. The importance of considering multivariate selection and trait correlations to properly predict evolutionary trajectories has long been stressed (Lande, 1979;Lande & Arnold, 1983;Phillips & Arnold, 1989;Arnold, 1992;Schluter, 1996).
Our results provide an illustration of this fact in the context of understanding and managing the response of populations to anthropogenic perturbations. An ubiquitous inability to respond to bivariate selection for both an earlier maturity and smaller body size in other species than medaka might contribute to explain apparent lack of phenotypic response to harvesting in some wild fish stocks (Devine & Heino, 2011;Silva et al., 2013;Marty et al., 2014), why the occurrence of harvest-induced evolution remains controversial (Borrell, 2013), or why examples of evolutionary rescue are so few in the wild (Carlson et al., 2014).

Hormonal response to selection
As a first approach to uncovering the molecular regulation of adaptive life-history evolution in medaka, we used RT-qPCR measurement of candidate genes in the pituitary. We specifically targeted genes known to play a central role in the regulation of somatic growth and reproduction. In teleosts, growth hormone (GH) is a pleiotropic pituitary hormone related not only to somatic growth rate (Reinecke et al., 2005;Canosa et al., 2007) but also to maturation, osmoregulation and stress (Le Gac et al., 1993;Rousseau & Dufour, 2007). Its action can be direct by stimulating cell growth in skeletal muscles expressing its receptor, or indirect by stimulating insulin-like growth factors expression and secretion in liver (Rousseau & Dufour, 2007). The luteinizing (LH) and follicle-stimulating (FSH) hormones are known to stimulate steroidogenesis and gametogenesis in teleosts as in other vertebrates (Zohar et al.,

2010).
We expected pituitary expression of LH, FSH and GH to be altered in response to selection in AS and IS medaka. Instead, despite we did observe the expected phenotypic response to selection in IS medaka we could not detect any significant neuroendocrine response, suggesting that pituitary gene expression may not be a major target of size-and maturity-dependent selection in the medaka. Instead, IS selection may have primarily acted downstream of the pituitary such as at the pituitary hormone receptor level, on IGF and its receptors (Duffy et al., 2013) (Rousseau & Dufour, 2007). Interestingly, pituitary activity of the somatotropic and gonadotropic axes were highly positively linked in medaka, suggesting that somatotropic and reproductive hormones were synergistic in their effects on medaka development. Future transcriptomic approaches on central and peripheral tissues will provide a deeper understanding of the molecular regulation of response to anthropogenic selection in our medaka lines.

Conclusions
Asymmetric response to anthropogenic selection in medaka is a warning signal that calls for increasing research efforts to assess life-history evolvability in the wild. A crucial line of work in achieving this goal will consist in accurately measuring the multivariate components of selection that act on correlated life-history traits, both in the wild and in laboratory experiments. The other key element of this effort will rely on developing diagnosis tools to evaluate potential for (and signature of) adaptive response to anthropogenic selection. In the future, comprehensive approaches melting wide-spectrum candidate genes, transcriptomics and genome scans of experimentally-selected and wild populations will probably be needed to finely decipher the molecular architectures that regulate the adaptive evolution of life histories and ultimately support the maintenance of biodiversity and ecosystem productivity. 18 Acknowledgments. We are grateful to Prof. Kiyoshi Naruse (NIBB, Okazaki, Japan) for his support in obtaining and maintaining medaka from Toyohashi. We are also thankful to Prof. Finn-Arne Weltzien for providing several primer sequences for our candidate genes. We thank the people who helped us at the laboratory: Beatriz Decencière, Julien Hirschinger, Alice Lamoureux, Alexandre Macé and Yohann Chauvier.
Funding. This work has benefited from technical and human resources provided by CEREEP-Ecotron IleDeFrance (CNRS/ENS UMS 3194) as well as financial support from the Regional Council of Ile-de-

Author contributions
EE, ALR, GM and SD designed the study. CR, DC, AM, SA maintained the fish and performed selection, breeding, phenotyping and data collection. CR and GM performed mRNA measurements. EE and ALR made final analyses and finalized paper writing.

Ethical statement
The protocols used in this study were designed to minimize discomfort, distress and pain of animals, and were approved by the Darwin Ethical committee (case file #Ce5/2010/041).

Data archiving statement
All data and codes used in this paper will be archived online.

Competing interests
The authors declare no competing financial interests.     Represented phenotypes, differentials and gradients were averaged across generations from F 1 to F 7 .
Selection gradients were scaled to phenotypic units by multiplying gradient estimates by phenotypic variances.   have applied small-vs. large-sized harvesting to a standard mixture of Daphnia magna clones during a 150 day period. Clones exposed to small-harvesting evolved rapid somatic growth through small size classes and delayed maturation, while clones exposed to large-harvesting evolved a slow growth through small size classes and earlier maturation. Computation of reproductive values showed that evolution resulted in a redistribution of reproductive investment towards size classes that were not harvested.  exposed soil mites Sancassania berlesei to juvenile or adult harvesting during 70 weeks. In line with theoretical predictions , they found that juvenile harvesting induced evolution towards earlier maturation, while adult harvesting induced evolution towards delayed maturation. Interestingly, they also found that the amplitude of harvest-induced evolution was overwhelmed by evolution to delayed maturation in all treatments. This change was interpreted by authors as a response to the captive environment, in which density and competition for resources were increased compared to the natural environment from where mites were initially sampled.  applied small-vs. large-sized harvesting at 190 days postfertilization during 5 generations in six experimental populations of the Atlantic silverside Menidia menidia maintained in 700 L tanks (about 100 breeders/generation/population, random mating, unknown effective population sizes). The Atlantic silverside is an annual fish, and it was assumed that all individuals were mature at selection such that selection was imposed on somatic growth only.  found that the mean weight of fish evolved in the expected direction such that, by generation F 5 , an average fish aged 190 days postfertilization weighted 4.5 g in the small-harvested line but 2.5 g in the largeharvested line (3.5 g in the random-harvested line). These differences were due to differences in somatic growth rate and underlying traits (Walsh et al., 2006). Zebra fish were harvested at an age corresponding to 50% of mature fish in the random-harvested line and reproduced 14 days later. Response to selection was contingent upon both the trait considered and upon the direction of selection. Compared to the random-harvested line, the small-harvested line showed no change in juvenile somatic growth rate or asymptotic length and matured at a later age (but not size), while the large-harvested line showed no change in juvenile somatic growth rate but evolved lower asymptotic length and maturation at a smaller size (but not age).

Reverse transcription quantitative real-time polymerase chain reaction (RT-qPCR)
After sample homogenization by agitation (15 sec vortexing), total RNA were extracted according to the manufacturer's indications, suspended in 10 µL RNAse-free water, and treated with DNAse I (Dnase I recombinant RNAse-free, Roche Diagnostics). Then cDNA were produced from 5 µL of total RNA using RT Superscript III (RT Superscript III First Strand cDNA Synthesis Kit; Invitrogen, Life Technologies) and random hexamer primers (50 ng; Invitrogen, Life Technologies), at 50°C for 60 min after an initial step of 25°C for 10 min. Medaka specific primer sets for FSH were designed with primer3 software (Koressaar & Remm, 2007) on two successive exons, or on exon junctions. Gene specific primer sets for LHβ, FSHβ, GH and actin-β (used as reference gene to correct for technical noise), were previously designed (see Supporting Information, Table S1). Efficiency and amplification specificity were checked for each primer set. The sets with the highest efficiency were chosen for the following quantification experiment.
Messengers RNAs were assayed using Light Cycler system (LightCycler® device; Roche Diagnostics) with the LightCycler FastStart Master plus SYBR Green I kit (Roche Diagnostics) as recommended by the manufacturer, from 4 µL of diluted 1:10 cDNA samples and the specific primers concentrated at 500 nM (Eurofins). The PCR conditions were 95°C for 10 min followed by 50 cycles at 95°C for 5 sec, 60°C for 10 sec and 72°C for 5 sec.
Expression levels of mRNA for LHβ, FSHβ, GH and actin-β in each individual fish were measured in duplicate using the "relative quantification" method (Applied Biosystems User Bulletin #2 1 ). Briefly, the standard relationships between fluorescence and gene-specific sample RNA concentrations were constructed using a bulk RNA pool, hereafter dubbed "calibrator". The Lightcycler software estimated the number C q of quantification cycles needed to reach the inflection point (second derivative equal to 0) of fluorescence amplification for a series of 7 calibrator volumes. From this, the software estimated the intercept and slope parameters for the gene-specific, linear relationship between log10 calibrator volume and C q . These linear relationships were then used to predict sample-specific mRNA expressions in log10 calibrator volume units (i.e., "arbitrary" units) for LHβ, FSHβ, GH and actin-β from their 1 http://www3.appliedbiosystems.com/cms/groups/mcb_support/documents/generaldocuments/cms_040980.pdf 32 660 665 670 675 680 sample-specific C q . At each PCR run, a known amount of the calibrator plus a blank (water) were measured for C q so as to adjust for possible inter-run noise. Following standard practises, we used as input data relative pituitary gene expression, calculated as the natural logarithm of the ratio between mRNA expression for the interest gene and mRNA expression for actin-β (see model 5 below). In particular, this approach corrects for the effects of variability in pituitary size.

Statistical analyses
The aim of our statistical analyses was to estimate and test significance for the overall effect of selection on the selected phenotypes (somatic growth and maturity probability) and correlated lifehistory and physiological traits. We pooled data from generations F 3 to F 7 and treated generation as a random effect, because a visual appreciation of response to selection  showed that the divergence between the three selected lines somehow stabilized from generation F 3 .

Standard body length at 75 dph
We modelled response to selection as the effect of selected line on Sdl at 75 dph ( Sdl 75 ) of each individual i : where Sdl 75 is standard body length of mature fish at 75 dph, N is the normal distribution, We computed maturity "ogives" as: where y i is the maturity status of an individual fish i (0 or 1), Bern is the Bernoulli distribution of "success" (maturity) probability p , ln is the natural logarithm, Sdl is standard body length. Other subscripts or variables are as described above. By letting the effects of both Sdl and Age on p varying for each selected line this model captured potential effects of selection on both the intercept and slope of the PMRN.
We than computed the maturation probability m(a τ , s τ ) at each growth increment τ as described by : where o(a τ , s τ ) is age-and length-dependent maturity probability at the end of growth increment τ as predicted by model 2. We did so for slow, median and fast simulated growth curves (Van Dooren et al., 2005;Harney et al., 2013).
Finally, we computed line-specific PMRNs as the age and length combination (a t , s t ) at witch maturation probability reached 50%, i.e., as the age and length combination that satisfied the following condition (Van Dooren et al., 2005;Harney et al., 2013): (1−m(a τ , s τ ))=0.5 .
We estimated (a t , s t ) using 200 growth increments equally spread between ages 0 and 87 dph. Full error distribution for (a t , s t ) was obtained by iterating the procedure for each MCMC sample of the 34 715 720 725 parameter set in model 2.

Survival
We tested for differential mortality among the selected medaka lines using models of the form: where N (t +1) i is the number of individuals still alive at time t+1 in sibling family i , Bin is the Binomial distribution, Mo . inb is mother inbreeding coefficient, Fa. inb is father inbreeding coefficient. Other subscripts are as described above. We fitted separately four models for t to t+1 steps corresponding to the egg-larvae (egg-to-0 dph), larvae-juvenile (0-to-15 dph), juvenile-adult (15-to-60 dph) and adult-adult (60-to-75 dph) transitions. σ G 2 is as described above and ϵ is an overdispersion effect accounting for the fact that observed variance was larger than canonical variance of the Binomial distribution.

Size-specific fertility and fecundity
Our aim here was to test for possible effects of size-dependent selection on medaka size-specific fecundity. Counts F of clutch size per breeding pair i were zero-inflated Poisson-distributed and modeled as: where Pois is the Poisson distribution with mean (and variance) equal to the product of probability for a nonzero count 1−π with nonzero counts λ . This way, we were able to simultaneously test the effects of the predictors both on the probability for a breeding pair to be infertile ( p ) and on the fecundity of a fertile pair ( λ ).
offset is the natural logarithm of number of days during which eggs were collected (varied from 4 to 5 days), Mo. sdl is mother's standard body length (here ln-transformed to linearize the fecunditysize relationship), and ϵ is an overdispersion effect accounting for the fact that variance overwhelmed the mean in non-zero medaka egg counts. Other subscripts or variables are as described above. In this model, the effect of mother Sdl is accounted for, such that a significant effect of the selected line would indicate that size-dependent selection affects medaka fertility or fecundity beyond direct effects on Sdl.

Egg size
Egg size was automatically measured in ImageJ as individual egg perimeter Pm i in mm and modelled as follows: where the variables are as described above.

Incubation time
Incubation time for eggs from each breeding pair It i was computed as the time laps (days) between weighted mean data of spawning and weighted mean date of hatching for larvae collected from mother's 95 to 100 dph. We evlauted the effect of selection on It i in: where variables are described above.

Larval size at hatch
We model individual standard body length at hatch Sdl 0 as: Model 7, where variables are described above.

Hormonal profile
Measurements