Ecology

For this study we use the University of Victoria Earth System Climate Model (UVic ESCM) version 2.940,41,42. The UVic ESCM is an intermediate-complexity earth system model with a resolution of 1.8(^circ ) latitude by 3.6(^circ ) longitude and 19 ocean depth levels. The surface ocean level is 50 m deep. The model contains a two dimensional energy moisture-balance model of the atmosphere, as well as representations of sea ice, ocean circulation and sediments, and terrestrial carbon. The particular biogeochemical version used here includes three phytoplankton functional types, namely diazotrophs (DZ), mixed phytoplankton (PH), and small phytoplankton and calcifiers (CO)43. The model pre-industrial climate has been previously described43, as has its response to business-as-usual atmospheric (hbox{CO}_2) forcing44. The following sections describe the MP model. A model schematic is presented in Fig. 6.
Figure 6

Microplastic model schematic. Marine snow is produced as a fixed fraction of the free detritus (DET) pool. MP aggregates with this marine snow, entering the (hbox{MP}_A) (marine snow entrained MP) pool. (hbox{MP}_A) held in aggregates sinks at the aggregate rate, with a fraction reaching the seafloor considered to be lost from the ocean. Detrital remineralisation releases the (hbox{MP}_A) from marine snow aggregates at the rate of detrital remineralisation. MP is also grazed by zooplankton and excreted into a pellet-bound (hbox{MP}_Z) pool. Pellet-bound (hbox{MP}_Z) sinks and is released back to the free MP pool at the rate of detrital remineralisation, but some is also lost at the seafloor. Details on the biogeochemical aspects of the model are previously described43.

Full size image

Model description
The base model43 was modified in order to quantify the roles of two of the three theorised biological export pathways on MP (aggregation in marine snow and zooplankton ingestion; for now, we neglect an explicit representation of biofouling29). We distinguish between detritus that becomes faecal pellets, and the physical aggregation of marine snow, by introducing a new faecal pellet tracer to divert 50% of zooplankton particulate losses into a separate detrital pool27. For simplicity, this new pellet detrital pool has the same sinking parameterisation as the original detritus. Using the same sinking rates for both detrital classes produces ocean biogeochemistry that is identical to the previously published versions of the model. In the model, plastic particles only interact passively with marine snow (they do not, for example, modify aggregate sinking rates), but they interact actively with zooplankton grazing (described below). Plastic particles have been observed to both increase and decrease the sinking rates of marine snow21,22 and decrease the sinking rates of faecal pellets20,38, but for simplicity and as a first approximation we neglect these effects in our model.
Three MP compartments are introduced; “free” (unattached) microplastic (MP), microplastic aggregated in marine snow ((hbox{MP}_A)), and microplastic in zooplankton faecal pellets ((hbox{MP}_Z)). All MP are considered to represent particles within a biologically active size range, but this size (and the particles’ composition) is never made explicit. These assumptions could bias modelled MP towards polymer types favoured by generic zooplankton34 and particle sizes in the lower end33,34 of the defined range of microplastic size. However, our parameter sensitivity testing in the Supplemental Information tests various fractional uptake rates that implicitly consider size and particle composition in how biologically “active” the MP pool is. As with all model ocean tracers, microplastic concentrations (MP) vary according to:

$$begin{aligned} frac{dmathrm {MP}}{dt} = T + S(mathrm {MP}) end{aligned}$$
(1)

With T including all transport terms and S representing all source minus sink terms. The source and sink terms for free microplastic are:

$$begin{aligned} S(mathrm {MP}) = Emis – S(mathrm {MP}_A) – S(mathrm {MP}_Z) + w_p frac{delta mathrm {MP} times F_R}{delta z} end{aligned}$$
(2)

Microplastic is emitted to the ocean (Emis) along coastlines and major shipping routes using a scaling against regional (hbox{CO}_2) emissions (a dataset provided with the standard UVic ESCM version 2.9 package download), in order to approximate degree of industrialisation and population density in this first version of this model. The rate of emission is a proportion of the total annual plastic waste generation ((F_T))45. For now, abiotic degradation of macroplastics as a source of microplastics to the ocean is neglected to keep the model simple and focus on biological transport. The MP then exchanges with the marine snow ((hbox{MP}_A)) and zooplankton faecal pellet ((hbox{MP}_Z)) pools. A fast particle rising rate ((w_p)) of 1.9 cm per second46 is prescribed to a fraction ((F_R)) of the free MP in each grid cell below the surface level as an approximation of positive buoyancy. An alternative approach would be to assign a uniform rise rate to all MP particles, and to subject the value of the rise rate to sensitivity testing. However, a weakness of this alternative approach is that the many types of plastic in the ocean have different characteristic buoyancies, which could produce unique particle pathways18. In this alternative approach it would be more appropriate to explicitly simulate multiple MP types in the model (which we sought to avoid in this first modelling effort for the sake of simplicity). Nevertheless, we conducted a sensitivity test using several different rise rates, and the effect of reducing the mean rise rate was similar to reducing the fraction assigned a rise rate.
In the current model version there are no abiotic breakdown rates (i.e., photo-degradation39) or respiration losses47 removing MP from circulation.
MP is modelled to aggregate in marine snow as:

$$begin{aligned} S(mathrm {MP}_A) = A_{upt} – A_{rel} – w_Dfrac{delta mathrm {MP}_A}{delta z} end{aligned}$$
(3)

MP particles are taken up ((A_{upt})) via a Monod function applied to the rate of marine snow formation (sources of detritus; (D_A) in nitrogen units, multiplied by an aggregation fraction, (F_A)) in order to approximate an increased likelihood of MP/marine snow encounter with increasing MP concentrations that approaches a level of saturation at high MP concentrations:

$$begin{aligned} A_{upt} = frac{mathrm {MP}}{k_P + mathrm {MP}} times source(D_A) times F_A end{aligned}$$
(4)

The uptake constant ((k_P)) is subjected to sensitivity testing, as is the fraction of marine snow aggregation ((F_A)). In this parameterisation, the aggregation of MP in marine snow represents the net uptake of MP into aggregates by both aggregation and biofouling processes. Biofouling occurs mostly in the upper 50 m35, which is the entire surface layer grid cell in our model. The entrainment-release cycle of biofouling is implicit in our parameterisation via the microbial loop, which is temperature-dependent. Sensitivity testing of the (k_P) and (F_A) parameters therefore represent testing of the net aggregation due to non-zooplankton biological aggregation effects. MP is released ((A_{rel})) from marine snow at the rate of detrital remineralisation ((mu _D)). This rate is temperature-dependent and results in higher rates of release in the low latitudes.

$$begin{aligned} A_{rel} = mu _D mathrm {MP}_A end{aligned}$$
(5)

A particle sinking term ((w_D)) applies to marine snow-associated MP, and has the same value as sinking detritus. The base unit of all MP tracers is number of plastic particles. As a first approximation we assume that all marine snow aggregates forming from free detritus have the characteristic of diatom aggregates (8.8 (upmu )g C per aggregate48). Model detritus in mmol N is converted to mmol C using Redfield stoichiometry, which is then converted to (upmu )g C to calculate the maximum number of aggregates. The maximum number of aggregates is then multiplied by the aggregation fraction (F_A), to calculate (hbox{MP}_A) source and sink rates. MP is conserved for all MP tracers when surface flux balances sedimentary loss rate. What fraction of MP particles reaching the seafloor via aggregate and faecal pellet ballasting are returned to the water column ((F_B)) is tested. For simplicity and as a first approximation, detritus ballasted by calcite, and calcite43, are assumed to not aggregate with microplastic.
Similarly, for MP associated with zooplankton, sources and sinks are:

$$begin{aligned} S(mathrm {MP}_Z) = P_{upt} – P_{rel} – w_Dfrac{delta mathrm {MP}_Z}{delta z} end{aligned}$$
(6)

The calculation of MP particle ingestion rate ((P_{upt})) is the same as for other food sources37. A grazing preference ((psi _{MP})) for MP is subjected to sensitivity testing. This sensitivity testing implicitly examines effects such as biofouling altering the grazing preference of zooplankton for MP. It is assumed that 100% of ingested MP will be egested as faecal pellets and released ((P_{rel})) to the “free” MP pool at the rate of faecal pellet remineralisation, with no plastic remaining in the gut and no plastic being metabolised by the zooplankton. Ingesting MP also results in a reduced zooplankton carbon uptake rate19, with implications for primary and export production (although, Redfield ratios are conserved). Pellet-bound (hbox{MP}_Z) is considered to sink at the rate of faecal pellets ((w_D)).
Plastic is eaten by zooplankton in this model. The Holling II grazing formulation37 is extended to include MP. Grazing of MP ((G_{MP})) is calculated as:

$$begin{aligned} begin{aligned} G_{MP}&= mu _Z^{max} times Z times mathrm {MP}times R_{M:P}times R_{F:MP}times R_{N:F}times psi _{MP}\&quad times ,(psi _{CO}CO+psi _{PH}PH+psi _{DZ}DZ+psi _{Detr_{tot}}Detr_{tot}\&quad +,psi_{Z}Z+psi _{MP}mathrm {MP}times R_{M:P}times R_{F:MP}times R_{N:F} + k_Z)^{-1} end{aligned} end{aligned}$$
(7)

The maximum potential grazing rate ((mu _Z^{max})) is scaled by zooplankton population (Z) and MP availability (MP), and weighted by a food preference ((psi _{MP})), total prey (CO, PH, DZ, (hbox{Detr}_{{tot}})), and Z representing the food sources described in44 and a half saturation constant for zooplankton ingestion ((k_z)). Grazing preferences must always sum to 1 in the model, so sensitivity testing of (psi _{MP}) requires that all grazing preferences must also be adjusted. This is done by varying (psi _{MP}) but requiring (psi _{DZ}) always be set to 0.1 (on the basis that diazotrophs are a poor food source, and to minimize disruption to the nitrogen cycle). The remaining allowance is equally split by the other (psi ) terms. The calculation occurs in N units, so MP is first converted to grams of MP using the MP particle-to-mass conversion of 236E3 tonnes MP = 51.2E12 particles MP ((R_{M:P}))4. It is assumed that 1 g MP will roughly replace 1 g of food (at Redfield ratios; (R_{N:F}) is the conversion from mol Food to mol N) in the zooplankton’s diet, and MP is thus converted to mmol N for the grazing calculation. However, we subject this ratio ((R_{F:MP})) to sensitivity testing. Zooplankton uptake of plastic is therefore:

$$begin{aligned} P_{upt} = frac{G_{MP}}{R_{M:P}times R_{F:MP}times R_{N:F}} end{aligned}$$
(8)

MP particles are released from faecal pellets via remineralisation, which occurs at the same rate as the remineralisation of aggregates:

$$begin{aligned} P_{rel} = mu _D mathrm {MP}_Z end{aligned}$$
(9)

Model forcing
The model was integrated at year 1765 boundary conditions (including agricultural greenhouse forcing and land ice) for more than 10,000 years until equilibration was achieved. From year 1765 to 1950, historical (hbox{CO}_2) concentration forcing, and geostrophically adjusted wind anomalies are applied. From 1950 to 2100 the model is forced with a combination of historical (hbox{CO}_2) concentration forcing (to 2000) and a business-as-usual high atmospheric (hbox{CO}_2) concentration projection RCP8.549,50. MP emissions start from 2 million metric tonnes in year 1950 (a total plastic waste generation estimate45), increasing at a rate of 8.4% per year. (hbox{CO}_2) and MP forcing is summarized in Fig. 7. It has been estimated that about 4% of total plastic waste generated enters the ocean30, but that the microplastic mass found at the sea surface represents only about 1% of the annual plastic input to the ocean4. We test a range of input fractions (see Table 2), after applying a mass conversion from tonnes to number of MP particles4. Using a considerable over-estimation of MP pollution rate also implicitly accounts for abiotic degradation of larger plastics.
Figure 7

Model forcing from years 1950–2100. Atmospheric (hbox{CO}_2) follows RCP8.5 (panel a). Plastic flux into the ocean is assumed to be some fraction of the total historical and projected plastic waste generation estimate (panel b), with a continuing rate of increase of 8.4% per year45, converted to MP particles using a mass conversion4. Previous estimates of actual total plastic mass flux into the ocean is only about 4% of the total plastic waste generation30, with the MP fraction being a small proportion of that.

Full size image

Table 2 Microplastic model parameters and range tested.
Full size table

Experimental setup
A 700-member Latin Hypercube51 was used to test the microplastic parameter space of the model using the forcing described in the previous section. While biological model parameters might also influence microplastic uptake and transport, we limited our initial tests to the new parameters introduced above. A range of values was prescribed to the parameters listed in Table 2, in which the parameter space was randomly sampled with a normal distribution. The objective was to see what can be learned about plastic accumulation in the ocean, when very little is known about plastic/particle interactions and basic processes are still poorly understood. An analysis of the full Latin Hypercube parameter search is provided as Supplemental Information.
We adopted an incremental approach to increasing model complexity. We started with a control Hypercube where biology was not allowed to take up plastic, in order to first test the physical parameters ((F_T) and (F_R), the fraction of total annual plastic produced entering the ocean as MP, and the fraction assigned a rise rate, respectively). One hundred simulations were performed in this configuration, with the results analysed in the Supplemental Information. We next included passive plastic aggregation in marine snow (MP plus the (hbox{MP}_A) tracer) in a 300 simulation Hypercube, spread across the (k_P) (marine snow uptake coefficient) parameter space (0–1, 1–100, 100–1000 particles (hbox{m}^{-3}), each with 100 Hypercube simulations) in a normal distribution. These 300 simulations explored the 5 relevant MP model parameters: (F_T), (F_R), (F_A) (marine snow aggregation fraction), (k_P), and (F_B) (fractional return to ocean at the seafloor). These results are also provided in the Supplemental Information. Finally, we added active zooplankton-associated plastic (MP, plus (hbox{MP}_A) and (hbox{MP}_Z) tracers) as a third 300-individual Hypercube set. This third Hypercube is similarly split across the (k_P) parameter space in a normal distribution, but with the addition of grazing parameters (psi _{MP}) (MP grazing preference) and (R_{F:MP}) (the food to MP substitution ratio; 7 parameters in total). More

Meaningless predictors and spatial dependence
The variogram range of the GRF predictors has a strong influence on model validations. Figure 6 shows the increase in the cross-validation (R^2) with increasing ranges for all three study sites. The maximum prediction accuracy is reached when the range of the structurally meaningless predictors (e.g. the GRFs) is similar to the range of the soil properties. This is also the reason why meaningless predictors, with respect to a structural relationship to a soil property, such as photographs of faces or paintings, can produce accurate evaluation statistics.
The effect is more stable when 100 GRFs are used, compared to only 10 GRFs for each range of spatial dependence. That is, when using 100 GRFs, relatively accurate model validations are already reached using lower scale predictors compared to using only 10 GRFs. However, cross-validation accuracy can be relatively high when the variogram ranges of the 10 GRFs are long and thus better resemble the effects of EDFs.
Figure 6

