Metapopulation structure modulates sexual antagonism

Abstract Despite the far‐reaching evolutionary implications of sexual conflict, the effects of metapopulation structure, when populations are subdivided into several demes connected to some degree by migration, on sexual conflict dynamics are unknown. Here, we used experimental evolution in an insect model system, the seed beetle Callosobruchus maculatus, to assess the independent and interacting effects of selection histories associated with mating system (monogamy vs. polygamy) and population subdivision on sexual conflict evolution. We confirm traditional predictions from sexual conflict theory by revealing increased resistance to male harm in females from populations with a history of intense sexual selection (polygamous populations) compared to females from populations with a history of relaxed sexual selection (monogamous populations). However, selection arising from metapopulation structure reversed the classic pattern of sexually antagonistic coevolution and led to reduced resistance in females from polygamous populations. These results underscore that population spatial structure moderates sexual selection and sexual conflict, and more broadly, that the evolution of sexual conflict is contingent on ecological context. The findings also have implications for population dynamics, conservation biology, and biological control.


Study system
The seed beetle Callosobruchus maculatus is a widely distributed pest of stored legumes. Adults are facultatively aphagous (they do not need to feed or drink) and the species is a capital breeder since the beetles obtain all of the resources needed for survival and reproduction during the larval stage (see below and Messina and Fry 2003). Mated females cement their eggs on the surface of the host beans. Upon hatching the larvae enter the seed, where they complete their development and pupate, before emerging as adults. The species has a short life cycle (around 22 days at 30•C) and high fecundity (a female typically produces several scores of eggs over the course of a few days). The species is a model system in the study of sexual conflict (Crudgington and Siva-Jothy 2000;Arnqvist et al. 2005;Gay et al. 2010;Cayetano et al. 2011;Zuk et al. 2014;Berger et al. 2016;Dougherty and Simmons 2017;Iglesias-Carrasco et al. 2018a;Sayadi et al. 2019;McNamara et al. 2020). The species is highly polygamous and the male intromittent organ is densely coated with thick spines. These spines determine to a large extent male reproductive success and reduce female fitness (Crudgington and Siva-Jothy 2000;Hotzy and Arnqvist 2009;Gay et al. 2010;Hotzy et al. 2012). We sourced the beetles from an outbred population (South Indian stock population) that presents substantial phenotypic and genetic variance (Fox et al. 2004;Bilde et al. 2008;Berg and Maklakov 2012;Rodriguez-Exposito 2018;Zajitschek et al. 2018;Canal et al. 2021). Beetles were kept in walk-in climate chambers (Fitoclima 10000 EHF, Aralab) at a constant 29 °C temperature with 40% humidity and a 12 h/12 h light/dark cycle, and were cultured in organic mung beans (Vigna radiata). Further details about the culturing of our stock population can be found elsewhere (Zajitschek et al. 2018;Canal et al. 2021). Importantly, the standard culturing conditions of C. maculatus resemble the conditions to which these animals have adapted for thousands of years (infestation of dry legume seed storages), and the conditions that natural populations experience nowadays, which makes this species a suitable laboratory model system (Berger et al. 2014).

