The Anthropocene has brought substantial change to ocean ecosystems, but whether this age will bring more or less marine disease is unknown. In recent years, the accelerating tempo of epizootic and zoonotic disease events has made it seem as if disease is on the rise. Is this apparent increase in disease due to increased observation and sampling effort, or to an actual rise in the abundance of parasites and pathogens? We examined the literature to track long-term change in the abundance of two parasitic nematode genera with zoonotic potential: Anisakis spp. and Pseudoterranova spp. These anisakid nematodes cause the disease anisakidosis and are transmitted to humans in undercooked and raw marine seafood. A total of 123 papers published between 1967 and 2017 met our criteria for inclusion, from which we extracted 755 host–parasite–location–year combinations. Of these, 69.7% concerned Anisakis spp. and 30.3% focused on Pseudoterranova spp. Meta-regression revealed an increase in Anisakis spp. abundance (average number of worms/ fish) over a 53 year period from 1962 to 2015 and no significant change in Pseudoterranova spp. abundance over a 37 year period from 1978 to 2015. Standardizing changes to the period of 1978–2015, so that results are comparable between genera, we detected a significant 283-fold increase in Anisakis spp. abundance and no change in the abundance of Pseudoterranova spp. This increase in Anisakis spp. abundance may have implications for human health, marine mammal health, and fisheries profitability.

Literature search and data extraction

To identify papers estimating anisakid abundance, we conducted a search in ISI Web of Science. We used the search string: TS = (anisak* or “herring worm” or “herringworm” or Pseudoterranova or whaleworm or “whale worm” or phocanema or “whale-worm”), on October 10th, 2017, which yielded 2,284 papers. We then used a systematic screening process to eliminate non-relevant papers. The screening process was performed by EAF, CAW, and KAD. We first screened titles for suitability, excluding papers that quantified anisakids in humans, birds, or marine mammals and titles that clearly indicated that the paper concerned a different suite of parasite species. After title screening, we retained 1,336 papers for the next stages of the screening process. We then screened abstracts to further focus the dataset, excluding publications that examined anisakids in humans, birds, or marine mammals, that performed experimental manipulation of parasites or hosts, or that were reviews. This process resulted in the retention of 576 papers, which were read in full to determine their eligibility for inclusion in the study. To qualify for final inclusion, papers were required to contain information on host species identity, parasite genus identity, location of host collection (either as a named place or latitude and longitude), collection year or year range, size of the host, how the hosts were examined for anisakids, and either prevalence and intensity of infection with an error estimate or abundance with an estimate of error. Ultimately, we extracted data from 123 papers, resulting in 755 data points (i.e., unique estimates of parasite abundance for a host species in a particular location at a particular time). Of these 755 data points, 69.7% described Anisakis spp. and 30.3% described Pseudoterranova spp.. For each data point, we extracted information on host species identity, collection location, year of collection, portion of the host examined (i.e., all, viscera, musculature), host examination method (i.e., standard visual assessment/candling, microscopy, UV light, acid digestion), parasite genus (i.e., Anisakis or Pseudoterranova), parasite prevalence (i.e., proportion of hosts infected), parasite intensity (i.e., average number of parasites per infected host), error associated with parasite intensity (i.e., standard deviation, range, standard error, or confidence intervals), parasite abundance (i.e., average number of parasites per host for all hosts), and error associated with abundance (i.e., standard deviation, range, standard error, or confidence intervals).

Data standardization

We standardized data prior to analysis. For example, to standardize host length, which was reported in different ways (i.e., standard length, total length, fork length) across studies, we used the standard linear conversion equation,

Length_standard = a + b*Length_reported

We obtained the a and b parameters for each fish species from FishBase (Froese and Pauly 2016), using the average a and b if multiple values were provided by FishBase.

Parasite abundance was quantified as the number of parasites per host, including uninfected hosts (Bush et al. 1997). When parasite abundance was not reported, but parasite intensity and prevalence were, we calculated parasite abundance by multiplying intensity and prevalence and propagating the error through in quadrature (the square root of the sum of squares). If the range of intensity or abundance was provided but the standard deviation was not, we estimated standard deviation by optimizing the negative binomial distribution for the dispersion parameter. We assumed that the maximum abundance or intensity was the 95^th quantile of the negative binomial distribution, since parasites often follow a negative binomial distribution (Shaw et al. 1998). With this assumption, we used the supplied mean intensity or abundance as the mean of the negative binomial distribution with the dispersion parameter unknown. We then used an optimization algorithm to estimate the dispersion parameter that best fit the supplied mean and 95th percentile. With the given mean and estimated dispersion parameter, we were able to calculate the error of the mean given that the variance of a negative binomial distribution is the sum of the mean and the squared mean divided by the dispersion parameter. We also converted other forms of error reported in the studies (i.e., standard error, confidence intervals) back to standard deviation. To convert from standard error to standard deviation, we multiplied the standard error by the square root of the sample size. To convert from a confidence interval, we took the difference from the upper bound of the confidence interval and the mean and then divided the difference by the appropriate z-score and multiplied by the square root of the sample size. For location information, we grouped data points based on major FAO fishing region (FAO 2008) using ESRI ArcGIS (ESRI 2011), to match our data to FAO major fishing region.

