Unidirectional response to bidirectional selection on body size. I. Phenotypic, life‐history, and endocrine responses

Abstract Anthropogenic perturbations such as harvesting often select against a large body size and are predicted to induce rapid evolution toward smaller body sizes and earlier maturation. However, body‐size evolvability and, hence, adaptability to anthropogenic perturbations remain seldom evaluated in wild populations. 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 bidirectional size‐dependent selection mimicking opposite harvest regimes. Specifically, we imposed selection against a small body size (Large line), against a large body size (Small line) or random selection (Control line), and measured correlated responses across multiple phenotypic, life‐history, and endocrine traits. As expected, the Large line evolved faster somatic growth and delayed maturation, but also evolved smaller body sizes at hatch, with no change in average levels of pituitary gene expressions of luteinizing, follicle‐stimulating, or growth hormones (GH). In contrast, the Small medaka line was unable to evolve smaller body sizes or earlier maturation, but evolved smaller body sizes at hatch and showed marginally significant signs of increased reproductive investment, including larger egg sizes and elevated pituitary GH production. Natural selection on medaka body size was too weak to significantly hinder the effect of artificial selection, indicating that the asymmetric body‐size response to size‐dependent selection reflected an asymmetry in body‐size evolvability. Our results show that trait evolvability may be contingent upon the direction of selection and that a detailed knowledge of trait evolutionary potential is needed to forecast population response to anthropogenic change.