Additional details on the propagation of selection lines
The selection of inoculated beans from each line to generate the next generation of animals was carried out randomly but to avoid unintended selection we followed a systematized plan for the propagation of the lines. This plan also ensured that we imposed the appropriate levels of sexual interactions envisaged by the experimental design: a) Non-structure and polygamy (NSPoly) lines: The 25 breeding couples per line were placed in a 750 ml plastic container with approx. 1600 beans (64 seeds per female). As for all selection lines, sexual interactions and egg laying took place during two days (see main text). For each line, the 50 breeders for the next generation were randomly selected from the virgin adults emerging from the 150 inoculated beans randomly collected from the population container.
b) Structure and polygamy (SPoly) lines: For each of the 5 subpopulations within each line we placed five breeding couples in a 140 ml plastic container with approx. 320 beans (64 seeds per female). Thirty inoculated beans per subpopulation (i.e., 150 per line) were randomly selected and individually kept. A random selection of 4 virgin males and 4 virgin females emerged from these beans in each subpopulation plus one virgin adult of each sex emerging from the beans in a different subpopulation (imposed migration, see above and Fig.  1) were designated the subpopulation breeders for the next generation. Thus, within each line metapopulation structure was mimicked (with gene flow among subpopulations) and individuals could engage in polygamous sexual interactions. c) Non-structure and monogamy (NSMono) lines: The 25 breeding couples per line were each placed in a 30 ml plastic container with 64 beans (i.e., 64 beans per female). Six inoculated beans per container (i.e., 150 beans per line, each bean containing a single egg) were randomly selected. The breeders for the next generation were then randomly sourced from the virgin adults emerging from the pool of 150 beans collected in each line. Thus, within each of these lines, monogamous matings were enforced and there was no population subdivision whatsoever (breeders originated from the pool of 150 beans). d) Structure and monogamy (SMono lines). In this treatment, as in SPoly lines, there were 5 breeding pairs for each of the 5 subpopulations in each line. However, monogamy was enforced because each male-female pair was housed in a 30 ml plastic container (with 64 beans each). Six inoculated beans per container were randomly collected, as in the NSMono lines; however, to impose the effects of spatial population structure, inoculated beans were not pooled at the population level (as in NSMono lines) but at the subpopulation level (n = 30 inoculated beans per subpopulation, n beans per line = 150). Like SPoly lines, a random selection of 4 virgin males and 4 virgin females emerged from the beans in each subpopulation plus one virgin adult of each sex emerging from the beans in a different subpopulation (imposed migration, see above and Fig. 1) were designated the subpopulation breeders for the next generation. Therefore, within each of these lines, metapopulation structure was mimicked (with gene flow among subpopulations) and the scope for sexual selection and sexual conflict was greatly reduced (individuals could only engage in enforced monogamous sexual interactions).

Generation and use of tester individuals
Each of the near-isogenic lines was started with one founder male-female pair from the stock population and the line was propagated through a full brother-full sister mating protocol that spanned many generations (>15, see specific details about the number of generations in each case below). In some cases the assays were carried out using tester individuals from such pure near-isogenic lines. In other cases, tester individuals were obtained after crossing two nearisogenic lines. These crosses were carried out to generate a standardized but heterozygous genetic background (see for instance Garcia-Gonzalez and Dowling 2015). We refer to the lines from these crosses as Standardized Heterozygous Lines (SHL). Crosses leading to SHL lines were established by housing together in a 750 ml plastic container 25 males and 25 females collected from two randomly chosen near-isogenic lines (12 beetles from each sex from one of the lines, and 13 beetles from each sex from the other line) two generations before carrying out the assays. These lines were kept with a non-overlapping generations protocol in which within-line matings (25 males and 25 females) were allowed. Likewise, when pure nearisogenic lines were used, two generations before the assays 25 pairs from each line were transferred to 750 ml containers and within-line matings allowed under a non-overlapping generations protocol. This punctual "mass-breeding" within-line protocol ensured the production of enough tester animals for the assays. Near-isogenic lines or SHLs were used depending on logistics, to obtain enough tester animals for the assays, but pure near isogeniclines and SHLs were never mixed. We thus standardized the interacting genetic background with which experimentally evolved individuals were measured. Since the relevant comparisons in our tests are those carried out across selection regimes within each specific assay, and for any given assay the tester strain used was identical across selection regimes, the type of sourcing for tester animals does not affect the conclusions in the study. Details of the tester individuals used in the different experiments are as follows: -Medium levels of sexual conflict assay (female longevity under variable female mating rates): Variation attributable to mate identity in female mating rates was minimized by using tester males that belonged to a same near-isogenic line or a SHL in each mating opportunity across all females. All females across all selection lines were exposed to the same line ID of tester males each mating opportunity, and the ID of the near-isogenic line or SHL used to generate the tester males were rotated across mating opportunities. Near-isogenic lines used had been cultured on a full sib mating protocol for 15 and 33 generations, for the assays carried out at generations 12 and 30, respectively.
-Low levels of sexual conflict assay (female lifetime reproductive success -LRS-and longevity after single mating): Tester individuals were drawn from a near-isogenic line generated after 35 generations of full-sib matings.
-Extreme levels of sexual conflict assay (female LRS and longevity under continuous exposure to intense male harassment): Tester individuals used in this experiment were collected from a near-isogenic line generated after 50 generations of full-sib matings.