Data analysis

To assess change over time in anisakid abundance, we developed two parallel meta-analytic regression models for the two genera of anisakid nematodes contained in our dataset (i.e., Anisakis and Pseudoterranova). For each model, we fourth-root transformed the abundance and standard deviation of the abundance to meet normality assumptions of the model. In our data, we had observations of anisakid abundance with associated error equal to zero. This arose if only a single host was examined, if reported anisakid abundance was zero, or if anisakid abundance was greater than zero, but all fish sampled had the same anisakid burden. Since meta-regression uses the inverse of variance to weight each observation, we corrected the variance to account for observed zero error. We chose to add 1 to every variance to prevent over-weighting (i.e., adding a small value like 0.000001 to every zero-error observation will over-weight these observations). Adding a very small variance correction would give many orders of magnitude more weight to observations with zero variance compared to any observation that had greater than 1 variance. By adding 1 to the variance, we prevent this issue and observations with zero error get full weight while all other observations receive progressively less weight. We included year of collection and host length as moderators in the meta-regression model. We also included host species, portion of the host examined nested in host species, FAO major fishing area (to account for geographic clustering of points), method for detecting parasites, and paper ID as random effects using the rma.mv() function in the package metafor (Viechtbauer 2010) in R. The final models took the form:

(Anisakis_abundance_ijkl)^1/4 or (Pseudoterranova_abundance_ijkl)^1/4 ~

host_length_ijkl +

year_ijkl +

(1 | host_species_ijkl / portion_of_fish_ijkl) +

(1 | FAO_region_ijkl) +

(1 | method_of_detection_ijkl) +

(1 | manuscript_ID_ijkl) ,

where the response variable_ijkl represents a measurement of parasite abundance from the ith study in the jth location at the kth time in the lth host species.

We were also interested in testing which fish species, geographic regions, detection techniques, and portions of host examined drove the patterns observed in the above meta-regression, so we performed four subanalyses. First, we sequentially excluded each host species and reran the model. We extracted an updated estimate of change over time and compared the estimate of the effect of time to the full model that included all host species. This allowed us to determine whether a single host species was driving the patterns we observed. We then performed the same analysis, but instead sequentially excluded each geographic region to determine whether a single geographic region was driving the patterns we observed. To test the influence of detection technique, we excluded studies that used techniques that came into common usage over the period of the study: digestion of fish tissues with acids and use of UV light for parasite detection. We sequentially excluded (1) all studies that used digestion, including digestion combined with other techniques, like microscopy and (2) all studies that used UV, including UV combined with other techniques. We then extracted an updated estimate of change over time and compared the estimate of the effect of time to the full model that included all detection techniques. Finally, to test the influence of portion of host examined, we sequentially excluded all studies that quantified anisakid burden in the (1) alimentary tract, (2) fillet, (3) viscera, and (4) whole body. We then extracted an updated estimate of change over time and compared the estimate of the effect of time to the full model that included all categories of portion of host examined.

This deposition includes several datasets: one raw meta-analytic dataset on the abundance of Anisakis spp. and Pseudoterranova spp. nematodes in studies downloaded from Web of Science (anisakid_data.csv), a dataset specifying the FAO Major Fishing Area of each collection site included in anisakid_data.csv (anisakid_fao.csv), a dataset containing the correct Latin names of all the host species included in anisakid_data.csv (host_data.csv), and two datasets containing the processed data that were used for meta-analysis and which arise from the raw dataset anisakid_data.csv: one for Anisakis spp. (processed_data_anisakis.csv) and one for Psuedoterranova spp. (processed_data_pseudoterranova.csv).

Note that a metadata file is provided for all of the data files (metadata.csv).

Data are available in this repository, and are also permanently archived in GitHub: https://github.com/wood-lab/Fiorenza_et_al_2020_GCB