Predicting the effects of climate change on the fertility of aquatic animals using a meta-analytic approach

Given that reproductive physiology is highly sensitive to thermal stress, there is increasing concern about the effects of climate change on animal fertility. Even a slight reduction in fertility can have consequences for population growth and survival, so it is critical to better understand and predict the potential effects of climate change on reproductive traits. We synthesised 1894 effect sizes across 276 studies on 241 species to examine thermal effects on fertility in aquatic animals. Our meta-analysis revealed that external fertilisers tend to be more vulnerable to warming than internal fertilisers, especially in freshwater species. We also found that increased temperature is particularly detrimental for gametes, and that under certain conditions, female fertility is more sensitive to warming than male fertility, challenging the prevailing view that males are more vulnerable. This work provides valuable new insights into the effects of temperature on fertility, with potential consequences for population viability. 

Methods

Literature Search and Data Collection 

This study was not pre-registered. To identify relevant studies for our meta-analysis, we used a large systematic map (Dougherty et al. 2024), which was created by searching published literature using the ISI Web of Science Core and following the ROSES Reporting Standards for Systematic Evidence Syntheses guidelines and the PRISMA guidelines (Haddaway et al. 2018, O’Dea et al. 2021,Supplementary Table 8). Studies included in this systematic map were peer-reviewed journal articles or book chapters that measured one or more reproductive traits at two or more temperatures in non-human animals. The systematic map includes 427 papers that contain data on the effects of temperature on fertility in aquatic animals. We excluded 151 of these from our meta-analysis for a variety of reasons (Supplementary Figure 1, Supplementary Table 9). For example, we excluded studies with data on gonad developmental stages and gonad gene expression, which have no clear link to fertility. We included gonad traits such as gonadosomatic index and gonad size, which have a clearer relationship with reproductive success. We also excluded studies examining seasonal temperature variations outside the breeding season, which was not relevant to our question. For studies with insufficient information for calculating effect sizes, we contacted the authors to request the missing data. Out of 37 authors, four responded, with one providing the necessary data. We thus had to exclude the remaining 36 studies (Supplementary Figure 1).  

Our final analysis included 276 papers, and each paper was assigned to one of the authors (AC, IG, EM, or CH). The data were then compiled, cross-checked, and validated by the first (AC) and last author (NP). A total of 1894 effect sizes were extracted and used in the final analysis. We extracted data on the effects of temperature on gamete traits (e.g., gamete size, quantity, or quality), gonad traits (e.g., gonad size, gonadosomatic index), and reproductive output (Supplementary Table 2). For each effect size, we recorded information on taxonomic family, class and phylum, fertilisation mode, which sex was exposed to the temperature change, whether temperature variation was natural or experimentally manipulated, the duration of exposure to the temperature change, and the developmental stage during which individuals were exposed to the temperature change. 

Effect size calculation 

Data were extracted from tables, text, or figures using WebPlotDigitizer (Rohatgi 2022). For studies using two or three discrete temperature treatments (e.g., low versus high temperatures), we extracted the mean and standard deviation and used these to calculate a biserial correlation coefficient and its sampling variance using the ‘escalc’ function in the ‘metafor’ package vers.4.6-0 (Viechtbauer 2010). The ‘escalc’ function uses equation 12 from Jacobs and Viechtbauer (2017) to calculate the sampling variance. To account for shared-control non-independence, we first identified all cases where a control group was compared to several treatment groups. We then reduced the weight of the control group in the analyses by dividing the sample size of such a group by as many times that group has been compared to a treatment group before calculating effect sizes.  

For studies where fertility was measured across a range of temperatures (i.e. where it was appropriate to treat temperature as a continuous variable), we calculated a Pearson’s correlation coefficient. For any studies that only provided test statistics, such as a chi-square test, we converted these to correlation coefficients, using the formulas provided in Borenstein et al. (2011). For both biserial correlation coefficients and Pearson’s correlation coefficients, a positive value indicated higher fertility at higher temperatures, whereas a negative value indicated a reduction in fertility at higher temperatures.  