Body size measurements
Elytron length was used as a proxy of body size (Eady et al. 2004;Rodriguez-Exposito 2018). Elytron length was measured from photographs of frozen beetles. Specimens were placed on top of a pedestal made of Blu Tack (Bostik, US)( Fig. A9) and images were taken using a stereomicroscope SteREO Discovery V8 (Carl Zeiss Microscopy GmbH, Germany) connected to an AxioCam Icc 1 camera (Carl Zeiss, Germany). The maximum length of the right elytron was measured on the images (Fig. A9) using the software ZEN 2 (blue edition, Carl Zeiss, Germany, 2011). Body size was measured on all male or female focal individuals. The size of the tester individuals was also measured in the single mating assays. The repeatability of body size measures, as calculated from subsets of individuals that were measured twice, is very high. All the experimental females in the female mating frequency assays were measured twice, and the repeatability of elytron length, calculated following Becker (1984), was 0.99 and 0.98 for generations 12 (n = 112 females) and 30 (n = 160 females) (all p < 0.0001). In the assays of sexual conflict following single mating the elytron length of 319 males and 317 females was measured twice and repeatabilities were again high (males: 0.99, females: 0.93; all p-values <0.0001). The repeatability of elytron length was additionally measured in a set of 320 males and 320 females at generation 34; repeatabilities were very high regardless of sex (all R> 0.99; all p < 0.0001). In light of these results, the rest of individuals were only measured once, and for those individuals with two measures of body size we used the mean value of these two measurements in the analyses.