| INTRODUC TI ON
Human activities often converge toward selecting against large-bodied individuals in animal populations, mainly through harvesting, habitat fragmentation, and climate warming (Edeline, 2016). In this context, the dynamics of wild populations may critically rely on their capacity to evolve in response to size-dependent selection.
Whether and how wild populations can respond to anthropogenic size-dependent selection have been mostly explored in the context of fisheries, which are often highly size-selective Kuparinen, Kuikka, & Merilä, 2009;Lagler, 1968;Law, 2000). Harvesting large-bodied individuals is predicted to induce adaptive evolution toward earlier maturation through reduced life expectancy and, at the same time, toward slower somatic growth through selection against a large body size at a given age (Heino, Díaz Pauli, & Dieckmann, 2015). Paradoxically, however, selection for an earlier maturation may also result in evolution of faster somatic growth, which allows for an earlier maturation (Diaz Pauli, Kolding, Jeyakanth, & Heino, 2017;Dunlop, Heino, & Dieckmann, 2009;Eikeset et al., 2016). This result highlights the importance of considering trait correlations and multivariate phenotypes in evolutionary biology.
In the wild, fishing has been associated with phenotypic changes toward earlier maturation at a smaller body size and/or toward slower growth rates (see reviews by Fenberg & Roy, 2008;Heino et al., 2015;Kuparinen & Merilä, 2007;Law, 2000;Trippel, 1995).
Yet, cases of stocks with no phenotypic response to fishing are also reported (Devine & Heino, 2011;Marty, Rochet, & Ernande, 2014;Silva, Faria, & Nunes, 2013), suggesting that harvested populations might not always be able to respond to harvest-induced selection.
Studies based on data from the wild, however, are often criticized for problems in measuring actual selection pressures (but see Carlson et al., 2007;Edeline et al., 2007;Kendall, Hard, & Quinn, 2009), in disentangling the effects on mean trait values of size-selective mortality versus evolutionary changes (Hairston, Ellner, Geber, Yoshida, & Fox, 2005), or in controlling for the confounding effects of phenotypic plasticity (Heino, Dieckmann, & Godø, 2002). Hence, there is still debate as to whether changes (or absence thereof) toward earlier maturation and slower somatic growth in exploited populations are genetic (Borrell, 2013) or are occurring rapidly enough to influence population dynamics and thus probability of population persistence (Diaz Pauli & Heino, 2014). 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, and to standardize environmental variation so that the effects of phenotypic plasticity are alleviated.
Size-selective experiments have been performed on model organisms such as Drosophila melanogaster (e.g., Partridge, Langelan, Fowler, Zwaan, & French, 1999), chicken Gallus gallus (Dunnington, Honaker, McGilliard, & Siegel, 2013), or mice Mus musculus (e.g., Macarthur, 1949). Often, selection is bidirectional, that is, is performed at random (Control line), against a small body size (Large line) and against a large body size (Small line, mimicking the effects of harvesting). Results from these experiments show that bodysize response to selection may sometimes be asymmetric, with either the Large or Small lines showing slower, or sometimes no or halted response to selection (Falconer & Mackay, 1996 and references therein; Dunnington et al., 2013;Lynch & Walsh, 2018 and references therein). Additionally, selection on body size may be associated with changes in other traits. For instance, selection for increased thorax length in D. melanogaster was associated with an increase in larval development time and no change in somatic growth rate, while selection for reduced thorax length was associated with reduced growth rate but no change in duration of larval development (Partridge et al., 1999). Similarly, experiments specifically designed to simulate harvesting on wild populations of model or nonmodel organisms have shown that size at age or size at maturity in populations subject to small-versus large-sized harvesting may (Amaral & Johnston, 2012;Cameron, O'Sullivan, Reynolds, Piertney, & Benton, 2013;Conover & Munch, 2002;Edley & Law, 1988;van Wijk et al., 2013) or may not (Uusi-Heikkilä et al., 2015) evolve in the direction imposed by selection (see the Discussion for a more detailed treatment of these harvest-simulating experiments). Hence, so far our knowledge of whether and how exploited populations can respond to size-selective harvesting remain limited.
To contribute filling this gap in our knowledge, we examined the ability of a wild population of medaka fish (Oryzias latipes) to respond to bidirectional size-dependent harvesting in the laboratory.
Specifically, we selected medaka randomly (Control line), against a large body size (Small line), and against a small body size (Large line) during 2.5 years (30 months, 6 medaka generations), measuring at each generation a total of 14 phenotypic, life-history, and neuroendocrine traits.
We made three specific predictions for medaka response to size-dependent selection: (a) Compared to the Control line, medaka from the Small line should evolve slower somatic growth rates. We predicted an opposite pattern in the Large medaka line. (b) Selection on body size has often been shown to induce correlated responses of reproductive traits and larval viability (e.g., Walsh, Munch, Chiba, & Conover, 2006). Therefore, we predicted that evolution of somatic growth in the Small medaka line should be paralleled by evolution toward increased reproductive investment, which may result in earlier maturation and/or higher fecundity at a given body size and/or larger egg sizes (Roff, 1992), and/or toward reduced size at hatch and larval K E Y W O R D S anthropogenic selection, body size, evolvability, fisheries, life history survival (Walsh et al., 2006). We predicted an opposite response in the Large medaka line. (c) The neuroendocrine control of vertebrate body growth and reproduction involves production of the growth hormone (GH), luteinizing hormone (LH), and follicle-stimulating hormone (FSH) in the pituitary (Rousseau & Dufour, 2007;Zohar, Munoz-Cueto, Elizur, & Kah, 2010). Hence, compared to Control line we predicted altered GH, LH, and FSH expression levels in the pituitary, with potentially opposite alteration patterns in the Small and Large medaka lines. Our results validate prediction (a), but in the Large medaka line only, because the Small line did not show any body-size response to selection. Prediction (b) was validated in the Large line, but only partially in the Small line that did not mature earlier but showed signs of increased reproductive investment. Finally, prediction (c) was mainly not supported since only the pituitary expression GH showed a marginally significant response to size-dependent selection.

| Fish origin and maintenance
Our start medaka population descended from 100 wild-caught individuals sampled from a single population 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 20 L 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 60 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 constant environmental conditions (common garden): 3 L aquariums connected to a continuous flow-through system ensuring good water quality, cycle of 14h of light-10h of darkness, and temperature 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 3,581), and live food (Artemia salina nauplii and/or Turbatrix aceti) manually delivered once a day, 5 days per week. These light, temperature, and food conditions provide optimal growth and maturation conditions to medaka (Kinoshita, Murata, Naruse, & Tanaka, 2009).