Influence of the variogram ranges of Gaussian Random Fields (GRF) on predictive accuracy. Blue: (R^2) values for models based on 100 GRFs; gold: (R^2) values for models based on 10 GRFs; red: variogram range of the corresponding soil property.

Full size image

We can therefore accept the first two hypotheses: using enough but totally meaningless predictors with similar or longer ranges of spatial dependence than that of the response variable, can result in models that produce high predictive accuracies, however, with zero descriptive accuracy. Both, GRFs and EDFs are spatial but not environmental predictors. So, we can “interpolate” with only a few but very large-scale random predictors. But when we interpolate, we cannot interpret. And when the predictors are not completely linear as EDFs, which is the case for GRFs with random meaningless variations between the sample locations, smoothing between the sample points, which is the concept behind interpolation, is not guaranteed. Hence, using GRFs or otherwise meaningless predictors is neither a structural model nor an interpolation model in the classical sense, and therefore must be rejected in principle.
Reference models and prediction accuracy
The prediction accuracies of the different models are presented in Fig. 7. Figure 8 shows the corresponding maps.
Figure 7

Modelling cross-validation accuracies (represented by the coefficient of determination, (R^2)) for all approaches and datasets (GMS: Gaussian mixed scaling, EDF: Euclidean distance fields, GRF: Gaussian random fields models with 100 or 10 predictors).

Full size image

Figure 8

Modelling results. The first row shows the reference models (GMS: Gaussian mixed scaling, EDF: Euclidean distance fields). The other three rows show the models based on 100 Gaussian random fields (GRF) with different variogram ranges.

Full size image

Generally, the GMS models produced the best results. This could be an effect of the number and structure of the predictors, which is larger compared to the merged restricted GMS + EDF dataset, and thus possibly an effect of overfitting to the dependence structure of the data23.
Models with EDFs generally performed better than models with the restricted GMS. Hence, the spatial dependence of soil properties cannot be described by the constrained GMS dataset. There is an increase in prediction accuracy when using the EDFs together with the restricted GMS dataset (except for the Meuse dataset), which could be an indicator for non-stationarity.
In terms of prediction accuracy, modelling with GRFs produced relatively high predictive accuracies (Fig. 7).
The information horizon—descriptive accuracy and relevance
Generally, the descriptive accuracy is strong if there is a causal relationship25 between the response and a predictor. The descriptive accuracy can also be high if there are associations among variables as usually inferred by statistical analysis, which can suggest potential causal relationships1.
Provided that the relevance of the predictors is given and the algorithm or method to generate the data is valid15, 17, the results of this study indicate that the descriptive accuracy is high if the range of spatial dependence of a predictor is equal or smaller than the range of the response variable. However, predictive accuracy can increase when predictors with longer ranges than the range of the environmental property are included, possibly due to non-stationarities in the environmental process4,26 or effects of anisotropy. Therefore, one should remove only those predictors from multiscale approaches15,17,27,28,29,30 with variogram ranges that are long with respect to the size of the study area and if their information content is below a certain minimum (Fig. 4).
If, on the other hand, predictors show ranges of about the diagonal length of the entire study area, then they resemble the properties of the EDFs (Fig. 5) and their ranges are too long to lend themselves to interpretation. In these cases, predictors behave indistinguishably from purely spatial predictors, although they might still be interpretable in some situations.
In summary, these results confirm our third hypothesis, and show that the primary information horizon is located somewhere between the range of the variogram of the soil property and a certain minimal variation of the predictors across the study site.
Beyond the information horizon—descriptive uncertainty and contextual complexity
We recently showed that when finer to coarser scales are successively removed from a set of all scales of a GMS modelling, prediction accuracy usually remains high, even if only the coarsest scales remain in the model4. In cases where the prediction accuracy decreases, we can assume that (i) structural information is lacking, that (ii) interpolation is not the appropriate method, or that (iii) the spatial dependence of the coarse scale GMS predictors are not suitable for interpolation, e.g. if all these coarse scale predictors only show a trend in one direction, for example in X direction only instead of X and Y direction.
We also found an increase in prediction accuracy beyond the range of the variogram of soil properties in GMS and similar approaches4,15,26, when successively adding coarser scales. There are two explanations. First, not all original terrain properties show exactly the same original scale or range, which is due to the convolution functions and the general approaches to calculate terrain properties (e.g. first and second order derivatives, i.e. slope and curvature). Second, this effect might be related to non-stationarity, where coarse-scale predictors can help to “divide” the study area into zones. We tested this here by combining the restricted set of GMS predictors, which are within the information horizon, with EDFs. In all cases prediction accuracy increased. Hence, there is obviously some spatial dependence present, resulting from predictor interactions on very coarse scales or long ranges, which are beyond the information horizon, inferable by the size of the study area. In these cases, there will be some uncertainty in the descriptive accuracy when using predictors that show information contents below a certain minimum of spatial variation.
Looking specifically at the three study sites we see complex soil property formation processes due to interactions of predictors at different scales. This contextual complexity has to be taken into account when interpreting environmental predictors beyond the information horizon, as discussed above.
In Piracicaba very coarse-scale predictors are important27. The soil formation system, however, is rather simple. It is based on rock formation, strike and dip, and subsequent erosion. In this case coarse scale terrain indicators for aspect are good proxies to differentiate between the two different types of parent material, even though they resemble properties of EDFs. In such cases partial dependence models should be applied to aide interpretation27.
The silt content in Rhine-Hesse is controlled by local silt translocation31, which occurred in the last glacial period of the Pleistocene epoch (Würm glaciation) and which was modulated by interactions of climate and terrain. This can be described in terms of a teleconnected system32 and can be mapped by terrain only, which then serves as a proxy for that system26,27. Similar to Piracicaba, interpretations of predictors with very large ranges and relatively low information content can be reasonable. However, the descriptive uncertainty is higher compared to predictors that fall within the information horizon.
The situation for Meuse is different due to a different dominant process system. The zinc content is driven by flooding events. Therefore, different and more relevant predictors, such as the distance to the river Meuse, should be used in this case9. We see that EDFs perform better compared to the mixed dataset (GMS restricted + EDF). This shows that the multiscale terrain predictors are not relevant, but represent noise, and can therefore serve at most as vague proxies. Another problem resulting in such an effect could be algorithmic issues related to feature selection within the Random Forests model, which in some specific cases might occur in relation to autocorrelated predictors33, or to effects due to fitting noise5.
Interestingly, in all cases the GMS models perform better compared to the mixed dataset (GMS restricted + EDF). This can be either due to a higher number of predictors in the GMS approach or relevant structural predictors beyond the information horizon.
Generally, the interpretation of environmental predictors beyond the edge of the information horizon needs specific care and is afflicted with more uncertainty. More

In this study, we generated a new time-series dataset to investigate the response of dark-synchronized MT+ cultures to SIP− filtrate (6 time points, 0–9 h, Fig. 2a) and to compare expression patterns with two existing datasets on the response of MT− to SIP+ [16, 18] (4 time points, 0–10 h, Fig. 2b). Here, we first report the results of a conventional DE analysis for each dataset separately. Next, we use these results to discover genes or biological processes related to sexual reproduction, either based on their predicted functional annotation or the current literature. Finally, we integrate the results of both the novel and publicly available datasets in an integrative analysis that aims at identifying key genes involved in either a single or both mating types.
Fig. 2: Transcriptional responses induced by SIP treatment of Seminavis robusta.

a Multidimensional scaling (MDS) plot for MT+ expression data (0 h–9 h), and (b) MDS plots of two MT− expression datasets (0 h–3 h and 10 h, respectively). Distances between samples in the MDS plot approximate the log2 fold changes of the top 500 genes. c Number of significant DE genes between control and SIP treated cultures for each time point in both mating types. Each dataset was analyzed on a 5% overall FDR (OFDR) level, i.e., the fraction of false positive genes over all rejected genes. The color represents mating type (MT+: red, MT−: blue) and the shade denotes direction of change (dark: downregulated after SIP treatment, light: upregulated after SIP treatment). d Number of significantly enriched GO terms on a 5% FDR level for each time point. The color represents mating type (MT+: red, MT−: blue). A high number of GO terms are discovered in the early time points for MT+, in comparison to the number of DE genes.

Full size image

General transcriptional response and identification of key SIP responsive genes
Multidimensional scaling plots of the RNA-seq data (Fig. 2a, b) showed that in both mating types the dark-to-light transition and time since illumination were the major drivers of gene expression change throughout the experiments. However, as time progresses, the effect of SIP becomes more pronounced. This is supported by DE analysis showing that the number of significant genes increased markedly at later time points (Fig. 2c). Overall, most genes are significantly DE in only a single time point on a 5% overall FDR (OFDR) level (Supplementary Fig. 1B). Combined, on a 5% OFDR level, 4037 genes are DE in response to SIP treatment for MT+, while 5486 genes are found to be DE in MT− in the first 3 h and 6079 genes after 10 h. The stronger response during the first 3 time points in MT− versus MT+ may be the result of the use of a chromatographic fraction of SIP+ for MT− while MT+ cultures were treated with a SIP− containing filtrate. Gene sets for many biological processes were significantly enriched in the conventional lists of DE genes: a total of 1081 and 740 enriched biological process terms were discovered for the response to sex pheromones in MT+ and MT−, respectively (Fig. 2d, Supplemenatry Fig. 2).
We discovered key genes by developing a statistical integrative analysis workflow that is capable of testing for equivalent, i.e., non-DE, expression between conditions. Coupling equivalence testing in one mating type with DE calls in the other allowed for the discovery of key genes exhibiting mating type specific responses to SIP, while DE calls in both datasets found key genes responsive in both mating types. This workflow revealed 52 key genes responding to SIP in both mating types (SRBs), 12 genes uniquely responding in MT+ (SRPs) and 70 genes uniquely responding in MT− (SRMs) (Fig. 3a, Supplementary Figs. 3, 4, 5). Similar to the conventional DE analysis, the response of MT− was more pronounced compared to MT+, likely due to technical differences such as different protocols for pheromone administration. Remarkably, while in MT− we discovered a comparable number of down- and upregulated SRMs, we only found upregulated SRPs and SRBs, possibly indicating that sexual processes induced by SIPs are mainly driven by the induction of key genes rather than the downregulation of inhibitory genes (Fig. 3a).
Fig. 3: Visualization and main results of the integrative workflow.

a Schematic representation of the integrative workflow indicating how SIP responsive genes with a shared response (SRBs) or mating type specific response (SRPs, SRMs) were discovered. Non-responsive genes consist of genes that are equivalently expressed after SIP treatment versus control, or that are very lowly/not expressed (filtered). A log fold change (LFC) cutoff of ±log(3) was used to define responsive (differentially expressed) genes and equivalent genes. At the right side of the panel, log2 fold changes of SRMs, SRPs and SRBs in both mating types are plotted. Each gene is plotted for the time point at which they are differentially expressed. Genes which were not expressed (“filtered”) in the non-responsive mating type are plotted as diamonds. The red horizontal lines represent the fold change cutoff used to determine equivalence and differential expression. The number of discovered genes is indicated in the top left corner of each plot. b Expression of a selection of SIP responsive genes (SRMs, SRPs, SRBs). For each gene, counts per million (CPM) are plotted as a function of time for both mating types. The data points correspond to gene expression of the replicates in each time point and the solid line represents the mean. Data points and lines are colored according to condition, i.e., black for control condition and orange for SIP treatment. c Expression of the five genes belonging to the gene family of SRP12 (Sro2882_g339270). Data are presented in the plots in the same manner as in (b).

Full size image

In what follows we will discuss in more detail the genes and pathways that are responding to SIP in both mating types or uniquely in one mating type and link these changes to physiological events in the mating process. In each section, we first discuss key genes highlighted by the integrative analysis (Fig. 3), after which we discuss results from the conventional DE analysis, focussing on selected biological processes (Fig. 4).
Fig. 4: Gene expression of Seminavis robusta genes involved in mating-related processes.

a Heatmap of genes related to mitotic and meiotic cell cycle progression, which are differentially expressed (DE) in both mating types in the conventional DE analysis. Each gene is plotted for control and SIP treated conditions in both S. robusta mating types. Genes are specified as row names and are scaled relative to the mean expression, amounting to counts per million (CPM) standardized to zero mean and unit variance for each gene in each mating type separately. Blue indicates low expression, while red indicates high expression. b Expression of genes related to diproline synthesis and reactive oxygen species (ROS) production, which are significantly DE in only one mating type in the conventional DE analysis. CPM are plotted as a function of time for both mating types. The data points correspond to gene expression of the replicates in each time point while the solid line represents the mean. Data points and lines are colored according to condition, i.e., black for control condition and orange for SIP treatment. P5CS = Δ1-pyrroline-5-carboxylate synthetase; P5CR = Δ1-pyrroline-5-carboxylate reductase.

Full size image