Effective population sizes (Ne)
Under our experimental protocol, slight variation in effective population sizes due to spatial structuring may arise due to sampling. Larger variation in Ne is, however, expected between polygamous and monogamous populations due to a series of reasons including sexual selection. Apart from the fact that in polygamous populations genetic hitchhiking of genetic polymorphisms in linkage disequilibrium with sexually selected alleles can lead indirectly to reduced Ne and genetic diversity (Snook et al. 2009), variance in male reproductive success due to variation in mating and fertilization success under sexual selection may reduce Ne in these populations compared to monogamous lines. In contrast, a higher reduction in Ne would be expected to arise in monogamous populations if infertile mating rates, due to either male or female infertility, were common. Infertile matings would reduce Ne in monogamous populations to a larger extent than in polygamous populations because polyandry would buffer the outcome of male infertility on female productivity (Garcia-Gonzalez 2004;Yasui and Garcia-Gonzalez 2016). For instance, a monogamous pair in which the male is infertile results in that neither the male nor the female leaves offspring, while under polyandry only high rates of infertility among males would result in the females having null reproductive success because the chances that all mates of a polyandrous female are all infertile would be low (Garcia-Gonzalez 2004). In other words, the more mates a female mates with, the less likely it is for the female to have complete reproductive failure (Yasui and Garcia-Gonzalez 2016). Our lines do not differ in the extent of infertility and the average infertile matings rate is 2.54% (Rodriguez-Exposito and Garcia-Gonzalez, unpublished; Rodriguez-Exposito 2018). This means that, on average, for monogamous lines, less than one pair of beetles out of the 25 pairs in each population would produce no offspring (probability that a monogamous female produces no offspring = 0.025; 0.63 females out of the 25 in each line). Assuming that all infertile matings are the result of male infertility, the probability that any polyandrous female produces no offspring would be, however, 0.00064 (i.e., only 0.16 females out of the 25 in each line), if she mates with two males, or 0.000016 (0.0004 females out of the 25 in each line), if she mates with three males. Assuming that male reproductive success is random the effective population size can be calculated as: where Nf is the number of breeding females and Nm is the number of breeding males (Crow and Kimura 1970;Falconer and Mackay 1996;Snook et al. 2009;Walsh and Lynch 2018). Incorporating female multiple mating and variance in male reproductive success, the effective population size can be calculated as (Rice and Holland 2005;Snook et al. 2009): where, Pi is the proportion of offspring sired by male P (Rice and Holland 2005). However, this formula, which has been widely used across studies, applies to situations in which each female is mated to a different set of males. Our case is more complex because in the polygamous populations each female can mate with different males but each male can also mate with different females. Furthermore, in some populations we have structure in regards to mating interactions. In our case, therefore, we recurred to Monte Carlo simulations to approximate effective population sizes. We used Poptools 3.06 (Hood 2008) to mimic the conditions in each selection regime, and to draw estimates of Ne for each selection regime. For each selection regime and each combination of conditions (scenarios) simulated, we calculated the number of sires and dams represented in a given generation. We then calculated Ne summing up both the number of sires and dams. The number of iterations within each scenario was set to 10000. The average and 95% CIs for Ne estimations for all the scenarios within each selection regime are shown below. The files with the simulation structure and sampling protocol from Poptools can be found on Dryad at https://doi.org/10.5061/dryad.r2280gbd9. Simulations needed to be adjusted for each selection regime to mimic the precise extraction protocol to obtain the breeding individuals followed in our selection experiment: a) NSPoly lines: Virgin C. maculatus females typically mate within seconds or minutes after encountering a male, and many singly-mated females that are given the opportunity to remate with a different male 24h after their first mating do so (Rodriguez-Exposito 2018). Accordingly, and given that females in NSPoly lines are continuously housed with plenty of males, in the simulations we conservatively assumed that over the course of 48 hours females would mate with at least two males. We simulated two scenarios here: double and triple matings. Our estimate of Ne is probably conservative as it is likely that mating rates are higher in the polygamous populations. The distribution of paternity for each male was based on the empirical estimate that last male sperm precedence when females mate to two or three males averages 0.80 in the species (Eady 1991;Eady 1994;Eady and Tubman 1996;Simmons 2001;Vasudeva et al. 2014;McNamara et al. 2016). Fecundity is kept constant at 40 eggs in these lines and all lines, based on previous findings regarding number of eggs produced the first days after mating (Zajitschek et al. 2018). In any case, we checked that reasonable variation around fecundity did not introduce substantial variation in the results (see below). A key factor that has been overlooked in previous studies but that we implement in our calculations is that mating interval matters for the calculation of Ne. This is so because even in a system with strong last male sperm precedence females would produce offspring of first male(s) in between matings. The production of offspring before all the known (or estimated) matings have been completed would affect Ne calculations to a great extent. We took into account this factor, and obtained estimates with or without correction for fecundity between matings. In absence of the correction the (unrealistic) assumption is that all matings take place at the same time. When the correction was applied we assumed that the first mating of females occurred at time 0 (virgins females typically mate within a few seconds or minutes after sexual encounter, as noted), and that each subsequent mating occurred 16 hours after the previous mating (i.e., for a double mating scenario the two matings would take place at time 0 and 16 h, and for a triple mating scenario the matings would take place at times 0, 16 and 32 h). Variation in mating interval does not introduce large variation in Ne (results not shown). For each case we calculated paternity shares according to the sperm precedence patterns in the system (see files with the simulation structure and sampling protocol at https://doi.org/10.5061/dryad.r2280gbd9). For instance, if a female mates with male A at time 0 h, male B at time 16 h, and male C at time 32 h, and assuming, for simplicity, constant oviposition rate over the first 48 hours (i.e, 0.83 eggs/hour for a total fecundity of 40 eggs), the distribution of paternity at the end of the 48 hours is calculated as follows: during the first 16 h (i.e., before the female's second mating) the female will lay 13 eggs (0.83 eggs per hour during 16 hours) from 1st male. For the remaining 27 eggs, 13 eggs will be laid during the next 16 hours (after the second mating and before the third mating). These 13 eggs would be sired according to Pn (paternity of the last male to mate) = 0.8; that is, 10 eggs would be sired by the second male and 3 eggs would be sired by the 1st male. The remaining eggs (14) would be laid after the female has mated with the three males, and again, they would be sired according to the sperm precedence pattern in the system. So, 11 eggs would be sired by the third male, 2 eggs from the second male, and 1 egg would be sired by the first male. In sum, the distribution of paternity for the first, second and third males, would be 17, 12, and 11 offspring, respectively. This illustrates that considering mating interval and fecundity in between matings is important. Failure to do so yields a quite different distribution of paternity. The distribution of paternity under the assumption that the eggs of a female are sired by three males according to Pn = 0.8 without considering that mating interval affects fecundity is 2, 6, and 32 offspring for the first, second and third males to mate with the female, respectively. The correction, however, does not imply a substantial change in Ne in our case (see Appendix Results). We also implemented in the simulations scenarios were infertile matings were considered and confirmed that their effects were trivial upon Ne variation (see the section: Extended results on effective population sizes). b) SPoly lines: Simulations here proceed as described above with the important exception that individuals are structured into 5 populations and that there is a 20% migration rate among subpopulations. We implemented in the simulations these features, and mimicked the empirical breeding protocol: thirty offspring are randomly extracted from each subpopulation, from which a random selection of 4 males and 4 females, plus another couple of adults randomly selected from a different subpopulation (migration) are the breeders at the subpopulation level for the next generation. This protocol is set for each of the 5 subpopulations within each line and the number of effective breeders at the line level is calculated. With these numbers, as above, we estimated Ne summing up the number of effective breeders, or using equation 1.
Ne in spatially structured populations can be also calculated following Wright's (1951) island model, where a metapopulation consists of d demes, each containing N individuals, and each deme contributes a fraction of its genes m (equal across demes) to the migrant pool, which is then distributed among the remaining demes (Walsh and Lynch 2018: Pp. 72-73): where Nd is the sum of the demic effective sizes. In our case, the structured populations consist of 5 demes, each containing a fixed number of individuals (10), and each deme contributes an equal fraction of its genes (0.2%, as 2 out of the 10 individuals in each deme migrate to a different deme) to the migrant pool. We can therefore calculate the effective population size in each deme following equation 3 as if there was no variance in male reproductive success: d = 5, Ne of each deme = 10, Nd = 50, m = 0.2, Metapopulation Ne = 54. This estimate is, however, unrealistic as it does not take into account variance in male reproductive success. We applied simulations in which we implemented variance in male reproductive success (as above, with Pn = 0.8) and calculated the sum of the demic effective sizes (Nd) in absence of migration to then apply equation 3. Using this method we estimated Ne with and without considering mating intervals and subsequent distributions of paternity (see above), and for the case of two or three mating partners per female, as above. c) NSMono lines: In these lines each monogamous female was mated to a unique monogamous male and each pair equally contributed (with the exception of occasional infertility) with 6 offspring to the pool of offspring (n = 150) that was randomly sampled to begin the next generation. Such random extraction from the pool of 150 offspring to seed the next generation inevitably means a slight reduction in Ne values compared to a situation where exactly two individuals are taken from each couple. Thus we also applied to these lines simulations to estimate Ne following the precise protocol we followed in our selection experiment. A scenario considering the effects of infertile matings (1 infertile male out of the 25 in each population; see above) was also inspected.