| Breeding design, pedigree, and fish numbers
Prior to starting selection, we bred medaka during two generations in the laboratory to alleviate maternal and grand maternal effects (a diagram of the experimental design is provided in Appendix 2).
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 mated 54 (F −1 ) and 56 (F 0 ) pairs, respectively (Appendix 2), to break any genetic structure or linkage disequilibrium that could remain from possible assortative mating in the wild population (Lynch & Walsh, 2018).
Each generation, eggs from each breeding pair were pooled for incubation 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 and 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, Therneau, & Schaid, 2014).
Offspring from multiple breeding pairs were never mixed in the same aquarium (not to break the pedigree), and the aquarium and the sibling family effects were confounded. Occasionally, a breeding pair produced many progeny that were spread across two different aquariums (118 breeding pairs, out of 375 pairs in total, produced two aquariums of progeny). Aquariums were randomly spread across two different racks such that the selected lines shared the same environmental conditions. Larvae were initially introduced in their aquariums at a controlled density of 19.6 ± 1.6, 19.2 ± 1.9, and 19.8 ± 1.0 (mean ± SD) larvae per aquarium in the Control, Small,

| Selection procedure
We proceeded with selection on the F 1 -F 7 generations from February 2014 to August 2016 (30 months, Appendix 2). A size-dependent selection differential was applied both on families at 60 dph and on mature individuals at 75 dph, an age at which 86% of the fish were mature on average (for dynamics of maturity in each line, see Le Rouzic et al., in press).
At 60 dph, we discarded families of less than 10 individuals to avoid confounding density effects on phenotypes. This procedure generated significant selection for a higher fecundity (overdispersed Bernoulli GLM, discarded ~ fecundity, p-value < .01) and for higher survival rate from egg to age 15 dph (p-value < .005), but not for a larger or smaller body length (p-value = .296). Among the remaining families, we kept 10 families at random (Control line) or that had the smallest (Small line) or largest (Large line) average standard body length.
At 75 dph, we individually selected breeders among mature fish based on their individual standard body length and precluded brothersister mating. Specifically, we kept in each family 4 mature fish (2 males and 2 females) that were paired with breeders from other families to form the subsequent generation (20 breeding pairs/line/generation, Appendix 2). We formed breeding pairs so as to minimize inbreeding using a computer resampling procedure (selection of the pairing pattern minimizing the median inbreeding coefficient). Assuming no inbreeding in F 1 , mean inbreeding rate was 9.6% (±1.9 SD) by F 7 . This corresponds to an average effective population size ("inbreeding effective numbers" sensu Crow & Kimura, 1970) of N e = 30.2.
Each generation, selection was performed on 636 fish on average (212 fish/line), and the selection procedure resulted in keeping on average 12% of individuals per line (number of breeders/ total number of fish before selection at 75 dph). We calculated the resultant selection differentials as the difference in maturity probability (i.e., proportion of mature fish) and standard body length after and before selection. Selection differentials across generations F 1 -F 6 for maturity probability and standard body length were +0.13 (0.12 SD) and +0.68 mm (0.18 mm SD) in the Control line, +0.10 (0.08 SD) and −1.06 mm (0.55 mm SD) in the Small line, and +0.13 (0.08 SD) and +2.05 mm (0.55 mm SD) in the Large line, respectively.

| Phenotyping
Eggs from each breeding pair were collected during a period corresponding to mother's 88-92 dph. Eggs were counted and photographed, and ImageJ was then used to measure their individual egg perimeters (9,795 eggs measured from F 1 to F 7 ). Hatched larvae were collected during a 5-day time window so as to synchronize hatching dates as much as possible. Birthdate was the median hatching date of each sibling family, and all siblings were thus assigned the same age.
At 0 (hatching), 15, 60, and 75 dph, each single individual was photographed, and then ImageJ was used to measure standard body length (from the tip of the snout to the base of the caudal fin, 16,808 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), which was a nondestructive proxy for the onset of maturity. All fish manipulations were performed after anesthesia in tricaine methanesulfonate (MS222), except at 0 and 15 dph when larvae and juveniles were manipulated with a Pasteur pipette and photographed in a droplet.

| Pituitary expression of candidate genes
An enzyme-linked immunosorbent assay (ELISA) is not available for medaka GH and ELISAs, in addition of being much less sensitive than reverse transcription quantitative real-time polymerase chain reaction (RT-qPCR), require plasma volumes that are too large to allow individual measurements in medaka. Hence, as a first approach to uncovering the molecular regulation of adaptive life-history evolution in medaka, we used RT-qPCR to measure mRNA levels of candidate genes in individual pituitaries. Specifically, we measured pituitary mRNA levels of β-subunits of gonadotropin hormones (LHβ and FSHβ) and GH. 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.
In each generation from F 1 to F 7 , 10 to 15 fish per line (233 fish in total) 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.
After sample homogenization by agitation (15 s vortexing), total RNA was 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 was 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. Medakaspecific 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 Appendix 1). 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 LightCycler 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 s, 60°C for 10 s, and 72°C for 5 s.
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). 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 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 practices, 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.

