Data collection and preparation—Dutch waters
Study specimens
Between 2006 and 2019, 1457 deceased harbour porpoises in The Netherlands were collected for post-mortem investigations and diet analysis, with the necropsies conducted following internationally standardized guidelines69. For this study focussing on female life history, we selected all females > 115 cm, as smaller animals may be maternally dependent and can be considered young of the year23. Cases in DCC5, which represent carcass remains, and of which reproductive organs were not assessable or present due to scavenging or incompleteness of the carcasses were excluded. The reproductive organs of the female porpoises > 115 cm (n = 328/1457) were macroscopically inspected to differentiate between immature and mature animals, with the presence of ovarian corporal scars used as indication of maturity70. Exact age was determined by assessing tooth growth layer groups (GLG) for a subsample of cases (n = 154), according to previously described methods71. Data can be found in STab. 13.
Health and nutritional status
Porpoises collected for post-mortem investigation were necropsied with the primary aim to determine the animals’ causes of death and their health status, with the quantity and quality of data and results strongly depending on carcass freshness and completeness as well as other logistical and financial factors69. For this study, we established three proxies based on the findings and metrics taken and assessed at necropsy: a proxy for health status and two proxies for nutritional status. For the proxy for health status, we assessed the cause of death among the mature females (n = 199/328) and divided all cases in two categories. In the first category, we placed all mature females which most likely died as a direct result of incidental bycatch (diagnosed based on the presence of encircling imprints or external incisions, the recent ingestion of prey, and the exclusion of other causes of death, for more details see IJsseldijk et al.63), as a direct result of a predatory attack (diagnosed based on the presence of large, sharp-edged mutilations with associated ante-mortem bite lesions bilaterally on the tailstock, extremities or on the head, for more details see Leopold et al.58) or as a direct result of another acute cause, such as sharp forced trauma or dystocia (obstructed labour, full-term foetus) which did not present signs of significant disease or debilitation. All other animals were placed in the second category, with these mature females displaying evidence of general and significant debilitation, including infectious disease (such as significant parasitism, bacterial, viral or mycotic infections) and/or emaciation. Cases that could not be grouped, mostly as a result of decomposition, were excluded from analyses that included this as a parameter.
The first proxy of nutritional status was based on the mean blubber thickness, measured during necropsies in a dorsoventral line on the left body flank just cranial to the dorsal fin, at three locations: dorsal, lateral, and ventral. Blubber thickness in small cetaceans has previously been shown to decrease during periods of fasting72,73 and this metric has been used as proxy of nutrition by others45,74,75,76. However, it should be noted that blubber thickness is not always a good reflection of individual health nor cause of death (e.g., animals dying of acute causes could also be debilitated62,63). There is uncertainty to what extent factors such as age, sex and season naturally influence blubber thickness, and this should be accounted for. Since we focus our analyses on mature females, no further correction for age and sex was done. However, to correct for season, we modelled the mean blubber thickness as a function of Julian date using a generalized additive model (GAM) to allow a smooth effect of the predictor variable (Julian date). This captures the sinus-shaped seasonal variation in blubber thickness which naturally occurs as a result of changing water- and air temperature25,72 (SFig. 5). The residuals of that model were thereby indicative of an adult females’ nutritional status independent of season, and hence they were used as the proxy for nutrition (referred to in the main text as: nutritional status using corBT, Model 1 in STab. 1).
The second proxy of nutritional status used the categorical variable “nutritional condition” (NCC), which is assigned during necropsies as good, moderate, or poor. Animals in good NCC generally presented a convex outline on a cranial perspective, no signs of muscle atrophy (abundant skeletal musculature) and presented signs of visceral fat. Animals in moderate NCC generally did not have a fully round outline on a cranial perspective, showed possible signs of muscle atrophy and did not present visceral fat. Animals in poor NCC generally had a concave outline on cranial perspective, with visible aspects of vertebrae and/or scapula externally, an hollow appearance caudal to the skull and signs of muscle atrophy (based on IJsseldijk et al.69). Since this categorial differentiation is collinear with the first established proxy of nutritional status (SFig. 5), it was not used in the same modelling procedures. Therefore, models were run twice, first with corBT and secondly with NCC (for an overview see STab. 1).
Pregnancy rate and foetus size
The pregnancy rate (PR) was calculated as the proportion of pregnant females in the total sample of mature females (following e.g.,70,77). Pregnancy rates were also calculated separately for the animals in the two different health status categories (see above). To avoid missing the presence of very small, early embryos, samples from the period of conception (June–August23) as well as samples from the period of calving (May–June23) were excluded in the PR calculations. All foetuses were measured during necropsy (of the dam) and a proportion of these were also weighed.
Mean energetic density of diets
As a measure of the quality of prey species constituting the diet of harbour porpoises necropsied in The Netherlands, we calculated the mean energy density of their diet (MEDD). Prey were identified from stomach contents, mostly from otoliths; for each individual prey that could be identified, the fresh mass was estimated (using78 and following29) The energy density (ED) is defined as the energy per kilogram of wet weight of prey8,79. ED values for all prey species encountered were taken from the literature (STab. 7). If for a given prey species no value for ED could be found, the ED of a comparable species (mostly same genus), or the mean value of its family, was used. For species for which multiple ED values were available, values were averaged. ED values reported in kcal were multiplied by 4.184 to convert to kJ (following e.g.80). To calculate the mean ED of the diet for a group of porpoises (MEDD, kJ·g−1, see Table 1) we used:
$$MEDD=frac{1}{sum_{i=1}^{n}{M}_{i}}sum_{i=1}^{n}({M}_{i}*{ED}_{i})$$
(1)
where i is the prey species and M the reconstructed prey mass in grams (following8). The reconstructed prey mass per species is multiplied by the species-specific ED and the energy sum is divided by the total mass of all prey, resulting in the MEDD.
Data analyses and statistical models—Dutch waters
Data were explored prior to analyses following Zuur et al.81,82. Data exploration and analyses were performed using R version 3.6.383, with packages ggplot2, grid, gridExtra, rsq, glmTMB, mgcv and ggpubr. Several statistical models were developed (for referencing in the text see overview in: STab. 1).
Influences on foetus size
To identify which variables influence foetus size, we firstly identified the best measure for foetus size. A Generalized Linear Model (GLM) for foetus length and weight was fitted (Model 2), with weight only available for a subset of all foetuses (n = 34). This model indicated a close relationship between length and weight (R2 of 0.8 for foetus length as a function of mass, SFig. 5), and foetus length was therefore used as representative for foetus size in the subsequent analysis, to increase sample size. GLMs with a Gaussian distribution were used (Model 3). The model selection tested for covariates and their influence on foetus length, with the predictor variables: Julian date to account for foetus length which increases throughout gestation, total length of the mother, health status of the mother, nutritional status of the mother. Interactions between length of the mother and her nutritional status were included following data exploration. Only cases with complete observation of all parameters were included (n = 43). A backwards model selection approach was applied with the drop1 function from the R language used to assess which model terms could be excluded83. The best fitting model was selected using Akaike’s Information Criterion (AIC), which provides a relative measure of the goodness of fit of statistical models. Model validation was done to identify potential violations of model assumptions by inspection of normalized residuals and assessment of residual probability plots. Likelihood profile confidence intervals (95%) and odds ratios of the most optimal model were calculated. Models were run twice, first using the first proxy of nutritional status based on blubber thickness corrected for season (corBT) and secondly using the nutritional condition category (NCC), taken into account blubber thickness, visceral fat and muscle mass (for full descriptions, see above).
Influences on pregnancy
To identify which variables influence pregnancy, we firstly coded all mature, pregnant females as 1 and all mature, non-pregnant females as 0. Next, GLMs with a binomial error distribution and logit link were used (Model 4) to test the influence of included covariates on the likelihood of pregnancy. Only cases with complete observations of all parameters were included (n = 65). The predictor variables included in the saturated model were age, year to assess temporal variance, month to assess seasonal variance, health status (proxy, categorical), and nutritional status. Interactions were added following data exploration: between health and nutritional status, between the health status and year and health status and month. Model selection, validation and interpretation was conducted following the protocol previously described above. Models were run twice, first using the first proxy of nutritional status based on for season corrected blubber thickness (corBT, numerical) and secondly using the nutritional condition category, taking into account blubber thickness, visceral fat and muscle mass (NCC, categorical) (for full descriptions, see above).
Age at sexual maturity
The age at sexual maturity (ASM), or age at 50% maturity, was determined using binomial logistic regression models. Maturity, coded as 1 for mature females and 0 for immature females, was modelled as a function of age (in years) to assess ASM (n = 154, Model 5). The model was fitted using a binomial error distribution and logit link, as is appropriate for binary data and the ASM was estimated by calculating the negative of the slope over the intercept.
Assessment of porpoise life history and environmental condition globally
The life history response variables assessed were PR and ASM, which were obtained from 17 different studies. The earliest study was conducted between 1941 and 1943, but the majority of the studies were performed between 1980 and 2019 (including the present study). The environmental predictor variables used were quality of diet, expressed as mean energy density of diet (MEDD), cumulative human impact (CHI) with data on climate change, fishing, land-based pressures, and other human activities, and lastly chemical pollution expressed as polychlorinated biphenyls (PCBs). Fifteen diet studies were used, ranging from 1985 to 2019 (including this study). One comprehensive study was used to obtain the CHI information for the year 2008. A total of 21 studies reporting PCB levels in harbour porpoises, conducted in the period 1971–2019 (including the present study) were collated. Details below and in STab. 1.
Life history
For PR the following were tabulated: (1) the number of pregnant females out of the total number of mature females in each study, (2) the determined conception period and whether this was accounted for in the calculation of the PR, (3) the method to assess pregnancy, which was either based on the presence of a foetus or presence of a corpora lutea (CL), and (4) the source of the specimens: either directly from fisheries, strandings including trauma cases, or a combination thereof (STab. 3). For ASM we provide: (1) how ASM was assessed in each study, and (2) the standard error (SE) or confidence interval (CI), if reported (STab. 4).
Energy density of prey
A literature search was performed for diet studies from stomach contents of porpoises from or near the study areas where PR and ASM were determined. When multiple diet studies were available the study was selected that best corresponded to the time frame at which PR and ASM were calculated. For the diet studies which reported the reconstructed prey mass in grams we used formula (1) (STab. 8). When the prey mass was reported as a percentage of relative abundance in terms of estimated biomass of prey (%M), we multiplied %M by the ED of the prey species and divided the total %M (STab. 9), using:
$$MEDD=frac{1}{sum_{i=1}^{n}{mathrm{%}M}_{i}}sum_{i=1}^{n}({mathrm{%}M}_{i}*{ED}_{i})$$
(2)
For the studies where the %M was presented in a bar chart, we measured the %M using digital callipers.
Cumulative human impact
An ecosystem-specific, multiscale spatial model containing high resolution data on the intensity of human stressors and their impact on marine ecosystems was developed by Halpern et al.12,14 as part of their Ocean Health Index project. CHI values derived by this model are based on fourteen stressors related to human activities from four primary categories: (1) land-based drivers, including nutrient pollution runoff, organic chemical pollution runoff (pesticides), direct impact of humans (density of coastal human populations), and light; (2) five types of (commercial) fishing, including commercial demersal destructive, commercial demersal non-destructive high bycatch, commercial demersal non-destructive low bycatch, pelagic high bycatch, pelagic low bycatch, and artisanal; (3) climate change, including sea surface temperature, ocean acidification and sea level rise; and (4) shipping. Extensive descriptions of these drivers are published in the methods and supplementary material of Halpern et al.12,14, including information on the origin and validation of the data.
For this study we used the global CHI dataset that is publicly available via the Knowledge Network for Biocomplexity14. Data on CHI was based on the year 2008. We extracted the CHI scores for each of our study areas at ~ 1 km2 resolution and calculated the min, max, mean and median values. To do so, we defined our study areas using the standard georeferenced marine regions as published under the Flanders Marine Institute84. In most cases we combined two or more regions from the database to get full coverage of the study area, but for the study areas where the marine regions did not provide full coverage, we used a manually created polygon. The list of regions is given in STab. 15. For areas with more than one life history study (Denmark and The Netherlands) we used the newest studies since these provided the better match to the time of the CHI score calculation (STab. 2).
Chemical pollution
Polychlorinated biphenyls were not included in the list of organic polluters by Halpern et al.14. However, PCBs have been specifically associated with reproductive impairment in many marine mammal species16,17,18,30,50, therefore the correlation with life history parameters for this industrial organic pollutant was assessed separately. Data was retrieved from the International Whaling Commission’s (IWC) ‘POP Contaminants Trend Explorer’ tool, hosted on the portal of the Sea Mammal Research Unit (SMRU, University of St. Andrews, Scotland). This tool is established under the IWC Scientific Sub-Committee on Environmental Concerns (IWC SC/68A 2019) as part of the IWC Pollution 2020 Initiative and includes data from scientific publications from the 1970s–2000s34,43. The database was provided by the tool manager and included data restricted to adult males, to reduce the bias of biotransfer of chemicals, which occurs during gestation and lactation in females35. The tool reports PCB concentrations in blubber, which is the most commonly assessed tissue in marine mammals for studying the burden of the highly lipophilic and stable PCB compounds47. PCB concentrations that were measured in porpoises in the same areas from which life history parameters were verified with the literature and included. In addition, the literature was searched for PCB analyses of harbour porpoises published in the 2010s, as well as own institutional databases, and data added to align time frame, where possible, with time frame of conducted life history studies (STab. 2).
The presentation of concentrations of pollutants was based on either wet weight (ww) or lipid weight (lw). To allow comparison, the datapoints need to be converted to one common unit, with lw most frequently reported. Studies reporting only ww or dry weight were not included. Studies reporting ww and percentage of lipids (%lipids) were converted to lw, using:
$$lw=frac{ww}{%lipids}*100$$
(3)
The datapoints were converted to mg/kg lw for all studies and the mean ∑TotalPCB is reported per area.
The variance of the sum of congeners reported ranged from ∑6PCBs up to ∑99PCBs, with several older studies reported Aroclor mixtures. Data per congener was however largely not available in literature. We therefore present two mean ∑PCB datapoints: firstly including all studies regardless of the sum of congeners or mixtures (referred to as PCB1), and secondly limited to studies reporting ∑17-99PCBs (referred to as PCB2).
Statistical models for global assessment
For the analyses we restricted to study areas with complete observations of the environmental conditions to compare models. A GLM fitted with a binomial distribution and logit link was used to determine the effect of environmental conditions on pregnancy rates (Model 6). The response variable was the number of pregnant females (Npreg) in the total number of females (Ntotal) (grouped binomial data, STab. 3). Since the differences between study areas can be large because of unknown effects, an individual normal random effect for area was added on the logit scale. Another GLM was conducted to determine the effect of the three environmental conditions on age at sexual maturity (Model 7) fitted with a Gaussian distribution and weighed by sample size (Ntotal) (STab. 4). This model was applied twice using two individual predictor functions: first, with the predictor variables MEDD, CHI and PCB1 and secondly with the predictor variables MEDD, CHI and PCB2. The latter restricted the analyses to a smaller number of study areas due to missing data, but it reduced some of the bias because of very small (< ∑17PCBs) or very large (Aroclor mixture) reported ∑PCB datapoints.
For all models, a backward stepwise model selection process using the drop1 function was conducted and the best models explaining the life history parameters pregnancy rate (Model 8) and age at sexual maturity (Model 9) were identified by assessment of the AIC. The logistic regression model yields log-odds ratios. This represents the odds (number of pregnant per non-pregnant) of being pregnant compared to non-pregnant in each study area, while accounting for differences in sample size. We also obtained the confidence intervals for the estimated log-odds ratios of the most optimal models.
Ethics statement
The animals described in this study were free-living harbour porpoises which died of natural causes. No consent from an Animal Use Committee is required, as the animals described in this study were not used for scientific or commercial testing. Consequently, animal ethics committee approval was not applicable to this work.
Source: Ecology - nature.com