d) SMono lines:
In the simulations to ascertain Ne in SMono lines we introduced, as for the SPoly lines, spatial structure according to the empirical protocol. Here each monogamous male-female pair contributed with 6 offspring to a pool of 30 offspring per subpopulation, from which the breeders (n = 5 pairs of breeders per subpopulation) for the next generation were randomly sampled (taking into account that 1 pair of breeders migrates from one population to another). This is simulated for each of the 5 subpopulations within each line and the number of effective breeders at the level of the line is calculated. As for SPoly lines, we also calculated Ne using the island model formula (equation 3), and for the different combination of conditions (variation in mating frequency, infertility and oviposition in between matings).

Full details for statistical analyses Modelling approach:
In the analysis of fitness traits after single mating or continuous male-female cohabitation we had data on the evolution of male harm (fitness of tester females mated to focal males), and data on the evolution of female resistance (fitness of focal females mated to tester males). The response variables in this group of assays were female longevity and LRS. We also measured baseline longevity data in virgin individuals from the different selection lines. We run a series of complementary analytical approaches on the response variables. First, since the variables are count variables (days alive, number of offspring) we run analyses using Generalized Linear Mixed Models (GLMMs) with a Poisson error distribution. These models were fitted using the function glmer in the package lme4 (Bates et al. 2015). To account for issues of data dispersion in our GLMM we applied quasi-likelihood by correcting the standard errors of the estimates by the dispersion factor and then recomputing Z-and p-values accordingly. Second, complementary and confirmatory Linear Mixed Models (LMMs) were fitted on the dependent variables using the function lmer in the package lme4 (Bates et al. 2015). To improve the validation of the LMM models we squared LRS. Both LMMs and GLMMs yielded similar results in general, and we show here the results pertaining to LMMs because the validation of these models was superior.
Fixed effects structure of models: In general, as predictors, the models included the mating system experimental evolution treatment (monogamy/polygamy), metapopulation experimental evolution treatment (yes/no), and the interaction of these two factors, as fixed effects. Body size data pertaining to the individuals in the analyses were included as covariates in all the analyses, and in the analyses of responses in male harm and female resistance we also included the second order interactions between male or female body size and treatments because male body size might be associated with male harm and female body size with female resistance (see results and Pitnick and Garcia-Gonzalez 2002), and the influence of body size on harm/resistance could be expected to differ according to treatment. For instance, larger male body size should be associated with increased harm to females but only for males from polyandrous lines; larger male body size in monandrous males could, however, expected to increase female longevity if larger males transfer larger ejaculates (Czesak and Fox 2003;Iglesias-Carrasco et al. 2018b). To account for reproduction-survival relationships (see text), female longevity was included as a covariate in the analyses of LRS, while LRS was included in the models where female longevity was the response variable. In the data from continuous male-female cohabitation there was slight variation in the age that tester females entered the male harm assays (mean ± SE age of females: 1.75 ± 0.07, range 1-5 days old), and thus female age was also entered, as a control covariate, in these models. Female age was less variable in the assays where females from the selection lines were cohabiting with tester males (female resistance assays: mean ± SE age of females: 1.05 ± 0.02, range 1-2 days old), but age was also included in the models. In all models all the covariates were mean-centered (Schielzeth 2010).
In the assays in which experimental females were given controlled opportunities for remating only female longevity (number of days alive after entering the assay) was measured as in these assays females were not allowed to lay eggs. All females mated at least once and no female died before completing the mating opportunities period. Tough mating frequency does not differ among selection regimes at the time of the test (Rodriguez-Exposito and Garcia-Gonzalez, unpublished; Rodriguez-Exposito 2018), female mating frequency was included as a control covariate in the model (mean mating frequency = 1.94, SE = 0.04, range 1-4, n = 262, 73% of females mated more than once). There was very little variation in female age at the start of the assay (mean ± SE = 1.48 ± 0.04, range 1-3, n = 262), but age was nonetheless entered as another control covariate in the models, as well as generation. Predictors in these models included the two experimental evolution treatments, female body size, and the second order interactions involving these three predictors.
Random effects structure of models: Line ID (unique code for each of the 16 selection lines) was included as a random effect in all models. Models included covariates for which interactions with the experimental evolution treatments were deemed of interest (e.g., see discussion on body size, or on the relationship between LRS and longevity above) and so we fitted by-line ID random intercept and random slopes models to avoid inflation of type I error (Schielzeth and Forstmeier 2009;Barr et al. 2013;Bates et al. 2015;Arnqvist 2020). The correlations between intercept and slopes were also included in the models (Barr et al. 2013). Nevertheless, in complex models (e.g., including 2 fixed factors plus three covariates) there were occasional issues of convergence. In such cases, the models were run without including the correlation between random intercept and random slopes, to simplify the random effects structure of the model. This decision was based on the results of Barr et al. (2013) informing that when maximal linear mixed effects models fail to convergence, removing random correlations does not affect the anti-conservativity of the model. Models excluding random correlations are indicated in the table legends when applicable. For covariates that were included in the model as control predictors (e.g., age of female when entering the assay) random slopes were not modelled, as interactions between these control predictors and the treatments were not present in the model (they were not the focus of the analyses or justified by the data), and this procedure is expected not to inflate type I errors (Barr et al. 2013).