| Data analyses
The aim of our statistical analyses was to estimate and test for an overall effect of the selected lines on traits, pooling data from generations F 3 to F 7 and treating generation as a random effect. An archive containing datasets and scripts to reproduce analyses can be downloaded at https://doi.org/10.5061/dryad.0cfxp nw05.

| Standard body length at 75 dph
We modeled response to selection as the line effect on standard body length at 75 dph (Sdl75) of each individual i: effects of the generation (F 3 -F 7 ) treated as a normally distributed random effect, selected line (Small, Large, and Control), and sex (I, M, or F), respectively. Age is age in dph coded as a continuous variable, and Inb is individual inbreeding coefficient. Finally, 2 is residual variance and 2 G is the variance of the normally distributed generation effect.

| Probabilistic maturation reaction norms
We visualized the effect of anthropogenic selection on the maturation process using probabilistic maturation reaction norms (PMRNs). This approach was developed to account for the plastic effects of juvenile somatic growth rate on the maturation process, such that a shift in the maturation reaction norm may be interpreted as an evolutionary shift in maturation (Heino & Dieckmann, 2008;Heino et al., 2002). PMRNs classically account for the effects of age and body length on maturation, but they may also be "higher dimensional" to account for the effects of body mass or individual somatic growth rate (e.g., Morita & Fukuwaka, 2006). Here, however, we neither weighed individual medaka nor followed individual growth trajectories. Therefore, we used classical age-and length-dependent PMRNs, which have been demonstrated to be as efficient as higher-dimensional PMRNs to detect evolutionary trends (Dieckmann & Heino, 2007).
For each medaka line, we computed age-and length-dependent PMRNs, defined as the age-and length-dependent 50% probability for an immature medaka to initiate maturity (as informed by the onset of secondary sexual characteristics), using the methods of Barot, Heino, O'Brien, and Dieckmann (2004) and Van Dooren, Tully, and Ferrière (2005). The methods consisted first in computing maturity "ogives" as: where y i is the maturity status of an individual fish i (0 or 1), is the Bernoulli distribution of "success" (maturity) probability p, and In is the natural logarithm. Other subscripts or variables are as described above.
By letting the effects of both and on p varying for each selected line, this model captured potential effects of selection on both the intercept and slope of the PMRN.
Second, we computed the maturation probability m(a , s ) at each growth increment as described by Barot et al. (2004): 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 simulated slow, median, and fast growth curves (Harney, Van Dooren, Paterson, & Plaistow, 2013;Van Dooren et al., 2005).
Finally, we computed line-specific PMRNs as the age and length combination (a t , s t ) at witch maturation probability reached 50%, that is, as the age and length combination that satisfied the following condition (Harney et al., 2013;Van Dooren et al., 2005): We estimated m(a , s ) for 200 growth increments equally spread between ages 0 and 87 dph. Full propagation of error distribution for m(a , s ) was obtained by iterating the procedure for each

Markov chain Monte Carlo (MCMC) sample of the 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, is the Binomial distribution, Mo. inb is mother inbreeding coefficient, Fa. inb is father inbreeding coefficient, and is an overdispersion effect accounting for the fact that observed variance was larger than canonical variance of the Binomial distribution. 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.

| 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 (1) (3) breeding pair i were zero-inflated Poisson-distributed and modeled as: where 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 ( ).
is the natural logarithm of number of days during which eggs were collected (varied from 4 to 5 days), Parent. sdl is average parent standard body length, and is an overdispersion parameter.
Other subscripts or variables are as described above. In this model, the effect of parent 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
Individual egg perimeter i in mm was modeled as follows: where the variables are as described above in Equations 1 and 3.

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

| Larval size at hatch
We modeled individual standard body length at hatch Sdl0 as: where variables are described above.

| Hormonal profile
where W is the Wishart distribution, R is a scale matrix (diagonal matrix of dimension j), and = j denotes degrees of freedom (Lunn, Jackson, Best, Thomas, & Spiegelhalter, 2012). In practice, R, which is supplied as data, contains on the diagonal and 0s in nondiagonal entries.