Responses to SIP conserved in both mating types
Integrative analysis reveals key genes responsive in both mating types
A large fraction (22/52) of SRBs, i.e., key genes with a strong response to SIP in both mating types, lack any functional annotation and homology to sequenced genomes of other diatoms (Supplementary Table 2), suggesting that the molecular mechanisms underlying early mating are highly species-specific. The remaining 30 SRBs can be linked to energy metabolism, ROS signaling and meiosis, amongst others. Pyruvate kinase (Sro373_g129070) and isocitrate dehydrogenase (Sro492_g153950), respectively involved in glycolysis and the citric acid cycle, are strongly upregulated in both mating types (Fig. 3b), suggesting an increased energy demand. Interestingly, pyruvate kinase is also upregulated during gametogenesis in the brown alga Saccharina latissima [26] and the parasite Plasmodium berghei [27]. In addition, two enzymes from the pentose phosphate pathway (PPP) are among the SRBs: transketolase (Sro524_g159900) and transaldolase (Sro196_g083630) (Supplementary Fig. 3). The PPP generates NADPH, a reductive compound needed in various metabolic reactions and involved in detoxification of ROS by regenerating glutathione [28, 29]. Furthermore, one SRB encoding a heme peroxidase (Sro1252_g256250) exhibited strong upregulation upon SIP treatment (Fig. 3b). Upregulation of heme peroxidases was also reported during sexual reproduction in other eukaryotes, e.g., mosquitoes (Anopheles gambiae) [30] and fungi [31, 32]. Heme peroxidases promote substrate oxidation in various metabolic pathways and are essential for the detoxification of ROS [33], suggesting that ROS signaling plays a role in the response to SIP, as seen in the green algae Volvox carteri, where high ROS levels trigger sex [34]. Finally, a highly expressed SRB encodes a transmembrane protein containing an Epidermal Growth Factor-like (EGF-like) domain (Sro65_g036830, Fig. 3b), with potential orthologs encoded in pennate and centric diatoms including P. tricornutum and T. pseudonana (BLASTp, E  More

1.
Brazelton, C. An overview of global blueberry production in 2016. World Blueberry Statistics & Intelligence Report Global Summary. International Blueberry Organization. Preprint at https://www.internationalblueberry.org/library/ (2017).
2.
Farneti, B. et al. Exploring blueberry aroma complexity by chromatographic and direct-injection spectrometric techniques. Front. Plant Sci. 8, 617 (2017).
PubMed  PubMed Central  Article  Google Scholar 

3.
Anand, S. P. & Sati, N. Artificial preservatives and their harmful effects: Looking toward nature for safer alternatives. Int. J. Pharm. Sci. Res. 4, 2496–2501 (2013).
Google Scholar 

4.
Bennett, J. W. & Klich, M. Mycotoxins. Clin. Microbiol. Rev. 16, 497–516 (2003).
CAS  PubMed  PubMed Central  Article  Google Scholar 

5.
Ochmian, I., Kozos, K., Chelpinski, P. & Szczepanek, M. Comparison of berry quality in highbush blueberry cultivars grown according to conventional and organic methods. Turk. J. Agric. For. 39, 174–181 (2015).
CAS  Article  Google Scholar 

6.
Ochmian, I., Oszmiański, J., Jaśkiewicz, B. & Szczepanek, M. Soil and highbush blueberry responses to fertilization with urea phosphate. Folia Hort. 30, 295–305 (2018).
Article  Google Scholar 

7.
Serna-Jimenez, J. A., Quintanilla-Carvajal, M. X., Rodriguez, J. M., Uribe, M. A. & Klotz, B. Development of a combined temperature and pH model and the use of bioprotectants to control of Mucor circinelloides. Am. J. Food. Technol. 11, 21–28 (2016).
CAS  Article  Google Scholar 

8.
Mijowska, K., Ochmian, I. & Oszmiański, J. Rootstock effects on polyphenol content in grapes of ‘Regent’ cultivated under cool climate condition. J. Appl. Bot. Food Qual. 90, 159–164 (2017).
CAS  Google Scholar 

9.
Caruso, F. L. & Ramsdell, D. C. Compendium of blueberry and cranberry diseases. 1st ed. ISBN-0890541736 (American Phytopathological Society, February 15, 1995).

10.
Koyshibayev, M. & Muminjanov, H. Monitoring Diseases, Pests, and Weeds in cereal crops. Food and Agriculture Organizationof the United Nations. ISBN 978-92-5-109180-7. Preprint at https://www.fao.org/3/a-i5550e.pdf (2016).

11.
Drusch, S. & Ragab, W. Mycotoxins in fruits, fruit juices and dried fruits. J. Food Prot. 66, 1514–1527 (2003).
CAS  PubMed  Article  Google Scholar 

12.
Fernández-Cruz, M. L., Mansilla, M. L. & Tadeo, J. L. Mycotoxins in fruits and their processed products: Analysis, occurrence and health implications. J. Adv. Res. 1, 113–122 (2010).
Article  Google Scholar 

13.
Tournas, V. & Katsoudas, E. Mould and yeast flora in fresh berries, grapes and citrus fruits. Int. J. Food Microbiol. 105, 11–17 (2005).
CAS  PubMed  Article  Google Scholar 

14.
Barrau, C., de los Santos, B. & Romero, F. Susceptibility of southern highbush and rabbiteye blueberry cultivars to postharvest diseases in Huelva, Spain. Acta Hortic. 715, 525–530 (2006).

15.
Wright, E. R., Rivera, M. C., Esperon, J., Cheheid, A. & Codazzi, A. R. Alternaria leaf spot, twig blight, and fruit rot of highbush blueberry in Argentina. Plant Dis. 88, 1383–1383 (2004).
CAS  PubMed  Article  Google Scholar 

16.
Frąc, M., Hannula, S. E., Bełka, M. & Jędryczka, M. Fungal biodiversity and their role in soil health. Front. Microbiol. 9, 707 (2018).
PubMed  PubMed Central  Article  Google Scholar 

17.
Bensch, K. et al. Common but different: The expanding realm of Cladosporium open access. Stud. Mycol. 82, 23–74 (2015).
CAS  PubMed  PubMed Central  Article  Google Scholar 

18.
Hole, D. G. et al. Does organic farming benefit biodiversity?. Biol. Conserv. 122, 113–130 (2005).
Article  Google Scholar 

19.
Erisman, J. W. et al. Agriculture and biodiversity: A better balance benefits both. AIMS Agric. Food. 1, 157–174 (2016).
Article  Google Scholar 

20.
Kućmierz, J., Nawrocki, J. & Sojka, A. Fungi isolated from diseased early green fruits and fruits of blueberry (Vaccinium corymbosum L.). Prog. Plant Prot. 53, 779–784 (2013).
Google Scholar 

21.
Teillard, F. et al. A review of indicators and methods to assess biodiversity – Application to livestock production at global scale. Livestock Environmental Assessment and Performance (LEAP) Partnership. FAO, Rome, Italy. Preprint at https://www.fao.org/3/a-av151e.pdf (2016).

22.
Sheoran, H. S., Phogat, V. K., Dahiya, R. & Gera, R. Long term effect of farming practiceson microbial biomass carbon, enzyme activities and microbial populations in different textured soils. Appl. Ecol. Environ. Res. 16, 3669–3689 (2018).
Article  Google Scholar 

23.
Bagyaraj, D. J. & Ashwin, R. Soil biodiversity: Role in sustainable horticulture. In Biodiversity in Horticultural Crops (ed. Peter, K. V.) 5–6 (Daya Publishing House, New Delhi, 2017).
Google Scholar 

24.
Reboud, X., Eychenne, N., Délos, M. & Folcher, L. Withdrawal of maize protection by herbicides and insecticides increases mycotoxins contamination near maximum thresholds. Agron. Sustain. Dev. 36, 43 (2016).
Article  CAS  Google Scholar 

25.
Benbrook, C. M. Trends in glyphosate herbicide use in the United States and globally. Environ. Sci. Eur. 28, 3 (2016).
PubMed  PubMed Central  Article  CAS  Google Scholar 

26.
Alshannaq, A. & Уu, J.-H. Occurrence, toxicity, and analysis of major mycotoxins in food. Int. J. Environ. Res. Public Health. 14, 632 (2017).
PubMed Central  Article  CAS  PubMed  Google Scholar 

27.
Petruzzi, L. et al. In vitro removal of ochratoxin a by two strains of Saccharomyces cerevisiae and their performances under fermentative and stressing conditions. J. Appl. Microbiol. 116, 60–70 (2014).
CAS  PubMed  Article  Google Scholar 

28.
Zhu, R. et al. Detoxification of mycotoxin patulin by the yeast Rhodosporidium paludigenum. Food Chem. 179, 1–5 (2015).
CAS  PubMed  Article  Google Scholar 

29.
Chi, Z. et al. Bioproducts from Aureobasidium pullulans, a biotechnologically important yeast. Appl. Microbiol. Biotechnol. 82, 793–804 (2009).
CAS  PubMed  Article  Google Scholar 

30.
Gindro, K. & Pezet, R. Purification and characterization of a 40.8-kDa cutinase in ungerminated conidia of Botrytis cinerea Pers: Fr. FEMS Microbiol. Lett. 171, 239–243 (1999).
CAS  PubMed  Article  Google Scholar 

31.
Jiang, H., Sun, Z., Jia, R., Wang, X. & Huang, J. Effect of Chitosan as an antifungal and preservative agent on postharvest blueberry. J. Food Qual. 39, 516–523 (2016).
CAS  Article  Google Scholar 

32.
Ochmian, I. & Kozos, K. Influence of foliar fertilisation with calcium fertilisers on the firmness and chemical composition of two highbush blueberry cultivars. J. Elem. 20, 185–201 (2015).
Google Scholar 

33.
Kozos, K. & Ochmian, I. The influence of fertilisation urea phosphate on growth and yielding bush of two highbush blueberry cultivars (V. corymbosum). Folia Pomer. Univ. Technol. Stetin., Agric., Aliment., Pisc Zootech. 37, 29–38 (2016).
Article  Google Scholar 

34.
Ochmian, I., Grajkowski, J. & Skupień, K. Influence of substrate on yield and chemical composition of highbush blueberry fruit cv. ‘Sierra’. J. Fruit Ornam. Plant Res. 17, 89–100 (2009).
CAS  Google Scholar 

35.
Jiang, Y., Li, Y., Zeng, Q., Wei, J. & Yu, H. The effect of soil pH on plant growth, leaf chlorophyll fluorescence and mineral element content of two blueberries. Acta Hortic. 1180, 269–276 (2017).
Article  Google Scholar 

36.
Pospieszny, H. Systemic Acquired Resistance (SAR) in integrated plant protection. Prog. Plant Prot. 56(4), 436–442 (2016).
Google Scholar 

37.
Fallahi, E., Conway, W. S., Hickey, K. D. & Sams, C. E. The role of calcium and nitrogen inpostharvest quality and disease resistance of apples. HortScience 32, 831–835 (1997).
CAS  Article  Google Scholar 

38.
Malakouti, M. J., Tabatabaei, S. J., Shahabil, A. & Fallahi, E. Effects of calcium chloride on apple fruit quality of trees grown in calcareous soil. J. Plant. Nutr. 22, 1451–1456 (1999).
CAS  Article  Google Scholar 

39.
Ochmian, I. et al. The feasibility of growing highbush blueberry (V. corymbosum L.) on loamy calcic soil with the use of organic substrates. Sci. Hortic. 257, 108690 (2019).
CAS  Article  Google Scholar 

40.
Chiabrando, V. & Giacalone, G. Anthocyanins, phenolics and antioxidant capacity after fresh storage of blueberry treated with edible coatings. Int. J. Food Sci. Nutr. 66, 248–253 (2015).
CAS  PubMed  Article  Google Scholar 

41.
Chiabrando, V., Peano, C., Beccaro, G., Bounous, G. & Rolle, L. Postharvest quality of highbush blueberry (Vaccinium corymbosum L.) cultivars in relation to storage methods. Acta Hort. 715, 545–551 (2006).
Article  Google Scholar 

42.
Shiono, M., Matsugaki, N. & Takeda, K. Structure of the blue cornflower pigment. Nature 436, 791–791 (2005).
CAS  PubMed  Article  ADS  Google Scholar 

43.
Ścibisz, I., Kalisz, S. & Mitek, M. Thermal degradation of anthocyanins in blueberry fruit. Zywnosc Nauka Technologia Jakosc 17, 56–66 (2010).

44.
Seeram, N. P. Berry fruits: Compositional elements, biochemical activities, and the impact of their intake on human health, performance, and disease. J. Agric. Food Chem. 56, 627–629 (2008).
CAS  PubMed  Article  Google Scholar 

45.
Wolfe, K. L. et al. Cellular antioxidant activity of common fruits. J. Agric. Food Chem. 56, 8418–8426 (2008).
CAS  PubMed  Article  Google Scholar 

46.
Ribera, A. E., Reyes-Diaz, M., Alberdi, M., Zuñiga, G. E. & Mora, M. L. Antioxidant compounds in skin and pulp of fruits change among genotypes and maturity stages in highbush blueberry (Vaccinium corymbosum L.) grown in southern Chile. J. Soil Sci. Plant Nutr. 10, 509–536 (2010).
Article  Google Scholar 

47.
Wojciechowska, R. Accumulation of nitrates and quality of horticultural products. Pub. Coperite, Kraków, Poland, 21–27 (2005).

48.
Rozek, S. Factors affecting the accumulation of nitrate in vegetable yields. Zesz. Nauk. AR Krak. Ses. Nauk. 71, 19–31 (2000).
Google Scholar 

49.
Regulation (EC) no 396/2005 of the European Parliament and of the Council of 23 February 2005 on maximum residue levels of pesticides in or on food and feed of plant and animal origin andamending Council Directive 91/414/EEC -OJ L70. 16.3.2005 (2005).

50.
Connor, A. M., Luby, J. J., Hancock, J. F., Berkheimer, S. & Hanson, E. J. Changes in fruit antioxidant activity among blue-berry cultivars during cold-temperature storage. J. Agric. Food Chem. 50, 893–898 (2002).
CAS  PubMed  Article  Google Scholar 

51.
Rodrigues, E. et al. Phenolic compounds and antioxidant activity of blueberry cultivars grown in Brazil. Food Sci. Technol. 31, 911–917 (2011).
Article  Google Scholar 

52.
Koca, I. & Karadeniz, B. Antioxidant properties of blackberry and blueberry fruits grown in the Black Sea Region of Turkey. Sci. Hortic. 121, 447–450 (2009).
CAS  Article  Google Scholar 

53.
McDougall, G. J. et al. Different polyphenolic components of soft fruits inhibit α-amylase and α-glucosidase. J. Agric. Food Chem. 53, 2760–2766 (2005).
CAS  PubMed  Article  Google Scholar 

54.
Lachowicz, S., Wiśniewski, R., Ochmian, I., Drzymała, K. & Pluta, S. Anti-microbiological, anti-hyperglycemic and anti-obesity potency of natural antioxidants in fruit fractions of saskatoon berry. Antioxidants. 8, 397 (2019).
CAS  PubMed Central  Article  PubMed  Google Scholar 

55.
Johnson, M. H., Lucius, A., Meyer, T. & de Mejia, E. G. Cultivar evaluation and effect of fermentation on antioxidant capacity and in vitro inhibition of α-amylase and α-glucosidase by highbush blueberry (Vaccinium corombosum). J. Agric. Food Chem. 59, 8923–8930 (2011).
CAS  PubMed  Article  Google Scholar 

56.
Hodges, D. M., Lester, G. E., Munro, K. D. & Toivonen, P. M. Oxidative stress: Importance for postharvest quality. Hort. Sci. 39, 924–929 (2004).
CAS  Article  Google Scholar 

57.
Okan, O. T. et al. Antioxidant activity, sugar content and phenolic p rofiling of blueberries cultivars: A comprehensive comparison. Not. Bot. Horti. Agrobot. Cluj. Napoca. 46, 639–652 (2018).
CAS  Article  Google Scholar 

58.
Raghvendra, V. S. et al. Chemical and potential aspects of anthocyanins—A water-soluble vacuolar flavonoid pigments: A review. Int. J. Pharm. Sci. Rev. Res. 6, 28–33 (2011).
CAS  Google Scholar 

59.
Silva, S., Costa, E. M., Coelho, M. C., Morais, R. M. & Pintado, M. E. Variation of anthocyanins and other major phenolic compounds throughout the ripening of four Portuguese blueberry (Vaccinium corymbosum L) cultivars. Nat. Prod. Res. 31, 93–98 (2017).
CAS  PubMed  Article  Google Scholar 

60.
Carbonaro, M., Mattera, M., Nicoli, S., Bergamo, P. & Cappelloni, M. Modulation of antioxidant compounds in organic vs conventional fruit (peach, Prunus persica L., and pear, Pyrus communis L.). J. Agric. Food Chem. 50, 5458–5462 (2002).
CAS  PubMed  Article  Google Scholar 

61.
Reque, P. M. et al. Cold storage of blueberry (Vaccinium spp.) fruits and juice: anthocyanin stability and antioxidant activity. J. Food Compos. Anal. 33, 111–116 (2013).
Article  CAS  Google Scholar 

62.
Schaefer, H. M., Rentzsch, M. & Breuer, M. Anthocyanins reduce fungal growth in fruits. Nat. Prod. Commun. 3, 1267–1272 (2008).
CAS  Google Scholar 

63.
Huang, H. T. Decolorization of anthocyanins by fungal enzymes. J. Agric. Food Chem. 3, 141–146 (1955).
CAS  Article  Google Scholar 

64.
Keppler, K. & Humpf, H. U. Metabolism of anthocyanins and their phenolic degradation products by the intestinal microflora. Bioorg. Med. Chem. 13, 5195–5205 (2005).
CAS  PubMed  Article  Google Scholar 

65.
Miyahara, M. et al. Potential of aerobic denitrification by Pseudomonas stutzeri TR2 to reduce nitrous oxide emissions from wastewater treatment plants. Appl. Environ. Microbiol. 76, 4619–4625 (2010).
CAS  PubMed  PubMed Central  Article  Google Scholar 

66.
Borowik, A., Wyszkowska, J., Gałązka, A. & Kucharski, J. Role of Festuca rubra and Festuca arundinacea in determinig the functional and genetic diversity of microorganisms and of the enzymatic activity in the soil polluted with diesel oil. Environ. Sci. Pollut. R. 23, 1–14 (2019).
Google Scholar 

67.
Connor, A. M., Luby, J. J. & Tong, C. Genotype and environmental variation in antioxidant activity, total phenolic content, and anthocyanin content among blueberry cultivars. J. Am. Soc. Hortic. Sci. 127, 89–97 (2002).
CAS  Article  Google Scholar 

68.
Krupa, T. & Tomala, K. Antioxidant capacity, anthocyanin content profile in ‘Bluecrop’ blueberry fruit. Veg. Crops Res. Bull. 66, 129–141 (2007).
Google Scholar 

69.
Kozos, K., Ochmian, I. & Chełpiński, P. The effects of rapid chilling and storage conditions on the quality of Brigitta Blue cultivar highbush blueberries (Vaccinium corymbosum L.). Folia Hortic. 26, 147–153 (2014).
Article  Google Scholar 

70.
Anonymous. Microbiology of food and animal feeding stuffs – Horizontal method for the enumeration of yeasts and moulds. Part 1: Colony count technique in products with water activity greater than 0.95. ISO 21527-1:2008. 1st edn (International Organization for Standardization, Geneva, Switzerland, 2008).

71.
Ulbin-Figlewicz, N., Jarmoluk, A. & Marycz, K. Antimicrobial activity of low-pressure plasma treatment against selected foodborne bacteria and meat microbiota. Ann. Microbiol. 65, 1537–1546 (2015).
CAS  PubMed  Article  Google Scholar 

72.
Sandhya Deepika, D. & Laxmi Sowmya, K. Bioindices of bacterial communities. Int. J. Curr. Microbiol. App. Sci. 5, 219–233 (2016).
Article  Google Scholar 

73.
Ochmian, I., Oszmiański, J., Lachowicz, S. & Krupa-Małkiewicz, M. Rootstock effect on physico-chemical properties and content of bioactive compounds of four cultivars Cornelian cherry fruits. Sci. Hortic. 256, 108588 (2019).
CAS  Article  Google Scholar 

74.
Podsedek, A., Majewska, I., Redzynia, M., Sosnowska, D. & Koziołkiewicz, M. In vitro inhibitory effect on digestive enzymes and antioxidant potential of commonly consumed fruits. J. Agric. Food Chem. 62, 4610–4617 (2014).
CAS  PubMed  Article  Google Scholar 

75.
Nickavar, B. & Уousefian, N. Evaluation of α-amylase inhibitory activities of selected antidiabetic medicinal plants. J. Verbrauch Lebensm. 6, 1915 (2011).
Article  Google Scholar 

76.
Arnao, M. B., Cano, A. & Acosta, M. The hydrophilic and lipophilic contribution to total antioxidant activity. Food Chem. 73, 239–244 (2001).
CAS  Article  Google Scholar 

77.
Brand-Williams, W., Cuvelier, M. E. & Berset, C. Use of a free radical method to evaluate antioxidant activity. LWT Food Sci. Technol. 28, 25–30 (1995).
CAS  Article  Google Scholar 

78.
Mijowska, K., Ochmian, I. & Oszmiański, J. Impact of cluster zone leaf removal on grapes cv. regent polyphenol content by the UPLC-PDA/MS method. Molecules 21, 1688 (2016).
PubMed Central  Article  CAS  PubMed  Google Scholar 

79.
Błajet-Kosicka, A., Twarużek, M., Kosicki, R., Sibiorowska, E. & Grajewski, J. Co-occurrence and evaluation of mycotoxins in organic and conventional rye grain and products. Food Control 38, 61–66 (2014).
Article  CAS  Google Scholar  More

1.
IPCC. Special Report: The Ocean and Cryosphere in a Changing Climate (2019). https://www.ipcc.ch/report/srocc/
2.
Huss, M., Lindmark, M., Jacobson, P., van Dorst, R. M. & Gårdmark, A. Experimental evidence of gradual size-dependent shifts in body size and growth of fish in response to warming. Glob. Chang. Biol. 25, 2285–2295 (2019).
PubMed  PubMed Central  Google Scholar 

3.
Daufresne, M., Lengfellner, K. & Sommer, U. Global warming benefits the small in aquatic ecosystems. Proc. Natl. Acad. Sci. USA 106, 12788–12793 (2009).
ADS  CAS  PubMed  Article  PubMed Central  Google Scholar 

4.
Free, C. M. et al. Impacts of historical warming on marine fisheries production. Science 363, 979–983 (2019).
ADS  CAS  PubMed  Article  PubMed Central  Google Scholar 

5.
Smoliński, S. & Mirny, Z. Otolith biochronology as an indicator of marine fish responses to hydroclimatic conditions and ecosystem regime shifts. Ecol. Indic. 79, 286–294 (2017).
Article  Google Scholar 

6.
Reed, T. E. et al. Responding to environmental change: Plastic responses vary little in a synchronous breeder. Proc. R. Soc. B 273, 2713–2719 (2006).
PubMed  Article  PubMed Central  Google Scholar 

7.
Ricker, W. Computation and interpretation of biological statistics of fish populations. Bull. Fish. Res. Board Can. 191, 382 (1975).
Google Scholar 

8.
Nussey, D. H., Clutton-Brock, T. H., Albon, S. D., Pemberton, J. & Kruuk, L. E. B. Constraints on plastic responses to climate variation in red deer. Biol. Lett. 1, 457–460 (2005).
PubMed  PubMed Central  Article  Google Scholar 

9.
Morrongiello, J. R. & Thresher, R. A statistical framework to explore ontogenetic growth variation among individuals and populations: a marine fish example. Ecol. Monogr. 85, 93–115 (2015).
Article  Google Scholar 

10.
Nussey, D. H., Wilson, A. J. & Brommer, J. E. The evolutionary ecology of individual phenotypic plasticity in wild populations. J. Evol. Biol. 20, 831–844 (2007).
CAS  PubMed  Article  PubMed Central  Google Scholar 

11.
Paoli, A., Weladji, R. B., Holand, Ø & Kumpula, J. Early-life conditions determine the between-individual heterogeneity in plasticity of calving date in reindeer. J. Anim. Ecol. 89, 370–383 (2020).
PubMed  Article  PubMed Central  Google Scholar 

12.
Campana, S. E. & Thorrold, S. R. Otoliths, increments, and elements: Keys to a comprehensive understanding of fish populations?. Can. J. Fish. Aquat. Sci. 58, 30–38 (2001).
Article  Google Scholar 

13.
Campana, S. E. Accuracy, precision and quality control in age determination, including a review of the use and abuse of age validation methods. J. Fish Biol. 59, 197–242 (2001).
Article  Google Scholar 

14.
Black, B. A., Boehlert, G. W. & Yoklavich, M. M. Using tree-ring crossdating techniques to validate annual growth increments in long-lived fishes. Can. J. Fish. Aquat. Sci. 62, 2277–2284 (2005).
Article  Google Scholar 

15.
Morrongiello, J. R., Thresher, R. E. & Smith, D. C. Aquatic biochronologies and climate change. Nat. Clim. Chang. 2, 849–857 (2012).
ADS  Article  Google Scholar 

16.
Clutton-Brock, T. & Sheldon, B. C. Individuals and populations: The role of long-term, individual-based studies of animals in ecology and evolutionary biology. Trends Ecol. Evol. 25, 562–573 (2010).
PubMed  Article  PubMed Central  Google Scholar 

17.
Grønkjær, P. Otoliths as individual indicators: A reappraisal of the link between fish physiology and otolith characteristics. Mar. Freshw. Res. 67, 881–888 (2016).
Article  Google Scholar 

18.
Bonamour, S., Chevin, L. M., Charmantier, A. & Teplitsky, C. Phenotypic plasticity in response to climate change: The importance of cue variation. Philos. Trans. R. Soc. B Biol. Sci. 374, 20180178 (2019).
Article  Google Scholar 

19.
Guindre-Parker, S. et al. Individual variation in phenotypic plasticity of the stress axis. Biol. Lett. 15, 1–7 (2019).
Article  Google Scholar 

20.
van de Pol, M. Quantifying individual variation in reaction norms: How study design affects the accuracy, precision and power of random regression models. Methods Ecol. Evol. 3, 268–280 (2012).
Article  Google Scholar 

21.
van de Pol, M. & Wright, J. A simple method for distinguishing within- versus between-subject effects using mixed models. Anim. Behav. 77, 753–758 (2009).
Article  Google Scholar 

22.
Morrongiello, J. R., Sweetman, P. C. & Thresher, R. E. Fishing constrains phenotypic responses of marine fish to climate variability. J. Anim. Ecol. 88, 1645–1656 (2019).
PubMed  Article  PubMed Central  Google Scholar 

23.
Dingemanse, N. J., Kazem, A. J. N., Réale, D. & Wright, J. Behavioural reaction norms: Animal personality meets individual plasticity. Trends Ecol. Evol. 25, 81–89 (2010).
PubMed  Article  PubMed Central  Google Scholar 

24.
Hilborn, R., Quinn, T. P., Schindler, D. E. & Rogers, D. E. Biocomplexity and fisheries sustainability. Proc. Natl. Acad. Sci. 100, 6564–6568 (2003).
ADS  CAS  PubMed  Article  PubMed Central  Google Scholar 

25.
Schindler, D. E. et al. Population diversity and the portfolio effect in an exploited species. Nature 465, 609–612 (2010).
ADS  CAS  PubMed  Article  PubMed Central  Google Scholar 

26.
Brommer, J. E., Merila, J., Sheldon, B. C. & Gustafsson, L. Natural selection and genetic variation for reproductive reaction norms in a wild bird population. Evolution 59, 1362–1371 (2005).
PubMed  Article  PubMed Central  Google Scholar 

27.
Link, J. S., Bogstad, B., Sparholt, H. & Lilly, G. R. Trophic role of Atlantic cod in the ecosystem. Fish Fish. 10, 58–87 (2009).
Article  Google Scholar 

28.
Drinkwater, K. F. The response of Atlantic cod (Gadus morhua) to future climate change. ICES J. Mar. Sci. 62, 1327–1337 (2005).
Article  Google Scholar 

29.
Richardson, A. J. et al. Climate change and marine life. Biol. Lett. 8, 907–909 (2012).
PubMed  PubMed Central  Article  Google Scholar 

30.
Lorenzen, K. & Enberg, K. Density-dependent growth as a key mechanism in the regulation of fish populations: Evidence from among-population comparisons. Proc. R. Soc. B Biol. Sci. 269, 49–54 (2002).
Article  Google Scholar 

31.
Frater, P. N., Hrafnkelsson, B., Elvarsson, B. T. & Stefansson, G. Drivers of growth for Atlantic cod (Gadus morhua L.) in Icelandic waters—A Bayesian approach to determine spatiotemporal variation and its causes. J. Fish Biol. 95, 401–410 (2019).
PubMed  Article  PubMed Central  Google Scholar 

32.
Eikeset, A. M. et al. Roles of density-dependent growth and life history evolution in accounting for fisheries-induced trait changes. Proc. Natl. Acad. Sci. USA 113, 15030–15035 (2016).
CAS  PubMed  Article  PubMed Central  Google Scholar 

33.
Brander, K. M. The effect of temperature on growth of Atlantic cod. ICES J. Mar. Sci. 52, 1–10 (1995).
Article  Google Scholar 

34.
Sinclair, A. F., Swain, D. P. & Hanson, J. M. Disentangling the effects of size-selective mortality, density, and temperature on length-at-age. Can. J. Fish. Aquat. Sci. 59, 372–382 (2002).
Article  Google Scholar 

35.
Pálsson, ÓK. A review of the trophic interactions of cod stocks in the North Atlantic. ICES Mar. Sci. Symp. 198, 553–575 (1994).
Google Scholar 

36.
Pálsson, ÓK. & Bjrnsson, H. Long-term changes in trophic patterns of Iceland cod and linkages to main prey stock sizes. ICES J. Mar. Sci. 68, 1488–1499 (2011).
Article  Google Scholar 

37.
Denechaud, C., Smoliński, S., Geffen, A. J., Godiksen, J. A. & Campana, S. E. A century of fish growth in relation to climate change, population dynamics and exploitation. Glob. Chang. Biol. 26, 5661–5678 (2020).
Article  Google Scholar 

38.
Beverton, R. J. H. & Holt, S. J. On the Dynamics of Exploited Fish Populations (Fisheries Investigations, 1957).

39.
Stige, L. C. et al. Density- and size-dependent mortality in fish early life stages. Fish Fish. 20, 962–976 (2019).
Article  Google Scholar 

40.
Linehan, J. E., Gregory, R. S. & Schneider, D. C. Predation risk of age-0 cod (Gadus) relative to depth and substrate in coastal waters. J. Exp. Mar. Biol. Ecol. 263, 25–44 (2001).
Article  Google Scholar 

41.
Mattson, S. Food and feeding habits of fish species over a soft sublittoral bottom in the Northeast Atlantic: 1. Cod (Gadus morhua L.) (Gadidae). Sarsia 75, 247–260 (1990).
Article  Google Scholar 

42.
Bromley, P. J. Evidence for density-dependent growth in North Sea gadoids. J. Fish Biol. 35, 117–123 (1989).
Article  Google Scholar 

43.
Schopka, S. A. Fluctuations in the cod stock off Iceland during the twentieth century in the fisheries and environment. ICES Mar. Sci. Symp. 198, 175–193 (1994).
Google Scholar 

44.
Martino, J. C., Fwoler, A. J., Doubleday, Z. A., Grammer, G. L. & Gillanders, B. M. Using otolith chronologies to understand long-term trends and extrinsic drivers of growth in fisheries. Ecosphere 10, e02553 (2019).
Article  Google Scholar 

45.
Stephens, P. A. & Sutherland, W. J. Consequences of the Allee effect for behaviour, ecology and conservation. Trends Ecol. Evol. 14, 401–405 (1999).
CAS  PubMed  Article  PubMed Central  Google Scholar 

46.
Brander, K. M. Patterns of distribution, spawning, and growth in North Atlantic cod: The utility of inter-regional comparisons. ICES Mar. Sci. Symp. 198, 406–413 (1994).
Google Scholar 

47.
Pálsson, ÓK. & Thorsteinsson, V. Migration patterns, ambient temperature, and growth of Icelandic cod (Gadus morhua): Evidence from storage tag data. Can. J. Fish. Aquat. Sci. 60, 1409–1423 (2003).
Article  Google Scholar 

48.
Tanner, S. E. et al. Regional climate, primary productivity and fish biomass drive growth variation and population resilience in a small pelagic fish. Ecol. Indic. 103, 530–541 (2019).
Article  Google Scholar 

49.
Zhai, L. et al. Phytoplankton phenology and production around Iceland and Faroes. Cont. Shelf Res. 37, 15–25 (2012).
ADS  Article  Google Scholar 

50.
Heath, M. R. et al. Winter distribution, ontogenetic migration, and rates of egg production of Calanus finmarchicus southwest of Iceland. ICES J. Mar. Sci. 57, 1727–1739 (2000).
Article  Google Scholar 

51.
Björnsson, B., Steinarsson, A. & Árnason, T. Growth model for Atlantic cod (Gadus morhua): Effects of temperature and body weight on growth rate. Aquaculture 271, 216–226 (2007).
Article  Google Scholar 

52.
Björnsson, B. & Steinarsson, A. The food-unlimited growth rate of Atlantic cod (Gadus morhua). Can. J. Fish. Aquat. Sci. 59, 494–502 (2002).
Article  Google Scholar 

53.
Arnold, P. A., Nicotra, A. B. & Kruuk, L. E. B. Sparse evidence for selection on phenotypic plasticity in response to temperature. Philos. Trans. R. Soc. B Biol. Sci. 374, 20180185 (2019).
Article  Google Scholar 

54.
Imsland, A. K. & Jónsdóttir, ÓD. B. Linking population genetics and growth properties of Atlantic cod. Rev. Fish Biol. Fish. 13, 1–26 (2003).
Article  Google Scholar 

55.
Imsland, A. K. et al. A retrospective approach to fractionize variation in body mass of Atlantic cod Gadus morhua. J. Fish Biol. 78, 251–264 (2011).
CAS  PubMed  Article  PubMed Central  Google Scholar 

56.
Nussey, D. H., Clutton-Brock, T. H., Elston, D. A., Albon, S. D. & Kruuk, L. E. B. Phenotypic plasticity in a maternal trait in red deer. J. Anim. Ecol. 74, 387–396 (2005).
Article  Google Scholar 

57.
Enberg, K. et al. Fishing-induced evolution of growth: Concepts, mechanisms and the empirical evidence. Mar. Ecol. 33, 1–25 (2012).
ADS  Article  Google Scholar 

58.
Rayner, N. A. et al. Global analyses of sea surface temperature, sea ice, and night marine air temperature since the late nineteenth century. J. Geophys. Res. 108, 27–35 (2003).
Article  Google Scholar 

59.
Sólmundsson, J., Jónsdóttir, I. G., Ragnarsson, S. A. & Björnsson, B. Connectivity among offshore feeding areas and nearshore spawning grounds: Implications for management of migratory fish. ICES J. Mar. Sci. 75, 148–157 (2018).
Article  Google Scholar 

60.
Weisberg, S., Spangler, G. & Richmond, L. S. Mixed effects models for fish growth. Can. J. Fish. Aquat. Sci. 277, 269–277 (2010).
Article  Google Scholar 

61.
Smoliński, S. Sclerochronological approach for the identification of herring growth drivers in the Baltic Sea. Ecol. Indic. 101, 420–431 (2019).
Article  Google Scholar 

62.
van de Pol, M. et al. Identifying the best climatic predictors in ecology and evolution. Methods Ecol. Evol. 7, 1246–1257 (2016).
Article  Google Scholar 

63.
Bailey, L. D. & van de Pol, M. climwin: An R toolbox for climate window analysis. PLoS ONE 11, 1–27 (2016).
Google Scholar 

64.
Dingemanse, N. J. & Dochtermann, N. A. Quantifying individual variation in behaviour: Mixed-effect modelling approaches. J. Anim. Ecol. 82, 39–54 (2013).
PubMed  Article  PubMed Central  Google Scholar 

65.
R Core Team. R: A Language and Environment for Statistical Computing (2018).

66.
Bates, D. et al.lme4: Linear Mixed-Effects Models using “Eigen” and S4. R package version 1.1-12. https://cran.r-project.org/web/packages/lme4/lme4.pdf (2016). doi:https://doi.org/10.18637/jss.v067.i01.

67.
Wessel, P. & Smith, W. H. F. A global, self-consistent, hierarchical, high-resolution shoreline database. J. Geophys. Res. Solid Earth 101, 8741–8743 (2004).
Article  Google Scholar 

68.
Amante, C. & Eakins, B. ETOPO1 1 Arc-Minute Global Relief Model: Procedures, Data Sources and Analysis. NOAA Technical Memorandum NESDIS NGDC-24 (2009). https://doi.org/10.7289/V5C8276M More

Fractional bulk density concept
The first assumption is that soil particles with different sizes contribute to different porosities and water holding capacities in bulk soil. Based on a non-similar media concept (NSMC) defined by Miyazaki49, soil bulk density (ρb) is defined as

$$rho _{{{b}}} = frac{{{M}}}{{{V}}} = tau rho _{{{s}}} left( {frac{{{S}}}{{{{S}} + {{d}}}}} right)^{3}$$
(4)

where M is the mass of a given soil, V is the volume of bulk soil, ρs is soil particle density, and S and d are characteristic lengths of solid phase and pore space, respectively. The parameter τ is a shape factor of the solid phase, defined as the ratio of the substantial volume of solid phase to the volume S3. The value of τ is 1.0 for a cube and π/6 for a sphere. As pointed out by Miyazaki49, these characteristic lengths are not directly measurable but are representative lengths in the sense of the characteristic length in a similar media concept (SMC). Following the approach of NSMC represented by Eq. (4), we conceptually defined the volume of bulk soil as

$$V = frac{{mathop sum nolimits_{{{i } = { 1}}}^{{{n}}} {{m}}_{{{i}}} }}{{{rho }_{{{b}}} }} = frac{{{{m}}_{{1}} }}{{{rho }_{{{{b1}}}} }}{ + }frac{{{{m}}_{{2}} }}{{{rho }_{{{{b2}}}} }}{ + } cdot{mkern -4mu}cdot{mkern -4mu}cdot frac{{{{m}}_{{{n}}} }}{{{rho }_{{{{bn}}}} }}$$
(5)

where mi and ρbi are the solid mass and equivalent bulk density of the ith size fraction of soil particles, respectively. In this study, diameters of the first particle fraction and the last one were assumed to be 1 µm and 1000 µm, respectively8. This equation suggests that different particle size fractions are associated with different equivalent bulk densities due to different contributions of particle arrangement to soil pore space. As a result, the particles with the same size fraction could have different equivalent bulk densities in soils with different textures or after the soil particles are rearranged (e.g., compaction). Figure 4 provides a diagrammatic representation of such fractional bulk density concept for the variation of soil pore volume with soil particle assemblage.
Figure 4

Diagrammatic representation of the fractional bulk density (FBD) model. V and ρb are the volume of bulk soil and the bulk density of whole soil, respectively. mi, and ρbi refer to the solid mass and equivalent bulk density associated with the ith particle-size fractions, respectively.

Full size image

Calculation of volumetric water content
For a specific soil, Eq. (5) means

$${{V}}_{{{{pi}}}} left( { le {{D}}_{{{i}}} } right){{ } = { f}}left( {{{D}}_{{{{gi}}}} {{, M}}_{{{i}}} } right)$$
(6)

where Vpi(≤ Di) denotes the volume of the pores with diameter ≤ Di generated by soil particles with diametes ≤ Dgi in unit volume of soil. Mi is the cumulative mass percentage of the ≤ Dgi particles. Since the pore volume has the maximum value for a given bulk soil and the cumulative distribution of pore volume could be generally hypothesized as a sigmoid curve for most of the natural soils44,45, we formulated Eq. (6) using a lognormal Logistic equation,

$${{V}}_{{{{pi}}}} left( { le {{D}}_{{{{gi}}}} } right) = frac{{{{V}}_{{{{pmax}}}} }}{{1 + kappa left( {{{D}}_{{{{gi}}}} } right)^{{{{b}}_{{{i}}} }} }}$$
(7)

where Vpmax is the maximum cumulative volume of pores pertinent to the particles smaller than or equal to the maximum diameter (Dgmax) in unit volume of soil. In fact, here Vpmax is equal to the total porosity (φT) of soil. Vpi (≤ Dgi) is the volume of the pores produced by ≤ Dgi particles in unit volume of soil, and bi is a varying parameter of increase in cumulative pore volume with an increment of Dgi. By assuming a complete saturation of soil pore space, Eq. (7) changes into

$$theta_{{{i}}} left( { le {{D}}_{{{{gi}}}} } right) = frac{{theta_{{{s}}} }}{{1 + kappa left( {{{D}}_{{{{gi}}}} } right)^{{{{b}}_{{{i}}} }} }}$$
(8)

where θs is saturated volumetric water content calculated with

$$theta_{{{s}}} = left{ {begin{array}{*{20}l} {0.9varphi _{{{T}}} ,} hfill & {~rho _{{{b}}} < 1} hfill \ {varphi _{{{T}}} ,} hfill & {~~rho _{{{b}}} ge 1} hfill \ end{array} } right.$$ (9) $$varphi _{{{T}}} = frac{{rho _{{{s}}} - rho _{{{b}}} }}{{rho _{{{s}}} }}$$ (10) In the above equations, ρbis measured soil bulk density, and ρs is soil particle density (2.65 g/cm3). The empirical parameter κ in Eqs. (7) and (8) is defined as $${{kappa}} = frac{{theta_{{{s}}} - theta_{{{r}}} }}{{theta_{{{r}}} }}$$ (11) where θr is measured residual water content. In this study, θr is set as the volumetric water content at water pressure head of 15,000 cm. The empirical parameter bi is defined as $${{b}}_{{{i}}} = frac{{epsilon }}{{3}}{log}left( {frac{{{theta }_{{{s}}} {{ - omega }}_{{{i}}} {theta }_{{{s}}} }}{{{{kappa omega }}_{{{i}}} {theta }_{{{s}}} }}} right)$$ (12) with ε, a particle size distribution index, calculated with $${varepsilon }; = ;frac{{left( {{{D}}_{{{40}}} } right)^{{2}} }}{{{{D}}_{{{10}}} {{D}}_{{{60}}} }}$$ (13) where D10, D40, and D60 represent the particle diameters below which the cumulative mass percentages of soil particles are 10%, 40%, and 60%, respectively. The parameter ωi is coefficient for soil particles of the ith size fraction, with a range of value between θr/θs and 1.0. By incorporating soil physical properties, ωi can be estimated with $${omega }_{{{i}}} = frac{{{g}}}{{{{1 + kappa }}left( {{{lnD}}_{{{{gi}}}} } right)^{{lambda}} }}$$ (14) where g is regulation coefficient (1.0–1.2). We set it to be 1.2 in this study. λ is the ratio coefficient of particle size distribution fitted using the lognormal Logistic model, $$M_{i} = frac{{M_{T} }}{{1 + eta D_{{gi}} ^{lambda } }}$$ (15) where MT represent the total mass percentage of all sizes of soil particles, and η is a fitting parameter. We set MT = 101 in Eq. (15) for best fit of the particle size distribution. In this study, this continuous function was generated from the discrete data pairs of Dgi and Mi at cutting particle diameters of 1,000, 750, 500, 400, 350, 300, 250, 200, 150, 100, 50, 30, 15, 7.5, 5, 3, 2, and 1 μm. Considering the difference in the upper limits of particle sizes associated with existing datasets of Dgi and Mi, the particle size distribution with the upper limit of 2,000 μm for the Acolian sandy soil and volcanic ash soils in Table 2 was normalized to the case with the upper limit of 1,000 μm using Eq. (3). Table 2 Physical properties of soils used in the study. ρb is bulk density (g/cm3); θr is residual water content (cm3/cm3) at 15,000 cm water pressure head; ε is particle size distribution index. Full size table Calculation of water pressure head To estimate the capillary tube or pore diameter (Di in µm), which was composed of particles with the size of Dgi (µm), Arya and Paris19 developed an expression $${{D}}_{{{i}}} {{ } = { D}}_{{{{gi}}}} left[ {frac{{2}}{{3}}{{en}}_{{{i}}}^{{{{(1 - alpha )}}}} } right]^{{{0}{{.5}}}}$$ (16) where α is the empirical scaling parameter varying between 1.35 and 1.40 in their original model19, but was thought to vary with soil particle size in the optimized model of Arya et al.20. In Tyler and Wheatcraft's model22α is the fractal dimension of the pore. The parameter e is the void rate of entire soil and assumed unchanging with particle size. However, according to Eqs. (5) and (6), e in Eq. (16) should vary with particle size and be replaced by ei, which depends on soil particle sizes. ni is the number of particles in the ith size fraction with a particle diameter (Dgi in μm), assuming that the particles are spherical and that the entire pore volume formed by assemblage of the particles in this class is represented by a single cylindrical pore. The equation for calculating ni is given as19 $$n_{i} = frac{{6m_{i} }}{{rho_{s} pi D_{gi}^{3} }} times 10^{12}$$ (17) where mi is the mass of particles in the ith size fraction of particles. Assuming that soil water has a zero contact angle and a surface tension of 0.075 N/m at 25 °C, the minimum diameter of soil pore (Dmin) was taken to be 0.2 µm in this study, which is equivalent to the water pressure head of 15,000 cm according to Young–Laplace equation. We set this minimum pore size to correspond the minimum particle size (Dgmin = 1.0 µm). The FBD model might thus not apply well to porous media with pores smaller than 0.2 μm. As a result, Eq. (16) can be simplified into the following equation. $${{D}}_{{{i}}} { = 0}{{.2D}}_{{{{gi}}}}$$ (18) The equivalent capillary pressure (ψi in cm) corresponding to the ith particle size fraction can be calculated using $$psi_{{{i}}} = frac{{{3000}}}{{{{D}}_{{{i}}} }} = frac{{{15000}}}{{{{D}}_{{{{gi}}}} }}$$ (19) In Eq. (19), the maximum water pressure head (ψr = 15,000 cm) corresponds to θr and Dgmin (1 μm). The minimum water pressure head (ψ0 = 15 cm) corresponds to θs and Dgmax (1,000 μm). These assumptions were arbitrary and might not be appropriate for some soil types. But these values were used in the study because they approximated the practical range of measurements well. The resulting model of soil water retention Equations 8 and 19 formulate a FBD-based model for estimation of soil water retention curve. To simplify the computation, we incorporated the two equations into the following analytical form, $${theta }; = ;frac{{{theta }_{{{s}}} }}{{{1 + }left( {frac{{{theta }_{{{s}}} - {theta }_{{{r}}} }}{{{theta }_{{{r}}} }}} right)left( {frac{{15,000}}{{psi }}} right)^{{{b}}} }}$$ (20) with the parameter b obtained using $${{b}}; = ;frac{{epsilon }}{{3}}{log}left{ {frac{{{{(theta }}_{{{s}}} - {theta }_{{{r}}} {{)[ln(}}frac{{{15,000}{{.1}}}}{{psi }}{)]}^{{lambda }} - {{(g}} - {{1)theta }}_{{{r}}} }}{{{{g(theta }}_{{{s}}} - {theta }_{{{r}}} {)}}}} right}$$ (21) In Eq. (21), a water pressure head of 15,000.1 cm is employed to consecutively predict the soil water content until the water pressure head of 15,000 cm. Soil dataset Evaluation of the applicability of the proposed modeling procedure required datasets that included soil bulk density, residual water content, and soil particle size distribution covering three particle diameters (D10, D40, and D60) below which the cumulative mass fractions of particles were 10%, 40%, and 60%, respectively. In addition, measured water content and water pressure head were required for the actual retention curve in order to compare with the result of the FBD model. In this study, the soil water retention data of 30 different soils, measured by Yu et al.50, Chen and Wang51, Zhang and Miao52, Liu and Amemiya53, Hayano et al.54, and Yabashi et al.55 were used for model verification (Table 2). The data covered soils in China (such as black soil, chernozem soil, cinnamon soil, brown earth, fluvo-aquic soil, albic soil, red earth, humid-thermo ferralitic, purplish soil, meadow soil, and yellow earth) and soils in Japan (such as volcanic ash soil and acolian sandy soil). The USDA soil taxonomy of these soils was provided in Table 2. The 30 soils ranged in texture from clay to sand and in bulk density from 0.33 g/cm3 to 1.65 g/cm3, which covered a much wider range of soil bulk density than many of the existing models or pedotransfer functions56,57,58,59. Particle size fractions (Dgi) were chosen as the upper limit of the diameters between successive sieve sizes. For the data set in which particle density was not determined, 2.65 g/cm3 was used. More

Lethal toxicity
Living animals are constantly faced with various environmental stresses that challenge their daily lives. Cu is an essential metal that participates in normal physiological process of crustaceans, but several studies have shown that crustaceans are adversely affected when exposed to high concentrations of Cu. LC50 value represents a common point at lethal physiological response to toxicity, which has been well-documented in many crustaceans. For example, the 96-h LC50 value for shrimps of Exopalaemon carinicauda, Echinogammars olivii, Sphaeroma serratum, and Palaemon elegans was 0.712 mg Cu/L, 0.25 mg Cu/L, 1.98 mg Cu/L, and 2.52 mg Cu/L, respectively20,22. In addition, for paddy field crab Paratelphusa hydrodromus and freshwater crab, Barytelphusa cunicularis, the 96-h LC50 values recorded were 15.70 mg Cu/L and 215 mg Cu/L, respectively23,24. Likewise, in freshwater crayfish, Procambarus clarkia, the 96-h LC50 value reached 162 mg Cu/L25. These large variations in sub-lethal effects to Cu toxicity in crustaceans appear to be species specific. In our present study, the 96-h LC50 value for Cu exposure in C. dauricus was 32.5 mg/L, which is much higher than those of the most crustaceans, but this species seems relatively less tolerant to Cu, compared to P. clarkia. This difference may also be attributed to the various biotic and abiotic factors like age, sex, weight, salinity, and temperature, besides the species. For example, Taylor26 compared the 96-h Cu tolerance of Cambarus robustus in Pike Creek and in Wavy Lake and concluded the environment differences could affect population sensitivity to Cu toxicity.
Oxygen consumption rate
The effects of heavy metal on the respiratory rate of marine and estuarine organisms have been well documented. Spicer and Weber27 showed that heavy metal could cause respiratory impairment in crustaceans. The results obtained in the present study confirmed this previous finding. Both acute and sub-chronic Cu exposure induced significant inhibition of OCR in C. dauricus, with the maximum decreases of 48.4% and 57.9%, respectively, compared to the control. Similarly, a declined OCR by heavy metal has been observed in shrimps, including Penaeus indicus19, L. vannamei28, F. paulensis29, and E. carinicauda20, and crabs, including Uca annulipes, U. triangularis30, and Cancer pagarus27, as well as crayfish, P. clarkia25. The levels of inhibition of the respiration rate were mainly dependent on the exposure time and exposure concentration. Those authors assumed that the ultrastructural impairments of gill epithelium were related to the decrease in respiration rate, thereby affecting the oxygen carrying capacity of gills. Besides the cytological damage of gill, heavy metals also inhibit mitochondrial energy production, thereby affecting the key metabolic pathways. By contrast, an increased respiration rate has been detected in freshwater shrimp, Paratya curvirostris21, and lobster Homarus americanus31. The authors argued that it was attributed to an elevated rate of glycolysis, a mechanism of expenditure of energy reserves characteristic of a stress compensation process. In all, the changes of oxygen consumption level were mainly dependent on the time and concentration of exposure to heavy metals.
Ammonia excretion rate
Amino acids are the main sources of ammonia production in vivo. Crustaceans have the ability to regulate the concentration of intracellular free amino acids in order to deal with environmental stress32. In the present study, AER in either acute or sub-chronic Cu exposure showed a declining trend with increasing exposure concentration and time to Cu. A maximum decrease in AER of 79.4% and 70.06%, respectively, were observed respectively after exposure to 16.48 mg/L for 96-h and 2.06 mg/L for 14 days, in comparison to the control (Fig. 1B). In a similar manner, Chinni19 also reported a significant decrease in AER in post larvae P. indicus when exposed to Pb for 30 days. It assumed that such a decrease may be due to reduction in the metabolic rate or an interaction of heavy metal with the pathways for the production of ammonia-N. By contrast, elevations of ammonia excretion in response to heavy metals exposure were reported in other crustaceans. For example, an increase in AER was found in juvenile E. carinicauda after exposure to Zn and Hg20 and in F. paulensis after exposure to Cd and Zn29. It was considered that the gill function was impaired by the metal exposure, resulting in the dysfunction of ammonium excretion control; therefore, outflow of ammonia excretion from the hemolymph to ambient water induced an increased ammonia concentration in the water. In addition, no change in ammonia excretion rate was obtained in Paratya curvirostris after 96-h acute and 10-day sub-chronic Cd stress21. Therefore, the questions of the relationship between heavy metal exposure and ammonia excretion needs to be properly investigated.
Energy metabolism
O:N is a useful value for evaluating the characteristics of nutrients utilized by animals and can provide information on changes in energy substrate utilization under various environmental stresses33,34. Theoretically, pure protein catabolism will produce an O:N ratio of 835, and equal proportions of proteins and lipid results in an O:N of 2436. An O:N ratio higher than 24 indicates an elevation in lipid and carbohydrate metabolisms. In this study, in comparison with the controls, high values of O:N were obtained in individuals of C. dauricus exposed to Cu for 96 h and 14 days (Table 1, Fig. 1C). In generally, protein catabolism for energy is less efficient than lipid/carbohydrate catabolism. A species that relies on lipid and carbohydrate metabolism will likely be able to better meet energy demands of toxicant exposure than a species that principally metabolizes protein. The mean O:N ratio higher than 24 in acute Cu exposure and lower than 24 in sub-chronic exposure (Table 3 and Fig. 1C) indicated the differences in energy utilization strategy in response to two patterns of Cu stress. This could be a mechanism explaining the differences in energetic responses to Cu exposure in C. dauricus, relative to other crustacean species.
Tissues accumulation
Cu is an essential trace element for biological processes, particularly as a component of the respiratory pigment, hemocyanin. The body Cu concentration in decapod crustaceans can be regulated and does not accumulate until certain environmental threshold levels are achieved37. In addition, as an economic species of crustaceans and in relation to food quality and safety assessment, organ-specific accumulation data, especially for the muscle, are markedly required. In this study, tissue-specific bioaccumulation of Cu observed, and the Cu accumulation in hepatopancreas and muscles were highly dependent on water Cu concentration and exposure time (Fig. 2; Fig. 3A, 3B). Hepatopancreas is the organ most associated with the detoxification and biotransformation process and in direct contact with toxicants in water. The hepatopancreas, containing metal-binding protein, is the main target organ for regulating Cu level38. The maximum Cu accumulation was observed in hepatopancreas, which increased 12.7 folds and 31.6 folds after 4-day acute exposure to 16.48 mg Cu/L and chronic 14-day exposure to 4.12 mg Cu/L, respectively, this indicated that C. dauricus had a great potential for rapid accumulation of Cu in fresh waters. The greatest Cu accumulation occurring in hepatopancreas had been reported for the crayfish species, Astacus leptodactylus39 and Procambarus sp.40 as well as for the freshwater prawn, M. rosenbergii38. Although the hepatopancreas could regulate the Cu level in the animal’s body to avoid toxicity and deficiencies, the high level of external water Cu breaks down the regulation of Cu and causes continuous Cu accumulation, which might lead to the loss of muscular control and eventually, death, for crustaceans.
In this study, there was no significant time-dependent trend in the accumulation of Cu in the muscle between 7 and 14 days of Cu stress in the lower concentration of 2.06 mg/L (Fig. 3A), this suggests that C. dauricus was able to regulate Cu in the muscle to a fairly constant level under low Cu exposure concentrations. However, C. dauricus exposed to concentration of 4.12 mg Cu/L showed increased accumulation of Cu in the muscle and the equilibrium of Cu accumulation was not reached at 14 days, which might show that the high level of Cu in the external water breaks down the regulation of Cu and caused a continuous Cu accumulation, leading to its toxicity at high concentration. Similar result had been reported in Procambarus sp.40. The author found that Cu uptake reached a kinetic equilibrium within 10 days of exposure to 0.31 mg Cu/L in five organs (gills, ovaries, exoskeleton, hepatopancreas, and muscles), but Cu was rapidly accumulated in the organs of most Procambarus sp., especially in the hepatopancreas, when exposed to higher concentration of 0.38 mg Cu/L after the 14-d exposure test. However, muscle tissue, as the main edible portion, accumulates Cu at a relatively lower rate (Fig. 2; Fig. 3A) and this is important from the angle of human food quality and safety.
Conclusion
In this study, we observed that the acute and sub-chronic toxicity of Cu had a dramatic impact on the survival, oxygen consumption rate, ammonia excretion rate and bioaccumulation of C. dauricus. C. dauricus mainly took the strategies of inhibiting respiratory metabolism and shifting energy utilization to adapt to copper stress. The C. dauricus had higher concentration-dependent accumulation ability of copper. Our future work will focus on the metabolic characteristics of copper and other heavy metal from the angle of human food safety. Therefore, our studies provided basic information for further understanding of the toxicological responses of this species to trace metals. More

1.
Watson, J. E. M., Dudley, N., Segan, D. B. & Hockings, M. The performance and potential of protected areas. Nature 515, 67–73 (2014).
ADS  CAS  Google Scholar 
2.
Dudley, N. Guidelines for Applying Protected Area Management Categories (IUCN, 2008).

3.
Dudley, N. et al. The essential role of other effective area-based conservation measures in achieving big bold conservation targets. Glob. Ecol. Conserv. 15, e00424 (2018).
Google Scholar 

4.
Donald, P. F. et al. The prevalence, characteristics and effectiveness of Aichi Target 11′ s “other effective area-based conservation measures”(OECMs) in Key Biodiversity Areas. Conserv. Lett. 12, 12659 (2019).
Google Scholar 

5.
UN General Assembly. Transforming our World: The 2030 Agenda for Sustainable Development, 21 October 2015. A/RES/70/1 https://www.refworld.org/docid/57b6e3e44.html (accessed 11 November 2019).

6.
Convention on Biological Diversity. COP 10 Decision X/2: Strategic Plan for Biodiversity 2011–2020. http://www.cbd.int/decision/cop/?id=12268 (2011).

7.
UNEP-WCMC & IUCN. World Database on Protected Areas (WDPA). https://www.protectedplanet.net/ (UNEP-WCMC, 2019).

8.
UNEP-WCMC & IUCN. World Database on Other Effective Area-based Conservation Measures (WD-OCEM). https://www.protectedplanet.net/c/other-effective-area-based-conservation-measures (UNEP-WCMC, 2019).

9.
Lewis, E. et al. Dynamics in the global protected-area estate since 2004. Conserv. Biol. 33, 570–579 (2019).
Google Scholar 

10.
Klein, C. J. et al. Shortfalls in the global protected area network at representing marine biodiversity. Sci. Rep. 5, 17539 (2015).
ADS  PubMed  PubMed Central  Google Scholar 

11.
Venter, O. et al. Bias in protected-area location and its effects on long-term aspirations of biodiversity conventions. Conserv. Biol. 32, 127–134 (2018).
Google Scholar 

12.
Mouillot, D. et al. Global marine protected areas do not secure the evolutionary history of tropical corals and fishes. Nat. Commun. 7, 10359 (2016).
ADS  CAS  PubMed  PubMed Central  Google Scholar 

13.
Butchart, S. H. M. et al. Shortfalls and solutions for meeting national and global conservation area targets. Conserv. Lett. 8, 329–337 (2015).
Google Scholar 

14.
Christie, P. et al. Why people matter in ocean governance: incorporating human dimensions into large-scale marine protected areas. Mar. Policy 84, 273–284 (2017).
Google Scholar 

15.
Zafra-Calvo, N. et al. Progress toward equitably managed protected areas in Aichi target 11: a global survey. Bioscience 69, 191–197 (2019). This is the first large review of how well protected areas satisfy social equity metrics.
PubMed  PubMed Central  Google Scholar 

16.
Juffe-Bignoli, D. et al. Achieving Aichi biodiversity target 11 to improve the performance of protected areas and conserve freshwater biodiversity. Aquat. Conserv. 26, 133–151 (2016).
Google Scholar 

17.
Maron, M., Simmonds, J. S. & Watson, J. E. M. Bold nature retention targets are essential for the global environment agenda. Nat. Ecol. Evol. 2, 1194–1195 (2018).
Google Scholar 

18.
Geldmann, J. et al. Changes in protected area management effectiveness over time: a global analysis. Biol. Conserv. 191, 692–699 (2015).
Google Scholar 

19.
Di Minin, E. & Toivonen, T. Global protected area expansion: creating more than paper parks. Bioscience 65, 637–638 (2015).
PubMed  PubMed Central  Google Scholar 

20.
Gill, D. A. et al. Capacity shortfalls hinder the performance of marine protected areas globally. Nature 543, 665–669 (2017). This study compiles four years of data to assess capacity shortfalls and biodiversity outcomes from the management of 589 marine protected areas.
ADS  CAS  Google Scholar 

21.
Coad, L. et al. Widespread shortfalls in protected area resourcing undermine efforts to conserve biodiversity. Front. Ecol. Environ. 17, 259–264 (2019).
Google Scholar 

22.
Visconti, P. et al. Protected area targets post-2020. Science 364, 239–241 (2019).
ADS  CAS  Google Scholar 

23.
Barnes, M. D., Glew, L., Wyborn, C. & Craigie, I. D. Prevent perverse outcomes from global protected area policy. Nat. Ecol. Evol. 2, 759–762 (2018).
Google Scholar 

24.
IPBES. Summary for Policymakers of the Global Assessment Report on Biodiversity and Ecosystem Services of the Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services (IPBES secretariat, 2019). This report assesses the status of biodiversity and ecosystem services, their impact on human well-being and the effectiveness of conservation interventions.

25.
Dinerstein, E. et al. A global deal for nature: guiding principles, milestones, and targets. Sci. Adv. 5, eaaw2869 (2019).
ADS  CAS  PubMed  PubMed Central  Google Scholar 

26.
Noss, R. F. et al. Bolder thinking for conservation. Conserv. Biol. 26, 1–4 (2012).
Google Scholar 

27.
Wilson, E. O. Half-Earth: Our Planet’s Fight for Life (Liveright, 2016).

28.
O’Leary, B. C. et al. Effective coverage targets for ocean protection. Conserv. Lett. 9, 398–404 (2016).
Google Scholar 

29.
Bull, J. W. et al. Net positive outcomes for nature. Nat. Ecol. Evol. 4, 4–7 (2020).
Google Scholar 

30.
Mace, G. M. et al. Aiming higher to bend the curve of biodiversity loss. Nat. Sustain. 1, 448–451 (2018).
Google Scholar 

31.
Dinerstein, E. et al. An ecoregion-based approach to protecting half the terrestrial realm. Bioscience 67, 534–545 (2017).
PubMed  PubMed Central  Google Scholar 

32.
Spalding, M. D. et al. Marine ecoregions of the world: a bioregionalization of coastal and shelf areas. Bioscience 57, 573–583 (2007).
Google Scholar 

33.
UNEP-WCMC, IUCN & NGS. Protected Planet Report 2018 (UNEP-WCMC, IUCN and NGS, 2018). A biennial publication that reviews progress toward protected areas targets and goals.

34.
Rodrigues, A. S. L. et al. Global gap analysis: priority regions for expanding the global protected-area network. Bioscience 54, 1092–1100 (2004).
Google Scholar 

35.
IUCN. The IUCN Red List of Threatened Species. Version 2019-2 http://www.iucnredlist.org (accessed 10 September 2019) (2019).

36.
IUCN. A Global Standard for the Identification of Key Biodiversity Areas. Version 1.0 (IUCN, 2016).

37.
BirdLife International. World Database of Key Biodiversity Areas. www.keybiodiversityareas.org (accessed 20 June 2019) (2019).

38.
Jones, K. R. et al. The location and protection status of Earth’s diminishing marine wilderness. Curr. Biol. 28, 2506–2512 (2018).
CAS  PubMed  PubMed Central  Google Scholar 

39.
Allan, J. R., Venter, O. & Watson, J. E. M. Temporally inter-comparable maps of terrestrial wilderness and the last of the wild. Sci. Data 4, 170187 (2017).
PubMed  PubMed Central  Google Scholar 

40.
Watson, J. E. M. et al. The exceptional value of intact forest ecosystems. Nat. Ecol. Evol. 2, 599–610 (2018).
Google Scholar 

41.
Di Marco, M., Ferrier, S., Harwood, T. D., Hoskins, A. J. & Watson, J. E. M. Wilderness areas halve the extinction risk of terrestrial biodiversity. Nature 573, 582–585 (2019).
ADS  Google Scholar 

42.
Martin, T. G. & Watson, J. E. M. Intact ecosystems provide best defence against climate change. Nat. Clim. Chang. 6, 122–124 (2016).
ADS  Google Scholar 

43.
Griscom, B. W. et al. Natural climate solutions. Proc. Natl Acad. Sci. USA 114, 11645–11650 (2017).
ADS  CAS  Google Scholar 

44.
Soto-Navarro, C. et al. Mapping co-benefits for carbon storage and biodiversity to inform conservation policy and action. Phil. Trans. R. Soc. Lond. B 375, 20190128 (2020). This study combines multiple datasets to produce a new high-resolution map of global above- and belowground carbon stored in biomass and soil.
CAS  Google Scholar 

45.
Dargie, G. C. et al. Age, extent and carbon storage of the central Congo Basin peatland complex. Nature 542, 86–90 (2017).
ADS  CAS  Google Scholar 

46.
DeVries, T. & Weber, T. The export and fate of organic matter in the ocean: new constraints from combining satellite and oceanographic tracer observations. Glob. Biogeochem. Cycles 31, 535–555 (2017).
ADS  CAS  Google Scholar 

47.
Laws, E. A., D’Sa, E. & Naik, P. Simple equations to estimate ratios of new or export production to total production from satellite-derived estimates of sea surface temperature and primary production. Limnol. Oceanogr. Methods 9, 593–601 (2011).
Google Scholar 

48.
DeVries, T., Primeau, F. & Deutsch, C. The sequestration efficiency of the biological pump. Geophys. Res. Lett. 39, L13601 (2012).
ADS  Google Scholar 

49.
Henson, S. A., Sanders, R. & Madsen, E. Global patterns in efficiency of particulate organic carbon export and transfer to the deep ocean. Glob. Biogeochem. Cycles 26, GB1028 (2012).
ADS  Google Scholar 

50.
Roshan, S. & DeVries, T. Efficient dissolved organic carbon production and export in the oligotrophic ocean. Nat. Commun. 8, 2036 (2017).
ADS  PubMed  PubMed Central  Google Scholar 

51.
Lutz, M. J., Caldeira, K., Dunbar, R. B. & Behrenfeld, M. J. Seasonal rhythms of net primary production and particulate organic carbon flux to depth describe the efficiency of biological pump in the global ocean. J. Geophys. Res. Oceans 112, C10011 (2007).
ADS  Google Scholar 

52.
Magris, R. A. et al. Biologically representative and well-connected marine reserves enhance biodiversity persistence in conservation planning. Conserv. Lett. 11, e12439 (2018).
Google Scholar 

53.
Mendenhall, C. D., Karp, D. S., Meyer, C. F. J., Hadly, E. A. & Daily, G. C. Predicting biodiversity change and averting collapse in agricultural landscapes. Nature 509, 213–217 (2014).
ADS  CAS  Google Scholar 

54.
Harrison, H. B. et al. Larval export from marine reserves and the recruitment benefit for fish and fisheries. Curr. Biol. 22, 1023–1028 (2012).
CAS  Google Scholar 

55.
Johnson, D. W., Christie, M. R., Pusack, T. J., Stallings, C. D. & Hixon, M. A. Integrating larval connectivity with local demography reveals regional dynamics of a marine metapopulation. Ecology 99, 1419–1429 (2018).
Google Scholar 

56.
Saura, S., Bastin, L., Battistella, L., Mandrici, A. & Dubois, G. Protected areas in the world’s ecoregions: how well connected are they? Ecol. Indic. 76, 144–158 (2017).
PubMed  PubMed Central  Google Scholar 

57.
Saura, S. et al. Global trends in protected area connectivity from 2010 to 2018. Biol. Conserv. 238, 108183 (2019).
PubMed  PubMed Central  Google Scholar 

58.
Endo, C. A. K., Gherardi, D. F. M., Pezzi, L. P. & Lima, L. N. Low connectivity compromises the conservation of reef fishes by marine protected areas in the tropical South Atlantic. Sci. Rep. 9, 8634 (2019).
ADS  PubMed  PubMed Central  Google Scholar 

59.
Bergseth, B. J., Gurney, G. G., Barnes, M. L., Arias, A. & Cinner, J. E. Addressing poaching in marine protected areas through voluntary surveillance and enforcement. Nat. Sustain. 1, 421–426 (2018). This study uses a citizen science approach to estimate poaching rates inside 55 marine protected areas spanning seven countries.
Google Scholar 

60.
Jones, K. R. et al. One-third of global protected land is under intense human pressure. Science 360, 788–791 (2018).
CAS  Google Scholar 

61.
Costello, M. J. & Ballantine, B. Biodiversity conservation should focus on no-take marine reserves: 94% of marine protected areas allow fishing. Trends Ecol. Evol. 30, 507–509 (2015).
Google Scholar 

62.
Zupan, M. et al. Marine partially protected areas: drivers of ecological effectiveness. Front. Ecol. Environ. 16, 381–387 (2018).
Google Scholar 

63.
Spracklen, B. D., Kalamandeen, M., Galbraith, D., Gloor, E. & Spracklen, D. V. A global analysis of deforestation in moist tropical forest protected areas. PLoS ONE 10, e0143886 (2015).
CAS  PubMed  PubMed Central  Google Scholar 

64.
Herrera, D., Pfaff, A. & Robalino, J. Impacts of protected areas vary with the level of government: comparing avoided deforestation across agencies in the Brazilian Amazon. Proc. Natl Acad. Sci. USA 116, 14916–14925 (2019).
CAS  Google Scholar 

65.
Negret, P. J. et al. Effects of spatial autocorrelation and sampling design on estimates of protected area effectiveness. Conserv. Biol. https://doi.org/10.1111/cobi.13522 (2020).

66.
White, T. D. et al. Assessing the effectiveness of a large marine protected area for reef shark conservation. Biol. Conserv. 207, 64–71 (2017).
Google Scholar 

67.
Giakoumi, S. & Pey, A. Assessing the effects of marine protected areas on biological invasions: a global review. Front. Mar. Sci. 4, 49 (2017).

68.
Geldmann, J., Manica, A., Burgess, N. D., Coad, L. & Balmford, A. A global-level assessment of the effectiveness of protected areas at resisting anthropogenic pressures. Proc. Natl Acad. Sci. USA, 116, 23209–23215 (2019).
ADS  CAS  Google Scholar 

69.
Gray, C. L. et al. Local biodiversity is higher inside than outside terrestrial protected areas worldwide. Nat. Commun. 7, 12306 (2016). This controlled study shows how biodiversity outcomes from protected area management are mediated by different classes of land use.
ADS  CAS  PubMed  PubMed Central  Google Scholar 

70.
Kerwath, S. E., Winker, H., Götz, A. & Attwood, C. G. Marine protected area improves yield without disadvantaging fishers. Nat. Commun. 4, 2347 (2013).
ADS  Google Scholar 

71.
Speed, C. W., Cappo, M. & Meekan, M. G. Evidence for rapid recovery of shark populations within a coral reef marine protected area. Biol. Conserv. 220, 308–319 (2018).
Google Scholar 

72.
Caselle, J. E., Rassweiler, A., Hamilton, S. L. & Warner, R. R. Recovery trajectories of kelp forest animals are rapid yet spatially variable across a network of temperate marine protected areas. Sci. Rep. 5, 14102 (2015).
ADS  CAS  PubMed  PubMed Central  Google Scholar 

73.
Emslie, M. J. et al. Expectations and outcomes of reserve network performance following re-zoning of the Great Barrier Reef marine park. Curr. Biol. 25, 983–992 (2015).
CAS  Google Scholar 

74.
Campbell, S. J., Edgar, G. J., Stuart-Smith, R. D., Soler, G. & Bates, A. E. Fishing-gear restrictions and biomass gains for coral reef fishes in marine protected areas. Conserv. Biol. 32, 401–410 (2018).
Google Scholar 

75.
Mumby, P. J. et al. Trophic cascade facilitates coral recruitment in a marine reserve. Proc. Natl Acad. Sci. USA 104, 8362–8367 (2007).
ADS  CAS  Google Scholar 

76.
Boaden, A. E. & Kingsford, M. J. Predators drive community structure in coral reef fish assemblages. Ecosphere 6, art46 (2015).
Google Scholar 

77.
Lamb, J. B., Williamson, D. H., Russ, G. R. & Willis, B. L. Protected areas mitigate diseases of reef-building corals by reducing damage from fishing. Ecology 96, 2555–2567 (2015).
Google Scholar 

78.
Naidoo, R. et al. Evaluating the impacts of protected areas on human well-being across the developing world. Sci. Adv. 5, eaav3006 (2019).
ADS  CAS  PubMed  PubMed Central  Google Scholar 

79.
Zafra-Calvo, N. et al. Towards an indicator system to assess equitable management in protected areas. Biol. Conserv. 211, 134–141 (2017).
Google Scholar 

80.
Oldekop, J. A., Holmes, G., Harris, W. E. & Evans, K. L. A global assessment of the social and conservation outcomes of protected areas. Conserv. Biol. 30, 133–141 (2016).
CAS  Google Scholar 

81.
Giakoumi, S. et al. Revisiting “success” and “failure” of marine protected areas: a conservation scientist perspective. Front. Mar. Sci. 5, 223 (2018).

82.
Edgar, G. J. et al. Global conservation outcomes depend on marine protected areas with five key features. Nature 506, 216–220 (2014).
ADS  CAS  Google Scholar 

83.
Ban, N. C. et al. Well-being outcomes of marine protected areas. Nat. Sustain. 2, 524–532 (2019).
Google Scholar 

84.
Corrigan, C. et al. Quantifying the contribution to biodiversity conservation of protected areas governed by indigenous peoples and local communities. Biol. Conserv. 227, 403–412 (2018).
Google Scholar 

85.
Schleicher, J., Peres, C. A., Amano, T., Llactayo, W. & Leader-Williams, N. Conservation performance of different conservation governance regimes in the Peruvian Amazon. Sci. Rep. 7, 11318 (2017).
ADS  PubMed  PubMed Central  Google Scholar 

86.
Hoffmann, M. et al. The difference conservation makes to extinction risk of the world’s ungulates. Conserv. Biol. 29, 1303–1313 (2015).
Google Scholar 

87.
Watson, J. E. M. et al. Set a global target for ecosystems. Nature 578, 360–362 (2020).
ADS  CAS  Google Scholar 

88.
Stolton, S., Redford, K. H. & Dudley, N. The Futures of Privately Protected Areas (IUCN, 2014).

89.
IUCN WCPA. Guidelines for Recognising and Reporting Other Effective Area-based Conservation Measures (IUCN, 2019).

90.
Shabtay, A., Portman, M. E., Manea, E. & Gissi, E. Promoting ancillary conservation through marine spatial planning. Sci. Total Environ. 651, 1753–1763 (2019).
ADS  CAS  Google Scholar 

91.
Banks-Leite, C. et al. Using ecological thresholds to evaluate the costs and benefits of set-asides in a biodiversity hotspot. Science 345, 1041–1045 (2014).
ADS  CAS  Google Scholar 

92.
Schuster, R., Germain, R. R., Bennett, J. R., Reo, N. J. & Arcese, P. Vertebrate biodiversity on indigenous-managed lands in Australia, Brazil, and Canada equals that in protected areas. Environ. Sci. Policy 101, 1–6 (2019).
Google Scholar 

93.
Bennett, N. J. & Dearden, P. From measuring outcomes to providing inputs: governance, management, and local development for more effective marine protected areas. Mar. Policy 50, 96–110 (2014).
Google Scholar 

94.
Suchley, A. & Alvarez-Filip, L. Local human activities limit marine protection efficacy on Caribbean coral reefs. Conserv. Lett. 11, e12571 (2018).
Google Scholar 

95.
Cook, C. N., Valkan, R. S., Mascia, M. B. & McGeoch, M. A. Quantifying the extent of protected-area downgrading, downsizing, and degazettement in Australia. Conserv. Biol. 31, 1039–1052 (2017).
Google Scholar 

96.
Qin, S. et al. Protected area downgrading, downsizing, and degazettement as a threat to iconic protected areas. Conserv. Biol. 33, 1275–1285 (2019).
PubMed  PubMed Central  Google Scholar 

97.
Forrest, J. L. et al. Tropical deforestation and carbon emissions from protected area downgrading, downsizing, and degazettement (PADDD). Conserv. Lett. 8, 153–161 (2015).
Google Scholar 

98.
Golden Kroner, R. E. et al. The uncertain future of protected lands and waters. Science 364, 881–886 (2019). This study compiled data that are available globally on PADDD events.
ADS  CAS  Google Scholar 

99.
Roberts, K. E., Valkan, R. S. & Cook, C. N. Measuring progress in marine protection: a new set of metrics to evaluate the strength of marine protected area networks. Biol. Conserv. 219, 20–27 (2018).
Google Scholar 

100.
De Vos, A., Clements, H. S., Biggs, D. & Cumming, G. S. The dynamics of proclaimed privately protected areas in South Africa over 83 years. Conserv. Lett. 12, e12644 (2019).
Google Scholar 

101.
Costelloe, B. et al. Global biodiversity indicators reflect the modeled impacts of protected area policy change. Conserv. Lett. 9, 14–20 (2016).
Google Scholar 

102.
Pringle, R. M. Upgrading protected areas to conserve wild biodiversity. Nature 546, 91–99 (2017).
ADS  CAS  Google Scholar 

103.
Kuempel, C. D., Adams, V. M., Possingham, H. P. & Bode, M. Bigger or better: the relative benefits of protected area network expansion and enforcement for the conservation of an exploited species. Conserv. Lett. 11, e12433 (2018).
Google Scholar 

104.
Adams, V. M., Barnes, M. & Pressey, R. L. Shortfalls in conservation evidence: moving from ecological effects of interventions to policy evaluation. One Earth 1, 62–75 (2019).
Google Scholar 

105.
Coad, L. et al. Measuring impact of protected area management interventions: current and future use of the global database of protected area management effectiveness. Phil. Trans. R. Soc. Lond. B 370, 20140281 (2015).
Google Scholar 

106.
Hansen, M. C. et al. High-resolution global maps of 21st-century forest cover change. Science 342, 850–853 (2013).
ADS  CAS  Google Scholar 

107.
Venter, O. et al. Sixteen years of change in the global terrestrial human footprint and implications for biodiversity conservation. Nat. Commun. 7, 12558 (2016).
ADS  CAS  PubMed  PubMed Central  Google Scholar 

108.
Geldmann, J., Joppa, L. N. & Burgess, N. D. Mapping change in human pressure globally on land and within protected areas. Conserv. Biol. 28, 1604–1616 (2014).
Google Scholar 

109.
Wilkie, D. S., Bennett, E. L., Peres, C. A. & Cunningham, A. A. The empty forest revisited. Ann. NY Acad. Sci. 1223, 120–128 (2011).
ADS  Google Scholar 

110.
Volenec, Z. M. & Dobson, A. P. Conservation value of small reserves. Conserv. Biol. 34, 66–79 (2020).
Google Scholar 

111.
Nicholson, E. et al. Scenarios and models to support global conservation targets. Trends Ecol. Evol. 34, 57–68 (2019).
Google Scholar 

112.
Maron, M., Rhodes, J. R. & Gibbons, P. Calculating the benefit of conservation actions. Conserv. Lett. 6, 359–367 (2013).
Google Scholar 

113.
Schleicher, J. et al. Statistical matching for conservation science. Conserv. Biol. 34, 538–549 (2019).
PubMed  PubMed Central  Google Scholar 

114.
Ferraro, P. J. Counterfactual thinking and impact evaluation in environmental policy. New Dir. Eval. 2009, 75–84 (2009).
Google Scholar 

115.
Chandler, M. et al. Contribution of citizen science towards international biodiversity monitoring. Biol. Conserv. 213, 280–294 (2017).
Google Scholar 

116.
Convention on Biological Diversity. Long-Term Strategic Directions to the 2050 Vision for Biodiversity, Approaches to Living in Harmony with Nature and Preparation for the Post-2020 Global Biodiversity Framework. www.cbd.int/decision/cop?id=12268 (2018).

117.
Secretariat of the Convention on Biological Diversity. Global Biodiversity Outlook 4 (Secretariat of the Convention on Biological Diversity, 2014).

118.
McCarthy, D. P. et al. Financial costs of meeting global biodiversity conservation targets: current spending and unmet needs. Science 338, 946–949 (2012).
ADS  CAS  Google Scholar 

119.
Balmford, A. et al. Walk on the wild side: estimating the global magnitude of visits to protected areas. PLoS Biol. 13, e1002074 (2015).
PubMed  PubMed Central  Google Scholar 

120.
Waldron, A. et al. Reductions in global biodiversity loss predicted from conservation spending. Nature 551, 364–367 (2017).
ADS  CAS  Google Scholar 

121.
Murray, K. A., Allen, T., Loh, E., Machalaba, C. & Daszak, P. Emerging Viral Zoonoses from Wildlife Associated with Animal-Based Food Systems: Risks and Opportunities (Springer, 2016).

122.
Dobson, A.P. et al. Ecology and economics for pandemic prevention. Science 369, 379–381 (2020). 

123.
Burmester, B. Upgrading or unhelpful? Defiant corporate support for a marine protected area. Mar. Policy 63, 206–212 (2016).
Google Scholar 

124.
Larson, E. R., Howell, S., Kareiva, P. & Armsworth, P. R. Constraints of philanthropy on determining the distribution of biodiversity conservation funding. Conserv. Biol. 30, 206–215 (2016).
Google Scholar 

125.
Smith, T. et al. Biodiversity means business: reframing global biodiversity goals for the private sector. Conserv. Lett. 13, e12690 (2019).
Google Scholar 

126.
Elsen, P. R., Monahan, W. B., Dougherty, E. R. & Merenlender, A. M. Keeping pace with climate change in global terrestrial protected areas. Sci. Adv. 6, eaay0814 (2020).
ADS  PubMed  PubMed Central  Google Scholar 

127.
Poloczanska, E. S. et al. Global imprint of climate change on marine life. Nat. Clim. Chang. 3, 919–925 (2013).
ADS  Google Scholar 

128.
Bruno, J. F. et al. Climate change threatens the world’s marine protected areas. Nat. Clim. Chang. 8, 499–503 (2018).
ADS  Google Scholar 

129.
Schleuning, M. et al. Ecological networks are more sensitive to plant than to animal extinction under climate change. Nat. Commun. 7, 13965 (2016).
ADS  CAS  PubMed  PubMed Central  Google Scholar 

130.
Bonnot, T. W., Cox, W. A., Thompson, F. R. & Millspaugh, J. J. Threat of climate change on a songbird population through its impacts on breeding. Nat. Clim. Chang. 8, 718–722 (2018).
ADS  Google Scholar 

131.
Hoegh-Guldberg, O., Poloczanska, E. S., Skirving, W. & Dove, S. Coral reef ecosystems under climate change and ocean acidification. Front. Mar. Sci. 4, 158 (2017).

132.
Jones, K. R., Watson, J. E. M., Possingham, H. P. & Klein, C. J. Incorporating climate change into spatial conservation prioritisation: a review. Biol. Conserv. 194, 121–130 (2016).
Google Scholar 

133.
Green, A. L. et al. Larval dispersal and movement patterns of coral reef fishes, and implications for marine reserve network design. Biol. Rev. Camb. Philos. Soc. 90, 1215–1247 (2015).
Google Scholar 

134.
Krueck, N. C. et al. Incorporating larval dispersal into MPA design for both conservation and fisheries. Ecol. Appl. 27, 925–941 (2017).
Google Scholar 

135.
van Kerkhoff, L. et al. Towards future-oriented conservation: managing protected areas in an era of climate change. Ambio 48, 699–713 (2019).
Google Scholar 

136.
Ling, S. D. & Johnson, C. R. Marine reserves reduce risk of climate-driven phase shift by reinstating size- and habitat-specific trophic interactions. Ecol. Appl. 22, 1232–1245 (2012).
CAS  Google Scholar 

137.
Maxwell, S. L., Venter, O., Jones, K. R. & Watson, J. E. M. Integrating human responses to climate change into conservation vulnerability assessments and adaptation planning. Ann. NY Acad. Sci. 1355, 98–116 (2015).
ADS  Google Scholar 

138.
Bennett, J. R. et al. When to monitor and when to act: value of information theory for multiple management units and limited budgets. J. Appl. Ecol. 55, 2102–2113 (2018).
Google Scholar 

139.
Burgass, M. J., Halpern, B. S., Nicholson, E. & Milner-Gulland, E. J. Navigating uncertainty in environmental composite indicators. Ecol. Indic. 75, 268–278 (2017).
Google Scholar 

140.
Bennett, J. R. et al. Polar lessons learned: long-term management based on shared threats in Arctic and Antarctic environments. Front. Ecol. Environ. 13, 316–324 (2015).
Google Scholar 

141.
Hughes, T. P. et al. Global warming and recurrent mass bleaching of corals. Nature 543, 373–377 (2017).
ADS  CAS  Google Scholar 

142.
Bai, Y. et al. Developing China’s ecological redline policy using ecosystem services assessments for land use planning. Nat. Commun. 9, 3034 (2018).
ADS  PubMed  PubMed Central  Google Scholar 

143.
Hughes, A. C. Understanding and minimizing environmental impacts of the belt and road initiative. Conserv. Biol. 33, 883–894 (2019).
Google Scholar 

144.
Alamgir, M. et al. High-risk infrastructure projects pose imminent threats to forests in Indonesian Borneo. Sci. Rep. 9, 140 (2019).
ADS  PubMed  PubMed Central  Google Scholar 

145.
Azevedo, A. A. et al. Limits of Brazil’s forest code as a means to end illegal deforestation. Proc. Natl Acad. Sci. USA 114, 7653–7658 (2017).
ADS  CAS  Google Scholar 

146.
Simmonds, J. S. et al. Moving from biodiversity offsets to a target-based approach for ecological compensation. Conserv. Lett. 13, e12695 (2020).
Google Scholar 

147.
Spalding, M. D., Agostini, V. N., Rice, J. & Grant, S. M. Pelagic provinces of the world: a biogeographic classification of the world’s surface pelagic waters. Ocean Coast. Manage. 60, 19–30 (2012).
Google Scholar 

148.
NatureServe. Bird Species Distribution Maps of the World (BirdLife International, 2018).

149.
Hijmans, R. J., Cameron, S. E., Parra, J. L., Jones, P. G. & Jarvis, A. Very high resolution interpolated climate surfaces for global land areas. Int. J. Climatol. 25, 1965–1978 (2005).
Google Scholar 

150.
Pauly, D. et al. Sea Around Us Concepts, Design and Data. www.seaaroundus.org (2020).

151.
Ferraro, P. J. & Pressey, R. L. Measuring the difference made by conservation initiatives: protected areas and their environmental and social impacts. Phil. Trans. R. Soc. Lond. B 370, 20140270 (2015).
Google Scholar 

152.
Díaz, S. et al. Pervasive human-driven decline of life on Earth points to the need for transformative change. Science 366, eaax3100 (2019).
Google Scholar  More