Significance of effects and additional considerations:
The function Anova (car package) was employed to evaluate the statistical significance of fixed effects, using Type II (Type III in the presence of significant interactions) Wald Chi-square tests. Significance of effects was carried out on maximum likelihood models, while parameter estimates were calculated using Restricted Maximum Likelihood, as it has been recommended (Zuur et al. 2009). Significance of the fixed effects in the GLMMs was calculated using the z distribution, and in the case of quasi-likelihood estimation correcting the standard errors of the coefficients as noted above. Visual inspection of diagnostic plots (qqplots and plots of the distribution of the residuals against fitted values) was checked to validate the models. We report mean ± standard error values throughout. Final sample sizes are indicated in Table A1, and Table A5 shows the means for the response variables for the different selection regimes across experiments. All the analyses were run in R 3.5.1 (R Core Team 2018).

Extended results on effective population sizes (Ne)
Effective population sizes calculated using different methods (see above) ranged from around 41 to 51 across all selection regimes (Table A2). In the calculations of Ne in polygamous populations, implementing variation in the number of mating partners or accounting for fecundity between mating intervals has implications for effective population size, but the magnitude of the differences across scenarios is small (Table A3). Likewise, the effect of the observed rates of infertile matings on the estimated Ne is largely negligible (Table A4). It has to be noted that Ne estimations are not far from the number of breeding individuals. This is partly the result of three features of the study system and experimental procedures leading, in combination, to low variance in reproductive success (family size): i) infertility rates are extremely low in the system, with over 97% of singly mated females laying eggs (Rodriguez-Exposito and Garcia-Gonzalez, unpublished; Rodriguez-Exposito 2018), ii) two days for oviposition allow individual females to typically lay dozens of eggs (Zajitschek et al. 2018), and iii) random sampling of offspring in the selection experiment (and replicated in the calculations of Ne), ensuring even representation across families.