In total, we obtained 1326 effect sizes for marine species, 473 effect sizes for freshwater species, and 95 effect sizes for species that were neither freshwater nor marine (e.g., hypersaline or estuarine). This included 241 species from nine phyla with Chordata (38.2%) having the largest representation, followed by Arthropoda (20.7%) and Mollusca (17.2%) (Supplementary Figure 2). All species in our dataset were ectotherms. 

Data analysis 

All analyses were completed in R version 4.3.1 (R Core Team 2023), and figures were generated using the ‘orchaRd’ package vers. 2 (Nakagawa et al. 2023) and the ‘ggplot2’ package vers. 3.5.1 (Wickham 2016). These analyses and all figures use a combination of biserial correlation coefficients (where means and standard deviations were available, k=1583) and Pearson’s correlation coefficients (k=311). To construct a phylogenetic tree for the species included in our dataset, we used the tree for the 1,191 species in the full systematic map (Dougherty et al. 2024) and then used the ‘set.diff’ and ‘drop.tip’ functions in the ‘ape’ package vers. 5.7-1 (Paradis et al. 2004) to remove non-aquatic species (Figure 1). The variance-covariance matrix representing the phylogenetic relatedness between all species in our dataset was used as a random effect in all of our meta-analytic models described below (Paradis et al. 2004). We also included the following random effects in all models: ‘paper code’ to account for multiple effect sizes from the same study, ‘animal group ID’ to control for effect sizes from the same paper that used different groups of animals, ‘species ID’ to account for the same species being used across different papers, ‘shared control ID’ to account for effect sizes from the same paper which shared a control treatment, and ‘effect size code’ as a unit level effect to measure residual heterogeneity (Viechtbauer 2010).   

We first assessed the overall effect of temperature on fertility with a random effects model using the ‘rma.mv’ function in ‘metafor’. We estimated heterogeneity in our main effects model using I2 as an estimate of the proportion of variance explained due to differences between the levels of a random effect (Nakagawa and Santos 2012). We also display prediction intervals on figures as a measure of heterogeneity (Cooper et al. 2019). We then tested a range of moderators that may influence the effect of temperature on fertility, using multi-level meta-regression models with the same random effects described above. For each moderator, we calculated a marginal R2 to describe the percentage of heterogeneity explained by the moderator model compared to the main effects model. Statistical significance of moderators was determined using omnibus tests of meta-regressions with intercepts. We obtained mean effect size estimates (β) and 95% confidence intervals for each level of the moderator by removing the intercept from the meta-regression model. For statistically significant moderators with more than two levels, we then examined which factor levels differed from one another. We did this using the L argument in the ‘anova’ command to perform linear contrasts in meta-regressions without an intercept.  

We tested ten categorical moderators and four continuous moderators (Supplementary Table 1). We first looked at three categorical moderators to examine whether there were any differences in the effects of temperature on fertility that could be explained by the type of fertility trait being examined: (i) “sex of the trait” (male, female, both) and (ii) “type of fertility trait” (gamete, gonad or reproductive output). The category “both” for “sex of trait” refers to combined fertility traits, influenced by both male and female fertility changes, such as fertilisation success. Sex of trait differs from the sex exposed to temperature as some papers only exposed one sex but measured traits that could be affected by both sexes. For fertility traits that could not be assigned to a specific sex, such as fertilisation success or offspring number (k=479), we ran a model testing whether the effect on these fertility traits varied depending on whether only males, only females, or both sexes had been exposed to increased temperature (iii) “Sex exposed (for combined traits)”. This allowed us to test for additive or multiplicative effects of exposing both sexes to a higher temperature compared to just one sex.  

Five categorical moderators and one continuous variable were used to test whether variation in the effect of temperature on fertility could be explained by the organism’s physiological or life history traits: (i) “developmental stage of animal exposed” to compare effects due to the exposure of gametes, embryos and juveniles, or adults, (ii) “habitat type” to test for differences between freshwater and marine species, (iii) “fertilisation mode” to test for differences between internal and external fertilisers, (iv) “phylum” to test for different effects across taxonomic groups (v) “climate zone” to test for differences among species from tropical, subtropical, temperature, or polar regions, and (vi) each species’ thermal tolerance (CTmax). To examine the relevance of thermal tolerance to effects on fertility, we subtracted the elevated temperature organisms were exposed to from that species’ critical thermal limit (CTmax), using data for 13 species from four phyla in our dataset from the GlobTherm database (Bennett et al. 2018, Supplementary Figure 3).  As with the main dataset, the majority of the effect sizes (k=184 out of 250) came from Chordata species (n=6), and four Mollusc species contributed 47 effect sizes. The other species were Echinoderms (n=2, k=18) and Arthropods (n=1, k=1). In addition to examining the effects of each of these moderators individually, we looked at the interaction between fertilisation mode and habitat type to test whether the effects of temperature on internal versus external fertilisers differed between freshwater and marine environments. 

