Predicting the effects of climate change on the fertility of aquatic animals using a meta-analytic approach
Data files
Dec 12, 2024 version files 2.63 MB
-
Phylo_Data.csv
13.71 KB
-
README.md
22.92 KB
-
TEMP_with_EffectSizes.csv
1.35 MB
-
TEMP.csv
1.25 MB
Abstract
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.
README: Predicting the effects of climate change on the fertility of aquatic animals using a meta-analytic approach
Welcome to the project repository for "Predicting the Effects of Climate Change on the Fertility of Marine and Freshwater Animals Using a Meta-Analytic Approach." This README provides an overview of the project, explains how to run the code, and outlines the data used in the analysis.
Table of Contents
- Project Overview
- Files in the Repository
- Prerequisites
- Installation Instructions
- Running the Code
- Understanding the Outputs
- Troubleshooting
- References
Project Overview
This project explores the impact of climate change on the fertility of marine and freshwater animals using a meta-analytic approach. The analysis considers data from multiple studies, conducting meta-regressions with various moderators to identify factors influencing fertility changes.
Files in the Repository
Code for Predicting the effects of climate change on fertility.R
: This R script contains code to:
- Load the data.
- Calculate biserial correlations.
- Create a phylogenetic tree.
- Perform meta-analysis and meta-regressions and interaction models.
- Perform sensitivity analysis
- Generate orchard plots for each meta-regression.
- Generate summary graphs and plots
Functions for Predicting the effects of climate change on fertility.R
: **This R script contains code to:
- create R2 and I2 functions
- functions to edit the appearance of the meta-regression orchard plots and the interaction orchard plots
TEMP.csv
: This CSV file contains the data extracted from the papers including information from various studies such as species, temperature, location, and fertility metrics. without the biserial correlations used to calculate the biserial correlations and variance in R. Cells with "NA" mean the data was not available and empty cells mean that column is not applicable to that effect size.
TEMP_with_EffectSizes.csv
: *This CSV file contains the data used for the meta-analysis, including the biserial correlations and variance, along with a column of combined effect sizes and combined variances (ES and ES_VAR), which contains a combination of biserial correlation coefficients (where means and standard deviations were available, *k=1583) and Pearson’s correlation coefficients (k=311). Cells with "NA" mean the data was not available and empty cells mean that column is not applicable to that effect size.
Final_tree
: This is the full phylogenetic tree for the large systematic map (Dougherty et al. 2024), which we trimmed to create Figure 1.
Phylo Data.csv
: This CSV file contains the retained phylogenetic tree tips and their corresponding phylum and whether the species are marine, freshwater or other aquatic animals used to create the visual phylogenetic tree.
TEMP and TEMP_with_EffectSizes CSV File Headings and Descriptions
Heading | Description |
---|---|
Paper code | A unique identifier for each paper or study. |
Paper | The name of the paper or publication where the study is documented. |
Title | The title of the paper or study. |
Publication year | The year the paper or study was published. |
DOI | The Digital Object Identifier for the paper, providing a unique link to the publication. |
GENUS SPECIES1 | The genus and species of the organism studied. |
GENUS SPECIES | The genus and species of the organism studied repeated for the phylogenetic matrix. |
Tmax | Maximum temperature reached during the study for this effect size. Units: degrees Celsius. |
CTMaxMaxTemp | Critical thermal maximum, indicating the highest temperature tolerated by the organism minus the maximum temperature. Units: degrees Celsius. |
Phylum | The taxonomic phylum of the organism(s) studied. |
Sex exposed | The sex of the organism(s) exposed to the experimental conditions. |
Fertilisation mode | The mode of fertilization (e.g., internal or external). |
Stressor variation | Description of the type of stressor variation in the study (e.g., experimental, natural). |
Climate Zone | The preferred climate zone of the species exposed (e.g., tropical, temperate, polar, etc.). |
Life-stage of animal exposed | The development stage of the organism(s) exposed in the study (e.g., gamete, juvenile, adult). |
Constant or fluctuating | Indicates whether the temperature conditions were constant or fluctuating. |
Marine or Freshwater | Indicates whether the species was marine or freshwater or other (brackish, hypersaline). |
Minimum Temperature | The minimum temperature during the study for this effect size. Units: degrees Celsius. |
Maximum Temperature | The maximum temperature during the study for this effect size. Units: degrees Celsius. |
Magnitude of Change | The degree of temperature change during the study: maximum temperature - minimum temperature. Units: degrees Celsius. |
Duration of Temperature Change Days | The length of time (in days) over which the temperature change occurred. Units: Days |
Sex of Trait | The sex associated with the trait being studied (e.g., male, female, or both). |
Fertility Trait | The specific fertility trait being measured (e.g., ovary width, sperm velocity, hatchling number). |
Units | The units of the specific fertility trait being measured (e.g., ovary width, sperm velocity, hatchling number). |
Trait category | The category of the fertility trait (e.g., gonad, gamete, reproductive output). |
Effect size code | The specific code used to denote unique effect size. |
Experiment Code | A unique code for identifying the experiment and/or group of animals from which the data came. |
Source | The source of the data or experimental results within each paper. |
Mean_lowtemp | The mean measurement of the fertility trait at low temperature conditions. Units: depend on the specific fertility trait and are listed in the units column where available. |
Standard Error_lowtemp | The standard error for mean measurements at low temperature. |
SD_lowtemp | Standard deviation for mean low temperature conditions. |
ADJSDlow | Standard deviation for mean low temperature conditions adjusted using the adjusted lown to calculate where converting from standard error. |
Confidence Intervals_lowtemp 95% | The 95% confidence intervals for measurements at low temperature. |
Proportion_lowtemp | Proportion data for the low-temperature fertility trait where means and standard deviations are not available and only percentage/proportion data is available. |
n_lowtemp | The number of samples at low temperature. |
adjustedlown | The number of samples at low temperature adjusted by dividing the n_lowtemp by count_shared_control. |
Mean_hightemp | The mean measurement of the fertility trait at high temperature conditions. Units: depend on the specific fertility trait and are listed in the units column where available. |
Standard Error_hightemp | The standard error for measurements at high temperature. |
SD_hightemp | Standard deviation for high temperature conditions. |
ADJSDhigh | Standard deviation for low temperature conditions adjusted using the adjusted lown to calculate where converting from standard error. |
Confidence Intervals_hightemp 95% | The 95% confidence intervals for measurements at high temperature. |
Proportion_hightemp | Proportion data for the high-temperature fertility trait where means and standard deviations are not available and only percentage/proportion data is available. |
n_hightemp | The number of samples at high temperature. |
adjustedhighn | The number of samples at high temperature adjusted by dividing the n_hightemp by count_shared_control. |
Test | The type of statistical test used in the study. |
Final.Effect.Size | The type of statistical test used in the final analysis (Biserial or Pearsons). |
Test_value | The value derived from the statistical test. |
Significance | Indicates whether the result is statistically significant. |
Effect size direction | Indicates whether the effect size is positive or negative. |
Number.of.Animals | Total number of animals in the study, distinct from the total of the n_lowtemp or n_hightemp if discrepancies arise due to group or replicate comparisons. |
Total.Number.of.Replications | Total number of replications for the effect size. |
N_total | Total sample size at the individual level. |
N_totaladj | Total sample size used in the effect size calculations at the individual level (adjustedhighn + adjustedlown or Number.of.Animals where necessary). |
Pooled.Standard.Deviation | The standard deviation derived from pooling across all samples, using the adjusted sample sizes. |
Bias Correction Factor | A correction factor used to account for unequal sample sizes when calculating Hedges g effect sizes. Formula: 1 - (3 / (4*(adjustedhighn + adjustedlown - 2) - 9)). |
PreD | Hedges g effect size before negative trait correction. Formula: ((Mean_hightemp - Mean_lowtemp) / Pooled Standard Deviation) * Bias Correction Factor. |
D effect size | The standardized effect size using Hedges g times the negative trait correction. If a decrease in the trait has a positive effect on fertility, the effect size direction is reversed. |
D_VAR | The variance of Hedges g. Formula: ((adjustedhighn + adjustedlown) / (adjustedhighn * adjustedlown)) + (D effect size^2) / (2 * (adjustedhighn + adjustedlown - 2)). |
r | The correlation coefficient extracted from papers. |
r_VAR | The variance of the correlation coefficient. Formula: (1 − r^2)^2 / (N_totaladj − 1). |
Biserial_Correlation | Biserial correlations calculated in R using escalc. |
Biserial_Correlation_Variance | Biserial correlation variance calculated in R using escalc. |
ESPre | Combined column of Biserial correlations calculated in R and Pearsons correlations for effect sizes without means and standard deviations, before negative trait correction. |
Negative trait Correction | Indicates whether a correction was applied based on whether an increase in the fertility trait measured would have a positive or negative effect on fertility. |
ES | Biserial correlations calculated in R without means and standard deviations, multiplied by the negative trait correction. |
ES_VAR | The variance associated with the biserial and Pearsons correlations for the effect sizes. |
shared.control.id | A unique identifier of effect sizes that share a control value. |
count_shared_control | The number of times the control value in that row has been used. |
Phylo Data CSV File Headings and Descriptions
Heading | Description |
---|---|
TipLabel | The unique species as they appear on the phylogenetic tree tips. |
Phylum | Phylum of the individual species. |
Marine.or.Freshwater | The habitat type of each species and whether it is marine, freshwater, or other. |
## Prerequisites
To run the code in this repository, you need:
* R (version 4.0 or later)
* RStudio (optional but recommended)
* A complete set of R packages:
```r
# Install required packages
install.packages(c("ape" = "5.7-1", "bookdown" = "0.39", "car" = "3.1-2", "dplyr" = "1.1.4",
"forcats" = "1.0.0", "ggplot2" = "3.4.0", "ggpubr" = "0.6.0",
"ggtree" = "2.2.1", "metafor" = "3.0-2", "tidyverse" = "1.3.1"))
# Load pacman to simplify loading other packages
pacman::p_load(devtools, tidyverse, metafor, patchwork, R.rsp, emmeans)
# Install orchaRd from GitHub
devtools::install_github("daniel1noble/orchaRd", force = TRUE)
# Load libraries
library(readr)
library(ape)
library(orchaRd)
library(metafor)
library(car)
library(patchwork)
library(rotl)
library(ggplot2)
library(MCMCglmm)
library(Hmisc)
library(dplyr)
library(tidyr)
library(lme4)
library(MuMIn)
library(treebase)
library(viridis)
library(bookdown)
library(metadat)
library(stringr)
library(Matrix)
library(tidyverse)
library(forcats)
library(purrr)
library(tibble)
library(magrittr)
library(knitr)
library(writexl)
library(RColorBrewer)
Installation Instructions
Ensure R and RStudio are installed. Use the code block above to install all required packages and additional packages using pacman
. This includes installing orchaRd
from GitHub.
Running the Code
- Open the
Code for Predicting the effects of climate change on fertility.R
and for Functions forPredicting the effects of climate change on fertility.R
script in RStudio or any other R environment. - Ensure the path to the
TEMP.csv
andTEMP_with_EffectSizes.csv
file is correct. - Follow the instructions within the script to:
- Load the data.
- Calculate the effect sizes
- Create the phylogenetic tree.
- Perform meta-analysis and meta-regressions.
- Run functions
- Generate orchard plots.
- Execute the script to obtain the results.
Understanding the Outputs
- Phylogenetic Tree: A visualization of relationships among species in the dataset.
- Meta-Analysis Results: Summary statistics for the overall effect size and heterogeneity.
- Meta-Regressions: An analysis showing the impact of different moderators on effect size.
- Orchard Plots: Graphical representations of meta-regressions, illustrating relationships between moderators and fertility outcomes.
Troubleshooting
- If you encounter errors during package installation, ensure you have an active internet connection and appropriate permissions.
- Verify the path to
Predicting the effects of climate change on the fertility data.csv
is correct. - Ensure the R script and data file are in the same directory, or update the script with the correct path.
- If you experience errors related to package loading, double-check package installation and attempt to reinstall any packages causing issues.
References
For more information on methodologies and the R ecosystem, refer to:
- "Meta-Analysis in Ecology and Evolution" by Julia Koricheva, Jessica Gurevitch, and Kerrie Mengersen.
- R Documentation at r-project.org.
For further questions, issues, or suggestions, please contact the project maintainer or submit an issue on the repository's issue tracker.
Thank you for using this repository. I hope this README helps you successfully run the analysis and understand the results. Good luck with your work!
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.