Baseline longevity
Neither the mating system, nor the spatial structure or their interaction had a significant effect on the baseline longevity of either sex (Table A6, and see Maklakov et al. 2007). Longevity was significantly determined by body size in both sexes (Table A6, Fig. A4), with larger individuals living for longer regardless of the selection regime under which they had evolved.

Extended results from the extreme levels of sexual conflict assay
We measured female longevity and LRS under conditions of extreme intensity of sexual interactions, where each female was exposed to continuous lifelong male harassment and mating attempts from three males. In tests of the evolution of male harm (effects of male selection history on tester female fitness), there was a significant interaction between sexual selection and population subdivision selection histories on female LRS (see main text). In addition, female longevity was positively influenced by female body size (Table 4). The age of the female when entering the assay was also positively related to longevity, supporting the notion that females paid a reproductive cost when cohabiting with males (the older the female when entering the assay the longer she lived). Unlike the single mating situation, longevity and LRS were positively correlated. This is unsurprising as reversal of the longevity-offspring production trade-off is known to be context dependent in C. maculatus (Messina and Slade 1999;Messina and Fry 2003), and in this case, this reversal likely responds to the benefits that females obtain from multiple ejaculates (Savalli and Fox 1998;Zajitschek et al. 2018), including hydration benefits in this facultative aphagous capital breeder (Fox 1993;Arnqvist et al. 2005;Edvardsson 2007;Ursprung et al. 2009;Iglesias-Carrasco et al. 2018a). This survivalreproduction trade-off reversal, which is consistent across assays of evolution of male harm and female resistance under the conditions of male-biased sexual cohabitation, indicates that trade-offs between life-history traits in sexual conflict systems depends on a complex balance between the benefits and costs of sexual interactions.
The analysis of female resistance (effects of standardised males on focal female fitness) assessed after continuous sexual cohabitation showed a significant positive association between female age and LRS on longevity, and between female body size and longevity on LRS (Table 5). There were significant interactions between the experimental evolution treatments and female LRS on longevity (and between the mating system treatment and longevity on female LRS), and between metapopulation treatment and female body size on longevity (Table  5). We, however, refrain from drawing too many conclusions from these interactions, and we are in general cautious when interpreting the outcomes of our assays under conditions of continuous male-biased sexual cohabitation. Such continuous female exposure to male harassment and sexual interactions may impose extreme levels of sexual conflict that can be unrealistic in natural populations because of three reasons: extreme male-biased sex ratio, continuous lifelong sexual cohabitation, and highly reduced possibilities for females to escape male harassment (Zajitschek et al. 2018). The housing conditions in the continuous cohabitation assays, with three males per female in a confined space, may have actually compromised copulations because of the high levels of disruptions of mating attempts by competitor males under these conditions. It was not logistically possible to monitor mating rates in the continuous cohabitation tests and thus female mating frequency cannot be accounted for in the analyses. These assays may be therefore both less informative about the consequences of sexual interactions than the tests where mating frequency was taken into account (low and medium levels of sexual conflict assays). Tables   Table A1. Initial and final sample sizes relative to the number of individuals assayed in the analyses of longevity and lifetime reproductive success (LRS) across experiments. The missing data column indicates the number of individuals removed from the analyses due to sample deterioration or missing data. Data for females that did not produce any offspring due to absence of mating or LRS equal to zero were excluded (see text). Note that in the statistical analyses the unit of replication is the selection line (total number of replicated populations = 16).  Table A2.