We then examined the effects of a study’s experimental design, and the type of temperature exposure experienced by the organism, focusing on the duration of the temperature change, the magnitude of the temperature change, whether organisms experienced constant or fluctuating temperatures, and whether the temperature variation was experimentally manipulated or experienced under natural conditions. We also looked at whether the interaction between the magnitude of temperature change and the duration of the temperature change influenced the effect of temperature on fertility. 

Lastly, we examined sex-specific effects using a series of multi-level meta-analytic models with interaction effects. More specifically, we tested the interaction between focal sex (male or female) and (i) “type of fertility trait” (gamete, gonad or reproductive success), (ii) “fertilisation mode” (internal or external fertilisation), and (iii) “developmental stage of animal exposed” (gamete, early development, or adulthood). Alluvial plots show the distribution of effect sizes across the moderator levels for these interactions (Supplementary Figures 4-6). For statistically significant interactions, we conducted a leave-one-out analysis to assess the stability of the interactions and determine if the results changed significantly when excluding one group. This was done using the ‘metafor’ and ‘dplyr’ packages vers. 1.1.4 for data manipulation (Wickham et al. 2023). 

Publication bias tests 

Given the well-documented pattern that studies with statistically significant results or larger effects tend to be published earlier (Jennions and Møller 2002), we used a meta-regression with mean-centred publication year as a continuous moderator to test for time-lag bias using the same random effects described above. We also used the square root of the inverse of the sample size for each effect size as a moderator to test for small-study effects (Nakagawa et al. 2022).

Sensitivity analyses 

We chose to use effect sizes based on correlation coefficients rather than standardised mean differences for our analyses, because we were interested in examining the strength and direction of the relationship between temperature and fertility, rather than testing for differences between two groups. However, given that most of our effect sizes were based on means and standard deviations, we performed a sensitivity analysis to test whether our results are robust regardless of the type of effect size chosen. To this end, we re-ran all of our models using standardised mean differences (Hedges’ g effect sizes). Hedges’ g is more robust to unequal sampling and small sample sizes than Cohen’s d (Rosenberg et al. 2013). Correlation coefficients or test statistics were converted to Hedges’ g effect sizes using the equations provided by Borenstein et al. (2011).  

To ensure our results were robust given non-independence of effect sizes extracted from the same studies, we conducted cluster robust estimate of the variance models (RVE) for our moderator models using the ‘robust’ function in ‘metafor’ and compared them to the standard random effects models (Hedges et al. 2010).The RVE models provide a more conservative estimate of significance with larger standard errors and wider confidence intervals and can prevent Type 1 errors (Tanner-Smith and Tipton 2014). Given that tolerance of exposure to the thermal stressor (Bigelow 1921, Cerdá et al. 1998, Tang et al. 2000, 2000, Cerdá and Retana 2000, Armstrong et al. 2009, Rezende et al. 2014), we also conducted a sensitivity analysis using the duration of the temperature change as an additional fixed effect in our moderator models. These models only included 1,362 out of the 1,894 effect sizes in our full dataset, because we lacked information on the duration of temperature exposure for the remaining effect sizes – 28% of papers in our dataset did not report this information. 

Lastly, we ran two additional multi-level meta-regression models to examine whether there were any systematic differences depending on how our effect sizes were calculated. The moderator for the first model was the final effect size type (biserial correlation coefficient or Pearson’s correlation coefficient), and the moderator for the second one was the type of statistical test that the effect sizes were originally converted from (e.g. chi-square test, Spearman’s correlation). Both of these models included the same random effects as described above.