| Natural selection
Our datasets also allowed us to measure natural selection, which often opposes the effects of artificial selection (e.g., Carlson et al., 2007). In medaka in the laboratory, natural selection may act on the standard body length of the selected parents through affecting their reproductive success or through the survival of their progeny. We visualized these potential effects of natural selection using quadratic regressions of fitness components (namely fecundity, hatch rate, number of progeny reaching age 75 dph, and number of progeny kept as breeders for the next generation) on mean parental body length in linear models (Lande & Arnold, 1983). Specifically, for (4) , fecundity, the number of progeny reaching age 75 dph, and number of progeny kept as breeders for the next generation we used zero-

| MCMC parameter estimation
All models were fitted using MCMC in JAGS (Plummer, 2003) through the jagsUI R package (Kellner, 2019) in R 3.6.3 (R Core Team, 2020). We used weakly informative priors and, 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).
We tested the significance of effects from posterior parameter distributions using a test equivalent to a two-way 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, Meng, & Stern, 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 (Appendix 4), that is, by setting one level of each factor as a reference levels as is done by default in the R software.

| RE SULTS
Effect sizes for response to selection of the 14 measured traits are presented in Table 1 Our second prediction was that evolution of body size should be paralleled by evolution of correlated traits, and in particular of age and size at maturation, size-specific fecundity, egg sizes, size at hatch and larval survival. Only maturity probability at 75 dph responded as expected, and more sharply so in the Large than in the Small line (Table 1) In the Small medaka line, maturity probability at an average age and body length did not respond to selection (Model 2 in Appendix 4) and, accordingly, PMRNs for the Small and Control lines largely overlapped ( Figure 2). Noticeably, however, there were some signs of an increased reproductive investment in the Small medaka line: The length-corrected maturity probability decreased less fast with an increasing age than in the Control line (Model 2 in Appendix 4), and egg sizes increased (  However, these trends were not statistically significant (Table 1

| D ISCUSS I ON
We measured in the laboratory the realized evolvability of body size in response to size-dependent selection in wild-caught medaka fish. We show that medaka responded to selection for a large body size, but not to selection for a small body size. Before discussing this unexpected result, we start with a mini review of previous harvestsimulating experiments and how their results and designs compare to ours.

| Laboratory harvesting experiments
Size-selection experiments are a classic in evolutionary biology, and have been conducted multiple times on model organisms such as mice (e.g., Falconer, 1973;Macarthur, 1949), chicken (Dunnington et al., 2013), or drosophila (e.g., Hillesheim & Stearns, 1991;Partridge et al., 1999). More recently, problems with overexploitation have renewed the interest in size-selective experiments mimicking sizeselective harvesting. In a pioneering study, Edley and Law (1988) have  TA B L E 1 Effect sizes of bidirectional selection on body size on phenotypic, life-history, and neuroendocrine traits in medaka menidia maintained in 700 L tanks (about 100 breeders/generation/ population). The Atlantic silverside is an annual fish, and it was assumed that all individuals were mature at selection such that selection was imposed on body size only. Conover and Munch (2002) found that the mean weight of fish evolved in the expected direction and, by generation F 5 , an average fish aged 190 dpf weighted 3.5 g in the Control lines, 2.5 g in the Small lines, and 4.5 g in the Large lines.
These differences were due to differences in somatic growth rate and underlying traits (Walsh et al., 2006). (2012)  length but matured at a later age (but not size), while the Small line showed no change in juvenile somatic growth rate but evolved lower asymptotic length and maturation at a smaller size (but not age).

Amaral and Johnston
All the above-listed designs, and ours as well, imposed truncation selection on body size, which may or may not accurately reproduce the form of fishing-induced selection depending on the fishing gear.
Towed gears and long-lining catch all individuals above a threshold body size, and their effects are thus accurately simulated by truncation selection. In contrast, gillnets or traps selectively target medium-long individuals Kendall et al., 2009;Kuparinen et al., 2009;Lagler, 1968;Millar & Fryer, 1999), and thus generate at the same time disruptive selection and directional selection against a large body size Edeline et al., 2009) Another key feature of all previous laboratory harvesting experiments is that they used a mass-selection design with replication of the selected lines, but no control over effective population sizes, inbreeding rate or natural selection. To avoid these problems, we isolated selected pairs and raised their offspring in individual tanks, keeping track of the pedigrees along the experiment. This made it possible to control for the number of offspring per individuals, to maximize effective population sizes, to limit inbreeding throughout the selection procedure, and to measure natural selection. To our knowledge, this is the first time that such a high level of control is achieved in a size-selection experiment on fish.
However, because the number of individuals included in such an experiment is limited, line replication trades off with increasing effective population size N e . A large N e decreases genetic drift, limits the effect of linkage disequilibrium on selection limits, and delays the unavoidable increase in inbreeding (Hill & Robertson, 1966;Robertson, 1960); see, for example, Weber and Diggins (1990) for experimental evidence.
In particular, avoiding genetic drift and inbreeding is crucial when studying the evolution of correlated characters (Phillips, Whitlock, & Fowler, 2001). Therefore, we chose to derive three large-population lines (25 < N e < 30 in each) rather than replicating small-population treatments. A pedigree-based quantitative genetic analysis suggests that medaka trait dynamics in the Large line were not compatible with F I G U R E 3 Medaka endocrine response to bidirectional selection on body size. Pituitary mRNA levels for A: the luteinizing hormone (LH, β subunit), B: the follicle-stimulating hormone (FSH, β subunit), and C: growth hormone (GH) were standardized by actin β (ACT) levels and log-transformed. Dots represent raw data.

