Predation risk shapes the degree of placentation in natural populations of live‐bearing fish

Abstract The placenta is a complex life‐history trait that is ubiquitous across the tree of life. Theory proposes that the placenta evolves in response to high performance‐demanding conditions by shifting maternal investment from pre‐ to post‐fertilisation, thereby reducing a female’s reproductive burden during pregnancy. We test this hypothesis by studying populations of the fish species Poeciliopsis retropinna in Costa Rica. We found substantial variation in the degree of placentation among natural populations associated with predation risk: females from high predation populations had significantly higher degrees of placentation compared to low predation females, while number, size and quality of offspring at birth remained unaffected. Moreover, a higher degree of placentation correlated with a lower reproductive burden and hence likely an improved swimming performance during pregnancy. Our study advances an adaptive explanation for why the placenta evolves by arguing that an increased degree of placentation offers a selective advantage in high predation environments.

The model fit was assessed using a posterior predictive check on the predictions of embryo mass. log(EM i,j ) ∼ N (α pop×year + α mother + β 1,pop×year stage i,j + β 2,pop×year stage 2 i,j + β 3 mBF i + β 4 mSL i + β 5 where log(EM i,j ) corresponds to the ln-transformed mean embryo dry mass of the jth brood carried by the ith female. α pop×year corresponds to the year-specific intercept of a given population and α mother to the mother identity. β 1,pop×year stage j,i is the year-specific slope of a given population on the developmental stage of the jth brood carried by the ith female (second order polynomial), mBF i the arcsin square-root transformed proportion of body fat of the ith female, mSL i the standard length of the ith female, and σ the residual standard deviation.

Predation-specific Matrotrophy Indices
Predation-specific Matrotrophy Indices (MI's) were estimated using the Bayesian programming environment JAGS [8]. For this, ln-transformed embryo dry mass was fitted in a linear model as a function of the developmental stage of embryos (second order polynomial), high and low predation risk (i.e., piscivorous predator species present or absent), the proportion of maternal body fat, maternal standard length, the interaction between the developmental stage of embryos and the proportion of maternal body fat, and the interaction between the developmental stage of embryos and predation risk. For each population, the model estimates an additional year-specific intercept to correct for differences between populations in a specific year. Moreover, the model includes the mother identity as additional intercept to account for variation among females that is not accounted by maternal body fat and standard length. We used N(0,5 2 ) priors for all fixed effects. For mother identity, we used a N(0,σ 2 ) prior with the common variance σ 2 ∼ Inv-Gamma(0.01, 0.01), parameterized with shape and rate parameters. The year-specific intercepts and year-specific slopes on the developmental stage of the embryos (second order polynomial) were estimated using N(µ, σ 2 ) priors, with the common mean µ ∼ N(0,1 2 ) and the common standard deviation σ ∼ Student-t (0,∞) (0, 1 2 , 1), parameterized with mean, variance, and degrees of freedom. The MI for high and low predation females was subsequently calculated by dividing offspring dry mass at birth (developmental stage 45) by egg dry mass at fertilization (developmental stage 2) that were predicted for a given predation regime (i.e., piscivorous predator species present or absent) and for a female of overall average standard length (i.e., y standard length =53 mm) and proportion of body fat (i.e., y body fat =0.16). Three MCMC chains were run for 1,500,000 iterations, with a burnin of 500,000 and a thinning of 1000. Convergence was assessed by visual examination of the traces and by checking thatR < 1.01. The model fit was assessed using a posterior predictive check on the predictions of embryo mass. log(EM i,j ) ∼ N (α pop×year + α mother + intercept+β 1 stage i,j + β 2 stage 2 i,j + β 3 predation i + β 4 mBF i + β 5 mSL i + β 6 where log(EM i,j ) corresponds to the ln-transformed mean embryo dry mass of the jth brood carried by the ith female. α pop×year corresponds to the year-specific intercept of a given population and α mother to the mother identity. intercept is the overall intercept, β 1 stage j,i the developmental stage of the jth brood carried by the ith female (second order polynomial), predation i the predation risk experienced by the ith female, mBF i the arcsin square-root transformed proportion of body fat of the ith female, mSL i the standard length of the ith female, and σ the residual standard deviation.

Path analysis
Differences in reproductive allotment among populations could be due to effects on several lifehistory traits. For instance, reproductive allotment could be decreased by reducing brood size or superfetation, which in turn decreases the number of embryos (i.e., fecundity). Alternatively, reproductive allotment diminishes when producing smaller eggs at fertilization or offspring at birth. We used a path analysis to determine the contribution of each of these life-history traits to differences in reproductive allotment among predation regimes. In total, three (generalized) linear mixed effect models implemented in MCMCglmm [3] were used to estimate all paths.
In the first model, we estimated the direct effect of predation risk on z-standardized egg dry mass at fertilization (developmental stage 2) (ze), z-standardized offspring dry mass at birth (developmental stage 45) (zo), z-standardized average brood size for a given mother (zb), and the degree of superfetation (s) in a multivariate model: [ze, zo, zb, s] ∼ µ + Xβ + Zr, where µ is the intercept vector, β the vector of fixed effects, X the corresponding design matrix, and Z is a design matrix for additional random terms r. z-standardized egg dry mass at fertilization, offspring dry mass at birth, and the average brood size for a given mother are assumed to be Gaussiandistributed. Superfetation was formulated in a generalized linear mixed model framework using a log link for the Poisson-distributed response. Fixed effects include predation risk (i.e., piscivorous predator species present or absent), z-standardized proportion of maternal body fat, z-standardized maternal standard length, and the z-standardized developmental stage of the most-developed brood in the case of superfetation. Besides an error, the random terms include population and river identity accounting for spatio-temporal non-independence of observations. Moreover, the interaction between population identity and year quantifies variation of a given population between years. The multivariate framework allows for the covariance between the residuals of all responses. As priors, we used the default normal priors for the fixed effects with the expected value of 0 and variance 10 12 , and inverse-Wishart priors for the variances with the expected value of 1 and degree of belief of 3.002. The number of iterations was 5,500,000, with a burnin of 500,000 and a thinning of 5000.
In the second model, maternal fecundity was fitted as a function of superfetation and average brood size for a given mother, as changes in both brood size and superfetation will affect the number of embryos: where f i corresponds to the fecundity (i.e., number embryos) of the ith female. α corresponds to the overall intercept, and α pop , α river , and α pop×year to the random intercepts. zs i is the z-standardized degree of superfetation, zb i the z-standardized average brood size, zmBF i the z-standardized proportion of maternal body fat (second order polynomial), zmSL i the z-standardized maternal standard length, and zlatestStage i the z-standardized developmental stage of the most developed brood of the ith female. As priors, we used the default normal priors for the fixed effects with the expected value of 0 and variance 10 12 , and inverse-Wishart priors for the variances with the expected value of 1 and degree of belief of 0.002. The number of iterations was 1,500,000, with a burnin of 500,000 and a thinning of 1000.
The third model subsequently predicts the absolute dry reproductive allotment as a function of fecundity, egg dry mass at fertilization, and offspring dry mass at birth: where zra i corresponds to the absolute dry reproductive allotment of the ith female. α corresponds to the overall intercept, and α pop , α river , and α pop×year to the random intercepts. zf i is the z-standardized maternal fecundity (i.e., number embryos), ze i the z-standardized egg dry mass at fertilization (developmental stage 2), zo i the z-standardized offspring dry mass at birth (developmental stage 45), zmBF i the z-standardized proportion of maternal body fat, zmSL i the z-standardized maternal standard length of the ith female, and σ the residual standard deviation. As priors, we used the default normal priors for the fixed effects with the expected value of 0 and variance 10 12 , and inverse-Wishart priors for the variances with the expected value of 1 and degree of belief of 0.002. The number of iterations was 1,500,000, with a burnin of 500,000 and a thinning of 1000. For all models, convergence was assessed by visual examination of the traces and by checking that the autocorrelations of the parameter chain was less than 0.1. In addition, each model was re-fitted as a function of an intercept only (null model) to compare the deviance information criterion (DIC) of the full model against the DIC of the null model (∆DIC).

Measurements of water parameters
The water velocity was measured at each location to the nearest 0.01 m· s −1 with a Höntzsch Vane Wheel FA current meter (type ZS30 GFE md20 T/100-2/p10, Höntzsch Instruments, Waiblingen, Germany). Depending on the uniformity of the flow, the water velocity was taken at 9-17 incremental observation points across a transect of the stream. At each observation point, the mean water velocity was defined as the average of three repeated measurements at a height above the stream bed equal to 0.4 times the depth at that location. When the water depth exceeded 60 cm, the mean water velocity was calculated as the average between the velocities measured at 0.2 times the water depth and 0.8 times the water depth [5]. Each location was additionally characterized by measuring salinity (S) by using the ExStik II pH/conductivity/TDS meter (Extech Instruments, Nashua, USA), hardness (mg· L −1 ) with a titrimetric color-test kit (Merck KGaA, Darmstadt, Germany), ammonium concentration (NH + 4 ) (mg· L −1 ) as a proxy for the nitrogen loading of a stream by using a colorimetric ammonium test (Merck KGaA, Darmstadt, Germany), phosphate concentration (PO 3− 4 ) (mg· L −1 ), as phosphorus is an important determinant of primary production in freshwater ecosystems [12] by using a colorimetric phosphate test (Merck KGaA, Darmstadt, Germany), and dissolved oxygen (%) with the ExStik DO600 meter (Extech Instruments, Nashua, USA). Salinity and dissolved oxygen were measured 1-3 times at each location, whereas hardness, ammonium concentration, and phosphate concentration were measured only once.
Salinity was then predicted for each location at 25 • C in the following linear model by Maximum Likelihood: where S i,j corresponds to the jth measurement of salinity at location i, α to the overall intercept, location i to the ith location, temperature i,j to the jth measurement of water temperature at location i, and σ to the residual standard deviation.
Dissolved oxygen was predicted for each location at 13:00 hours and 25 • C in the following linear model by Maximum Likelihood: where oxygen i,j corresponds to the jth measurement of dissolved oxygen at location i, α to the overall intercept, location i to the ith location, temperature i,j to the jth measurement of water temperature at location i, time i,j to the z-standardized numeric day time of the jth measurement of dissolved oxygen at location i, and σ to the residual standard deviation.

Relationship between population-specific Matrotrophy Indices and environmental variables
Either by influencing egg mass at fertilization or offspring mass at birth, many population-specific factors are likely to contribute to the observed interpopulation variation in the degree of placentation. To quantify the potential effects of additional water parameters at each location on the degree of placentation, the population-specific Matrotrophy Indices (MI's) estimated in the model described in Equation 1 were fitted as a function of all measured environmental variables (predation risk, salinity, water velocity, hardness, NH + 4 , PO 3− 4 , and dissolved oxygen) in a linear mixed effect model using Restricted Maximum Likelihood. The continuous predictors were z-standardized in order to make them comparable [11]. Moreover, the model includes river identity as additional intercept to account for variation among rivers that is not accounted by the measured environmental parameters.
where zMI i corresponds to the Matrotrophy Index for the ith population in a specific year estimated in the model described in Equation 1. α river corresponds to the river-specific intercept and intercept to the the overall intercept. predation i is the predation risk (i.e., piscivorous predator species present or absent), zsalinity i the z-standardized water salinity (S), zhardness i the z-standardized hardness (mg·L −1 ), zwater velocity i the z-standardized mean water velocity (m·s −1 ), zNH + 4,i the z-standardized ammonium concentration (mg·L −1 ), zPO 3− 4,i the z-standardized phosphate concentration (mg·L −1 ), zO 2,i the z-standardized dissolved oxygen (%), and σ the residual standard deviation.

Embryo growth during gestation
Embryo growth during gestation was estimated as the exponential relationship between embryo dry mass and the developmental stage of embryos (second order polynomial) in a linear mixed effect model by Restricted Maximum Likelihood [4]. The model included maternal dry mass, and the interaction between the developmental stage of embryos and maternal dry mass as additional fixed effects (Fig.  S7).
Mother identity was fitted as random intercept to correct for pseudo-replication. Population, year, and river identity were fitted as random intercepts accounting for spatio-temporal non-independence of observations. Moreover, the interaction between population identity and year was fitted as random intercept to quantify the variation of a given population between years: log(EM i,j ) ∼ N (α + α mother + α pop + α year + α river + α pop×year + where log(EM i,j ) corresponds to the ln-transformed mean embryo dry mass of the jth brood carried by the ith female. α corresponds to the overall intercept, α mother , α pop , α year , α river , and α pop×year to the random intercepts. stage i,j is the developmental stage of the jth brood carried by the ith female, mM i the dry mass of the ith female, and σ the residual standard deviation.

Relationship between reproductive allotment, locomotor performance, and survival probability
Fleuren et al. 2019 [2] studied the locomotor performance in three placental live-bearing fish species (family Poeciliidae) that exhibit different levels of superfetation. Particularly, they used computervision based techniques to study changes in body shape (e.g. volume) and three-dimensional fast-start escape performance (e.g. maximum escape velocity) during pregnancy in Poeciliopsis turneri, Heterandria formosa, and Phalloptychus januarius. The slope between body shape and escape performance was not significantly different between the three species. Here we used the relationship between female volume and maximum escape velocity in Poeciliopsis turneri to predict the maximum escape velocity in Poeciliopsis retropinna. P. retropinna and P. turneri are similar regarding their degree of post-fertilization maternal provisioning and degree of superfetation. We derived female volume from female wet mass, assuming an uniform tissue density (ρ) of 1 g·cm −3 [10]. The relationship between female volume and maximum escape velocity is then given by: whereV max corresponds to the normalized maximum escape velocity (normalized for standard length SL), α to the overall intercept, and M ρSL 3 to the normalized female volume (normalized for standard length SL 3 ) derived from the tissue density ρ, female wet mass M, and the average standard length of a pregnant P. retropinna female. CR max is the maximum caudal peduncle curvature rate in the kinematic stage 2 of the fast-start escape response.
The estimated maximum escape velocity for a given wet mass was then used to predict the probability of evading the strike of a natural predator. Based on a study by Walker et al. 2005 [15] with Guppies (Poecilia reticulata), the probability of evading a predator (ω) was predicted using four parameters: (1) the initial distance between predator and prey (D pred ), (2) the average speed of the predator (ν pred ), (3) the evasion path of the prey relative to the strike path of the predator (θ pred ), and (4) the maximum velocity of the prey (vmax=V max · SL Guppy ) estimated in Equation 10:

Simulating selection
To illustrate that the strength of selection needed to explain observed rates of evolution, assuming that genetic drift is not involved, can be extremely weak, we have simulated the frequency p 2 of a high predation genotype (A 1 A 1 ), and the complementary frequency q 2 of a low predation genotype (A 2 A 2 ) over time (i.e., generations). The HP genotype is assumed to have a survival advantage (s) of 1.2%, and there is no mixing of the genotypes. The initial frequency (p 2 0 ) of the advantageous genotype in the population is assumed to be 0.01.
The initial genotype frequencies are given by: The fitness ω of the genotypes A 1 A 1 and A 2 A 2 are given by: The mean fitness of all the individuals in the population at time t is then given by: The genotype frequencies after selection can then be calculated by: 2 Supporting Results

Relationship between maternal traits and life-history
Both the proportion of maternal body fat and standard length are significantly associated with maternal fecundity (Table S9) Table S6; Fig. S2). In addition, egg dry mass at fertilization (mg) increases as a function of increased maternal standard length (t 117.432 =4.008, P<0.001; Table S2

Relationship between maternal traits and life-history
Maternal traits can profoundly influence the offspring phenotype and maternal life-history [7]. Consistent with previous findings in P. retropinna [4], larger females have greater fecundity, produce larger eggs at fertilization, and offspring at birth (Fig. S2). Furthermore, females that have more fat reserves produce fewer but larger offspring at birth, without investing more in egg size at fertilization. The negative correlation between the proportion of maternal body fat and fecundity may reflect a trade-off between offspring size and number; the production of large offspring may necessarily entail the production of fewer offspring owing to the limited size of the female body cavity [14].    Egg dry mass at fertilization (i.e., developmental stage 2) (n=129), b offspring dry mass at birth (i.e., developmental stage 45) (n=190), c proportion of egg fat at fertilization (i.e., developmental stage 2) (n=117), d proportion of offspring fat at birth (i.e., developmental stage 45) (n=190), e brood size (n=943), f maternal fecundity (i.e., number embryos) (n=411), g degree of superfetation (n=449), and h abortion incidence (n=463) (± 95% CI) as a function of the proportion of maternal body fat (left panels) and standard length (right panels) estimated in the models described in Table S2--S10. The models are predicted for a high predation female and account for maternal standard length (left panels) and the proportion of maternal body fat (right panels), which are kept constant at the overall population mean (i.e., y body fat =0.16, y standard length =53 mm).
In f and g, the developmental stage of the most-developed brood carried by the female is kept constant at the overall median (i.e., developmental stage 42.5). Data points correspond to the river-specific raw data. P-value is given at the top. . The potential mechanisms behind the increased post-fertilization maternal provisioning by larger and better-conditioned females, however, is unclear [4]. By contrast, maternal fecundity is proposed to increase with female size, as a consequence of more space available in the female's body cavity [13]. The physical constraint of the body cavity is also displayed by the populations where P. retropinna co-occur with P. dovii as the only predator species. Independently of the proportion of maternal body fat and standard length, these females produce very large offspring at birth (Fig. 2), but carry significantly fewer broods at different developmental stages (i.e., superfetation; Fig. 2). These females were mainly collected from a single stream (Rio Sucio) with exceptional high nitrogen loading (NH + 4 concentration measured in 2017 2.71 × larger than average), and probably nutrient-rich water. As a result, the females are very large (Table S20; Fig. S6), contain a large amount of body fat (Table  S21; Fig. S6), and invest more in offspring size at birth (Fig. 2). Here, we additionally show that the degree of superfetation does not correlate with the proportion of maternal body fat or size in P. retropinna (Table S6; Fig. S2).  Table S11--S16. a Salinity, b mean water velocity, c hardness, d ammonium concentration (NH + 4 ), e phosphate concentration (PO 3− 4 ), and f dissolved oxygen (± 95% CI). Data points (red: high predation; blue: low predation) correspond to the 'jittered' raw data. Sample size and P-value are given at the top. These data tentatively suggest that there are no obviously large differences in water quality parameters between study locations with and without predators. However, these results should be interpreted with care, because salinity and dissolved oxygen were measured only 1-3 times at each location, and hardness, ammonium concentration, and phosphate concentration were measured even only once. Water quality parameters are likely to vary during the day and throughout the year. A single or only a few measurements per location taken at different times during the day, on different days, and even in different years are therefore unlikely to accurately reflect the yearly mean local abiotic conditions experienced by the fish in each population. values between -1 and 1 [11]. Please note that salinity and dissolved oxygen were measured only 1-3 times at each location, and hardness, ammonium concentration, and phosphate concentration were measured even only once. The water quality parameters are likely to vary during the day and throughout the year. A single or only a few measurements per location taken at different times during the day, on different days, and even in different years are therefore unlikely to accurately reflect the yearly mean local abiotic conditions experienced by the fish in each population. Thus, care must be taken when interpreting the relationship between the MI and the water quality parameters. Still, these data show that independently of the measured water quality parameters, we still predict a higher degree of placentation in high predation populations.