Experiment/assays Response
Average effective population sizes and CIs (in brackets) in the different selection regimes, calculated following simulations (see text). The first row of numbers within each selection regime (combination of metapopulation structure by mating system experimental evolution treatments) indicates the Ne estimate calculated summing the number of effective breeders. The second and third rows indicate Ne calculated using the number of effective breeders (always after accounting for variation in male reproductive success) to feed equation 1 or equation 3 (island model), respectively. The estimates for polygamous populations are for the case of 3 mates per female and correcting paternity after accounting for fecundity between mating intervals (for variations see Table A5).
No metapopulation structure  Table A3.
Ne variation in polygamous populations based on number of mating partners (2 or 3) per female, and on the correction for oviposition between matings (i.e., accounting for fecundity during mating intervals). Note that the more realistic situation is where a correction for fecundity between matings is simulated, as indicated in the text. The first row of numbers within each selection regime (combination of metapopulation structure by mating system experimental evolution treatments) indicates the Ne estimate calculated summing the number of effective breeders. The second and third rows indicate Ne calculated using the number of effective breeders (always after accounting for variation in male reproductive success) to feed equation 1 or equation 3 (island model), respectively. Average and 95% CIs (in brackets) Ne estimations are shown. NSPoly: Non-structure and polygamy; SPoly: Structure and polygamy. Average effective population sizes and CIs (in brackets) are given.  NSPoly: Non-structure and polygamy; SPoly: Structure and polygamy. Note that females did not have access to oviposition substrate in the medium levels of sexual conflict assay and therefore longevity of females in this experiment is similar to longevity of virgin females (no sexual conflict assay) but it is sensibly higher than that in the other assays (low and extreme levels of sexual conflict) in which oviposition was allowed (Messina and Slade 1999;Messina and Fry 2003) Figure A1.

Appendix Figures
Outline of the experimental design to measure male harm and female resistance after single mating. Each virgin focal individual (in black) was paired with a single tester individual (in blue) of the opposite sex. Afterwards, the male was removed and the female was transferred to a small plastic container (30 ml) with beans to allow oviposition. The oviposition container was checked daily to measure female longevity. Adult offspring emerged up to day 29th after the day the female mated provided a measure for LRS (see text).

Figure A2.
Outline of the assays to measure female longevity when females were given multiple but controlled opportunities for mating/remating. Each female was placed daily (12 days for females from generation 12, 10 days for females from generation 30) with a standardized tester male (see text) in a glass vial for 15 minutes. At the end of that period the tester male was removed and the female was isolated in a 1.5 Eppendorf tube until the next mating opportunity 24 hours later. After all the mating opportunities were completed females remained in their isolation vials, which were checked once daily until female death. Day n* stands for the last day of the mating opportunities phase, and day x* stands for the day of female's death.

Figure A3.
Outline of the experiment to measure male harm and female resistance after lifelong cohabitation with individuals from the other sex. Male harm (A) was estimated in assays in which we measured longevity and LRS of tester females (in blue) housed with focal males (i.e., males from the selection lines; in black). Female's ability to resist the harm (B) was estimated in assays where longevity and LRS were assessed in focal females when they were housed with tester males (in blue). Each female (tester or focal) was housed with three males in a small plastic container (30 ml) with beans for the first 24 hours (day 1). On the second day, the individuals were transferred to a second oviposition container with new beans, where they remained for a week. Then, the individuals were transferred to a third oviposition container with new beans where they remained until female death. Adult offspring emerged up to day 29th after the day the experimental females were removed from the containers were counted, thereby providing a measure of LRS for each female (see text).

Figure A4.
Relationship between female body size and intrinsic female longevity (measured for virgin individuals) across females from selection regimes differing in selection associated with mating system (A), or differing in selection associated with metapopulation structure (B).
A B Figure A5. Effects of mating system (monogamy vs. polygamy) and metapopulation structure (no vs. yes) in the evolution of focal females' LRS after single mating. The interaction is not significant and there is an effect of population spatial structure on LRS, with females from structured populations producing more adult offspring along their lifetime than females from nonstructured populations (Table 2). Bars depict standard errors.

Figure A6.
Cost of reproduction: negative relationship between LRS and longevity across females from selection regimes differing in selection associated with mating system (A), or differing in selection associated with metapopulation structure (B), in assays in which females mated singly. Figure A7. Interaction of mating system (monogamy vs. polygamy) and metapopulation structure (no vs. yes) in the evolution of female resistance at generation 12 (left) and 30 (right) in assays with medium levels of sexual conflict (see text).

Figure A8.
Relationship between female body size and female longevity across females from selection regimes differing in selection associated with mating system (A) or differing in selection associated with metapopulation structure (B) in assays with medium levels of sexual conflict.

Figure A9.
Measurement of elytron length. The image shows the placement of the specimen on a fixed pedestal, and the dorsal view of an individual. Elytron length (in millimeters) was always measured on the right elytron. The measurement was always carried out parallel to the edge and connecting the most distant points of the right wing.