| Medaka phenotypic and life-history responses to bidirectional selection on body size
At the end of our experiment ( mating the fish 14 days after that 50% of the population reached maturity, a delay that was possibly not long enough to allow for 100% of the fish to reach maturity, in which case selection was applied both on body size and for maturity (similar to our own design). As discussed by Le Rouzic et al. (in press), available evidence suggests that response to such bivariate selection on correlated traits is often erratic.
In our experiment, lack of body-size response to selection in the Small medaka line could not be ascribed to an absence of artificial selection, which was strong and consistent, nor due to the counteracting effects of natural selection, which remained weak compared to the strength of artificial selection, nor due to inbreeding which was by F 7 identical among the random-and large-harvested lines. Instead, the absence of evolution in the Small medaka line suggests that medaka are at a lower evolutionary limit for body size.

| Medaka neuroendocrine response to bidirectional selection on body size
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 that stimulates not only somatic growth rate (Canosa, Chang, & Peter, 2007;Reinecke et al., 2005) but also maturation, and also mediates osmoregulation and the stress response (Le Gac et al., 1993;Rousseau & Dufour, 2007;Wendelaar Bonga, 1997).
We expected pituitary mRNA GH levels to be altered in parallel  (Gomez et al., 1999). Finally, the positive LH-GH correlation significantly increased in the Large medaka line, indicating that size-dependent selection may alter patterns of hormonal synergies. Future transcriptomic approaches on central and peripheral tissues will maybe provide a deeper understanding of the molecular regulation of response to size-dependent selection in the medaka.

| CON CLUS IONS
Inability of medaka to respond to selection for a smaller body size is a warning signal that calls for increasing research efforts to assess life-history evolvability in wild populations. 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 such as body size and maturity (Lande & Arnold, 1983;Le Rouzic et al., in press), 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 size-dependent, anthropogenic selection (Therkildsen et al., 2019).
In the future, comprehensive approaches melting wide-spectrum candidate genes, transcriptomics, and genome scans of experimentally and wild-selected populations will probably be needed to finely decipher the molecular architectures that regulate the adaptive evolution of life histories and that ultimately support the maintenance of biodiversity and ecosystem productivity.

ACK N OWLED G M ENTS
We are grateful to Prof. Kiyoshi Naruse (NIBB, Okazaki, Japan) for his support in obtaining and maintaining medaka from Kiyosu. 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.

E TH I C A L A PPROVA L
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). The committee also confirmed that our methods were performed in accordance with the relevant guidelines and regulations on animal research.

A PPE N D I X 2
Schematic diagram of the experimental design. Curves show medaka standard body length distributions at 75 days posthatch (dph). Shaded areas represent mature individuals and colored areas the mature individuals that were kept as breeders to form the next generation. Blue cubes represent the 3 L aquariums in which medaka were maintained. Summary of MCMC parameter estimates for models 1-8. We used an "effect" model parametrization in which the intercept is the mean value of the response variable (on the link scale) in the reference level of factor predictors. ΔX are estimates for the difference between mean value of the response variable in the factor level X and model intercept.