Oct 08, 2025
Version 1
Rice-fish coculture enhances agricultural sustainability and biological control V.1
- Nian-Feng Wan1
- 1Shanghai Key Laboratory of Chemical Biology, State Key Laboratory of Bioreactor Engineering, School of Pharmacy, East China University of Science and Technology, Shanghai 200237, China
Protocol Citation: Nian-Feng Wan 2025. Rice-fish coculture enhances agricultural sustainability and biological control. protocols.io https://dx.doi.org/10.17504/protocols.io.e6nvw4odwlmk/v1
License: This is an open access protocol distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited
Protocol status: Working
We use this protocol and it's working
Created: September 26, 2025
Last Modified: October 08, 2025
Protocol Integer ID: 228262
Keywords: Meta-analysis, rice-fish coculture, biodiversity gradient, PRISMA, comparing rice monoculture, findings of meta, rice monocultures to rice, analysis in rice, fish coculture, fish coculture system, agricultural sustainability, rice, meta, preferred reporting items for systematic review, systematic review
Funders Acknowledgements:
National Natural Science Foundation of China
Grant ID: 32472634
National Key R&D Program of China
Grant ID: 2023YFF1205101
Abstract
Protocols describing data describing data collection for a global meta-analysis comparing rice monocultures to rice-fish coculture applying PRISMA (Preferred Reporting Items for Systematic reviews and Meta-Analyses) guidelines. Supporting experimental work protocols to verify findings of meta-analysis in rice-fish coculture systems.
Guidelines
NA
Materials
Rice seedlings,
cylindrical plastic pot (44 cm top diameter, 37 cm bottom diameter, 52 cm height),
paddy field soil,
compound fertilizer (N: P2O5: K2O=15: 15: 15),
aquatic plant mud (Nitley, particle size 1-4mm),
Oujiang color common carp,
crucian carp,
Pardosa pseudoannulata spiders,
Nilaparvata lugens planthoppers.
Troubleshooting
Safety warnings
NA
Ethics statement
The protocol largely relates to a meta-analysis, but includes behavioral feeding experiment undertaken with fish to understand their preferences for different insect prey common in rice agriculture. The following ethics statement supports the use of fish in this experiment. Note, that while the East China University of Science and Technology does have an ethics committee it does not cover fish used as food. The ethics statement below therefore was not passed though this committee.
1. Justification for the use of animals
This experiment investigates the predatory behavior of two carp species, Oujiang color common carp (Cyprinus carpio var. color) and crucian carp (Carassius auratus), on brown planthoppers (Nilaparvata lugens) and wolf spiders (Pardosa pseudoannulata) in rice pot systems. Fish are used because they represent natural biological control agents in rice–fish systems, and no in vitro or computational alternative exists to study their foraging behavior under ecologically relevant conditions.
2. Number of animals and replication
A total of 40 individual fish will be used across all treatments (20 common carp - Cypinius carpia and 20 crucian carps - Carassius auratus). Observations are based on a 20-minute behavioral trials, replicated five times per condition (with 10, 20, or 30 planthoppers, or 3 male/female spiders, introduced per trial). The sample size is based on standard practice for behavioral ecology experiments and ensures adequate replication while minimizing the number of fish used. Replication was based on five treatments were: T1 = Fish (1 x common carp + 1x crucian carp) + 10 planthoppers; T2 = Fish + 20 planthoppers; T3 = Fish + 30 planthoppers T4 = Fish + 3 female spiders; T5 = Fish + 3 male spiders. Each of the five treatment was replicated four times for each of the two fish species. Each replicate observation was based on five separate individual behavioral observation.
3. Housing and acclimation
Fish will be acclimated for seven days in an indoor recirculating aquaculture system prior to trials, under controlled water temperature, aeration, and daily monitoring. They will be held in groups under appropriate stocking densities and checked daily for health and welfare.
4. Fasting procedure and justification
Fish will undergo a five-day fasting period before experimentation. The purpose of this fasting is to standardize hunger levels across individuals, ensuring natural foraging motivation and reducing variability in behavioral responses. A fasting period of this duration is commonly used in fish behavioral ecology and aquaculture research and is not harmful to healthy adult carp of the body weight used in this study (250–300 g). Both Oujiang color common carp and crucian carp are hardy species that tolerate short-term fasting without adverse effects, as they routinely experience natural periods of low food availability in field conditions. During fasting, fish will be monitored daily for health status, and any individuals showing abnormal behavior or signs of distress will be excluded from the experiment and returned to normal feeding.
5. Experimental procedure
Trials will be conducted in cylindrical pots containing soil and rice plants at the heading stage. Each pot will be filled with water (10 cm depth) and covered with aquatic mud substrate to maintain water clarity and reduce stress. A single fish will be introduced to each pot, together with a set number of planthoppers or spiders, and behavior will be recorded for 20 minutes using a video camera. No person will be present during recording to avoid disturbance. At the end of each trial, fish will be carefully removed and returned to their housing tanks. Surviving spiders and planthoppers will be collected and counted.
6. Animal welfare considerations
- No invasive procedures or harmful interventions will be performed.
- Trials are short (20 minutes), with minimal handling.
- Fish will be monitored for signs of stress before, during, and after trials.
- No fish will be sacrificed; all will be maintained in good condition following the experiments.
- Humane endpoints: if a fish displays abnormal behavior, injury, or distress (e.g., loss of equilibrium, erratic swimming), the trial will be terminated immediately and the animal will be removed.
7. Disposal of animals
At the conclusion of the study, all fish will be returned to holding tanks and maintained in standard aquaculture conditions.
Before start
NA
Study selection
Definition of rice-fish coculture and the number of added fish species
Here we considered “rice-fish coculture” as any increase in the number of fish species, relative to a rice monoculture (i.e., control group without fish). Namely, “rice-fish coculture” is a binary variable (zero or one), indicating whether fish species are present, irrespective of the number of fish species added. Here “number of added fish species” is a continuous variable describing the increase in fish species between the control (i.e., rice monocultures without fish species) and the rice-fish cocultures (i.e., rice-fish coculture with a number of fish species). When comparing rice monocultures with the treatment (i.e., rice-fish coculture, with≥1 fish species in rice fields), we ensured that comparisons were also made between the same trophic groups for rice, herbivores, rice diseases, and weeds.
We conducted a literature search on the Web of Science platform (i.e., Web of Science Core Collection, BIOSIS Previews, Derwent Innovations Index, KCI-Korean Journal Database, MEDLINE, Preprint Citation Index, ProQuestTM Dissertations & Theses Citation Index, and SciELO Citation Index), China Wangfangdata (https://wanfangdata.com.cn) and China National Knowledge Internet (www.cnki.net) (last accessed in late September 2025). We used the Boolean search string based on the “TOPIC” searching: [“rice-fish”] AND [“predat*” OR “parasit*” OR “wasp*” OR “natural enem*” OR “herbivor*” OR “pest” OR “biological control” OR “yield”]. In total, the search yielded 47808 articles. The titles and abstracts of this overall set of articles were first screened to determine whether the articles had measured a response variable relating to the effects of rice-fish co-culture on i) the growth, reproduction and quality of rice; ii) the growth and reproduction of invertebrate herbivores; iii) the rice disease damage; iv) the growth and reproduction of weeds; and v) the abundances of invertebrate predators and parasitoids. In total, 47361 articles were excluded at this stage of title and abstract screening. From the remaining 447 papers, we obtained our final selection of 110 articles based on the following criteria: i) the means, standard deviations (or standard errors) and sample sizes of the selected variables could be extracted from the text, tables, figures or supplementary information; ii) where pesticides and fertilizers were used, the amount and identity/type should be the same across all treatments (rice-fish co-culture, ≥1 fish species) and control (rice monoculture). In the cases where no standard deviations, standard errors, or confidence intervals were reported, we calculated standard deviations by Eq.6 which was referred as “All Cases” method, according to Nakagawa et al (2023). Articles were screened by N.F.W., S.S., M.L. and Y.C. Data were extracted from the articles by N.F.W., S.S., F.Y.P. and Q.L during which regular cross-checking was performed to ensure consistency in extracted effect sizes.
We used data giving the mean values of multiple sampling dates or years. If these mean values were not presented, we used the data of the latest sampling period (Wan et al., 2022a; Wang et al., 2025). For articles that covered more than one experimental location, we considered these experimental results separately (see locations in Fig. 1a). When numeric values were not provided directly, we extracted them from figures using the “GetData Graph Digitizer” 2.26. When we had questions in the process of extracting data, we contacted the authors directly. To avoid pseudoreplication of data, we excluded multiple comparisons conducted within a single experiment (Wan et al., 2020). Observations in rice monoculture were considered as control groups, while those in rice-fish coculture were considered as the treatment groups. When an article included different levels of fish species, measurements for the control groups without fish were compared to all other treatments with different levels of fish species and treated as independent paired observations. When extracting the data of rice yield in rice-fish co-culture systems, we removed the water ditch area without rice plants. For example, in an area of 500m2 rice-fish co-culture plots with 25m2 water ditch area (i.e., 5% of non-planted area of the total field), rice yield is 4.30 Mg•ha-1, and we re-calculate the rice yield by transforming it as 4.53 Mg•ha-1 (i.e., 4.30÷0.95=4.53).
A PRISMA (Preferred Reporting Items for Systematic Reviews and Meta-Analyses) describing the progress of data acquisition through the meta-analysis. It’s the records identified in terms of numbers, and subsequent processed on inclusions and exclusions to derive the final data set used to assess the effects of rice-fish co-culture on the performance of rice, pests and their natural enemies.
Predictor variables
Five categorical predictor variables and one continuous predictor variable were derived in association with each loge response ratio effect size (lnRR; see below for derivation).
Trophic group: A categorical variable indicates whether the organisms whose lnRR effect size responses were as primary producers or different classes of pests. Rice was the focus of the biodiversity-productivity relationship under investigation. Where lnRR effect sizes described pests, these were classified by an aggregate categorical variable (i.e., rice pests) including invertebrate herbivores, rice diseases and weeds. A further category classification distinguished between the sub-groups of invertebrate herbivores, rice diseases that infested rice and cause damage to rice, or weeds (i.e. pest plants in rice fields, competing for water, light and nutrients with rice). Moreover, an aggregate categorical variable (i.e. pests) including invertebrate herbivores, rice diseases and weeds.
Response category: Response categories include i) growth, reproduction and quality of rice (i.e., rice performance); ii) damage and reproduction of herbivores (i.e., herbivore performance); iii) abundances of predators and parasitoids (i.e., natural enemy performance); iv) rice disease damage (i.e., disease performance); and v) weed growth and reproduction (i.e., weed performance). We divided the lnRR effect size response variables into multiple categories nested within each trophic group. These were as follows:
Rice performance: This represents the principal metric used in the assessment of the diversity-productivity relationship. As such this was the metric exposed to the negative effects resulting from the activity of the pests. Note the effect of these pests on rice performance was directly assessed within the path analysis framework described below. Measures of rice performance included rice growth, reproduction and quality. Rice growth was assessed as growth rate and size characteristics (e.g., plant weight, plant height, culms height, aboveground biomass, dry matter biomass of rice plants, leaf area, leaf biomass, root length, number of tillers per plant). Rice reproduction included the yield of grains, reproductive plant parts and seeds per unit, and reproductive traits measured on individual plants (e.g., grain yield, number of panicles, panicle weight, thousand kernel weight, number of spikelets, filled grain number per panicle, seed-setting rate, and rice kernel set percentage). Rice quality included rice appearance, taste value, brown rice rate, head rice rate, milled rice rate, and protein content percentage, referring to The Ministry of Agriculture of the People's Republic of China (2010).
Pest herbivore performance: This included mollusk (i.e. snails) and invertebrate herbivore reproduction (e.g., herbivore population abundance and survival rate), and herbivory damage (e.g., percentage or grade of herbivore damage to rice). Rice disease performance: This included disease damage (e.g., disease index of rice sheath blight, incidence of rice blast, incidence of sheath blight disease, percentage of stripe blight in rice clusters, and severity of leaf blight disease). Weed performance: This was composed of weed growth (e.g., weed biomass, dry matter weight, and fresh weight), and weed reproduction (e.g., weed population/community abundance or density). Studies on plant-feeding nematodes and plant-feeding rodents were not found in this meta analysis. Overall metric of pest performance: This was an overall metric which integrated lnRR effect sizes for the performances of invertebrate herbivores, rice diseases, and weeds. This was done by binding the observations of herbivore, disease, and weed performances.
Climatic zone type: A categorical variable reflecting whether a certain study was performed in the temperate or tropical zones. Temperate zones ranged from 23.5°N to 66.5°N and from 23.5°S to 66.5°S, and the tropical zones ranged from 0−23.5°N and from 0−23.5°S. Data from greenhouse and indoor experiments were removed from models considering climatic predictors (Wan et al., 2022a).
Type of experimental study: This variable describes where the studies were used plot-based experiments where rice was planted directly into ground soil (e.g. in field-based plots) or pot-based experiments with rice planted into soil in above vessels. Experimental design is not or direct applied or policy interest, but has significant potential to result in artifacts in the analysis linked to the impacts of confined model systems in pots that may not be encountered under field conditions. By including this covariate we aim to assess the extent of such possible confounding effects. Common garden experiments with a few or several replicated plots and open field experiments in rice fields, were
considered as plot experiments. Pots, containers, bottles, trays, boxes and tankers experiments were considered as pot experiments. Greenhouse experiments were considered as plot or pot experiments. Indoor and laboratory experiments were considered as pot experiments.
Rice field types - organic or non-organic: A categorical variable indicating whether experiments employed rice-fish co-culture or rice monoculture systems, categorized as organic or non-organic. Organic systems strictly prohibited synthetic agrochemicals (including chemical pesticides and fertilizers), while non-organic systems may have used such inputs.
Number of added fish species over the control: A continuous variable describing the increase in fish species between the control (rice monocultures without fish) and the treatment containing an increased number of fish species relative to the control. When we compared the fish species of the control (i.e., rice fields without fish) with that of the treatment (i.e., rice-fish co-culture, ≥1 fish species), we confirmed that both the control and treatment were compared on the same trophic groups among the invertebrate herbivores, rice diseases, weeds and rice.
Meta analysis for rice yield
We derived a new dataset including paired observations of rice yield in rice-fish coculture vs. rice monoculture in each experiment. In total, among all the 113 studies, we extracted the data of 77 studies with paired 121 observations. Firstly, we calculated the increased percentage of rice yields. Increased percentage of rice yield = (rice yield in rice-fish coculture－rice yield in rice monoculture) ÷rice yield in rice monocultur×100%. Units of rice yield include ton·ha-1, kg·ha-1, Mg·ha-1, kg·plot-1, kg·666.7m-2, kg·m-2, kg·5 dot-1, g·m-2, kg·25 m-2, and q·ha-1.
Definition of effect size and its measures: To test the effect of fish species on the various groups (growth, reproduction and quality of rice, damage and reproduction of herbivores, abundances of predators and parasitoids, rice disease damage, weed growth and reproduction ; and moreover, an aggregate pest categorical variable including herbivores, rice diseases, weeds, predators and parasitoids), we calculated the effect size and lnRR of these groups. The first proposed formula as follows:
where m1 and m2 are the observed mean value in the treatment and control groups, sd1 and sd2 are the standard deviations (SDs) in the treatment and control groups, and n1 and n2 are the sample sizes in the treatment and control groups. m1, m2, sd1, sd2, n1 and n2 were extracted from original articles. Namely, mean value, SD and sample sizes of both treatment and control groups were included in our dataset to conduct meta-regression. This was done in order to deal with missing standard deviations (SDs) in dataset, using the approach of Nakagawa et al. (2023a), which is suited to accounting for missing SDs. Nakagawa et al. (2023a) proposed this new method weighting average coefficients of variation estimated from studies that do report SDs in the dataset. This is done by:
where CV=sd/m is the coefficient of variation; sd and n are the corresponding SDs and sample size, respectively; CV1i and CV1i are the CVs from the ith study (study: i=1,2,…,K). The proposed Eq.3 and Eq.5 can improve the accuracy and precision of the overall mean estimate. We can use Eq.5 and 6 to calculate the effect sizes and sample variances when SDs are missing and use Eq.1 and 2 to calculate the effect size and sample variances when SDs are not missing. This paper uses Eq.5 to calculate effect when the dataset includes missing SDs, whereas uses Eq.1 when the dataset do not include missing SDs. We used lnRR as the response variable for different models except for path analyses described below.
Meta-regression
Meta-regression (van Houwelingen et al., 2002) was applied to test whether the effect sizes of different trophic groups could be explained by increasing fish species and the various predictor variables. Specifically, we fitted three-level mixed-effects meta-regression models, using the R package metafor (version 3.8-1). The effect size metric lnRR was calculated using the function “lnrr_laj()” in R.file implemented using the package “func.R” developed by Nakagawa et al. (2023a). This was used to calculate the effect size metric lnRR for each observation as well as the unbiased sample variance estimates as defined under Eq.6 (using the function “v_lnrr_laj()” in R.file “func.R”). Fish species richness, trophic groups, trophic group response categories, climatic zones, study types, and rice field types, were included as moderators whose effects were assumed to be fixed. We also constructed a null model without any predictors. In all models, we treated ‘study identity’ as a random effect. To handle non-independence in effect sizes of each study, we added a random effect as a unique identifier for each effect size in every study (EsID), which allows true effect sizes to vary within studies, and to account for the within-study effect and quantify within-study heterogeneity. To ensure robust estimates and account for differences in precision across studies and effect sizes, we weighted effect sizes by the inverse of the sum of the variance-covariance matrix components. This approach explicitly incorporates the non-zero covariances that arise from correlations in sampling variance within the same source articles, as well as from random-effect variance (Nakagawa et al., 2023b). Phylogenetic correction for fishes was also undertaken to investigate the effects of fish species on trophic groups. To do this, fishes species phylogenies were included as a random effect with phylogenetic relatedness as part of the correlation structure. Specifically, we matched the fish species included in our analysis with the available synthetic tree in R package “rotl”, then the relationships between each matched fish species were returned so that they could be used to draw phylogenetic trees. To adjust for repeated measurement of control values, we assigned the argument “V” in “rma.mv()” function with the sampling variance-covariance matrix estimated as follows (Nakagawa et al., 2023b):
where the subscripts 1T and 2T represent the different treatment group who share the control value in a same group; the subscript 1C represents the 1st control group; n represents the number of observations, CV represents the coefficients of variance of observed value. Taking this approach, 24 fish species defined the correlation structure based on the fish phylogeny (Jin & Qian, 2019, 2022).
For each mixed-effects meta-regression model, we first fitted a base model by treating fish species and trophic group as the fixed effect terms. Second, the interactions between the trophic group and other predictor variables (climatic zones, study types, and rice field types) were also included in the model to assess whether model fit was improved, using a likelihood-ratio test (LRT). Third, the trophic group response category (nested within trophic group) and the interactive effects between the response category and predictors were also included in the model using a LRT to allow model comparisons (Supplementary Table 1). For example, the model with “trophic group + study type” was compared to the base model with just trophic group, and the model with “trophic group + trophic group × study type” would be compared to a model with “trophic group + study type” (Supplementary Tables 1 and 2). To examine whether the mean effect sizes of (added) fish species, response category and other predictors differed significantly from zero, we acquired estimations with their 95% confidence intervals, which were derived from the fitted meta-regression models. To assess the between-study heterogeneity in these models, I2 statistics were calculated (Ruppar et al., 2020; Nakagawa et al., 2023b).
Meta-regression models for model comparison and parameter evaluation
We performed meta-regression analysis using the R package metafor (version 3.8-1), in which we took log ratio response (lnRR) as the effect size measure (in the function “lnrr_laj()” in R.file “func.R” developed by Nakagawa et al (2023a). Unbiased sample variance estimates were constructed to evaluate the variances of lnRRs (in the function “v_lnrr_laj()” in R.file “func.R”). To avoid the potential bias induced by repeated control values, we set the sampling covariance matrix according to Eqs.7–9. We employed a likelihood ratio test to compare full model with null models to investigate the significance of the interactive effects between trophic groups or integrated trophic groups and various moderators. This refers to trophic groups (invertebrate herbivores, rice diseases, weeds, predators, parasitoid and rice), response categories (growth, reproduction and damage of herbivores; damage of rice diseases; weed growth, reproduction; abundance of predators and parasitoid; and growth, reproduction and quality of rice) and different
moderators (i.e. types of study, climatic zones, rice field types, and number of added fish species).
The null model corresponded to the equation:
The full model corresponded to the equation:
Where βij the jth effect size of study i, ri is the random effect represents the heterogeneity between studies or the phylogenetic relatedness of the fish species within ith study, and wi is the random effect indicates the heterogeneity of the jth effect size within ith study. In addition, εij is the random error term with variance equal to the lnRR’s sample variance estimate. In all cases trophic group response categories were nested within trophic groups. For example, invertebrate herbivore response categories (i.e., herbivore growth, herbivore reproduction, and herbivore damage) were nested in the invertebrate herbivore group. See above for description of these for each of the antagonists. Additionally, response categories of herbivore growth, herbivore reproduction, herbivore damage, rice disease, weed growth, weed reproduction, predator abundance and parasitoid abundance were nested into the overall pest grouping. Using this approach relevant response category effects could have an impact on the corresponding trophic group effect.
We applied the likelihood ratio tests (LRT) to test the null model with the full model to explore the interactions between the trophic group response categories (herbivore growth, herbivore reproduction, herbivore damage, rice disease, weed growth, weed reproduction, predator abundance, parasitoid abundance, rice growth, rice reproduction, rice quality, pest) and the different moderators. The null model was as follows:
The full model followed the equation:
Subsequently, to estimate the effects and their corresponding confidence intervals of each trophic group, trophic group response category and their associations with the moderators, the following models were used:
Phylogenetic covariance structures for fish
Phylogenetic relationships between the species have the potential to bias the meta-analysis as closely related species show greater commonality in their responses to increasing fish species affecting assumptions of independence. However, accounting for this within the current modelling framework required us to focus on only those lnRR effect sizes that were based on a comparisons of single species monoculture controls compared to a higher species treatment. For this sub-set of the data, it would be possible to identify where phylogenetic similarities between species in monocultures affected the overall community response to further increases in fish species, either for rice or pest performance. From the 113 articles we obtained 791 observations from 466 articles had a single species control from which the lnRR effect size was derived. Within this dataset the control groups included 24 fish species. We generated a phylogeny for these fish species based on a mega-tree of 74,531 species of extant vascular plants representing the largest phylogeny for vascular plants as implemented in rotl (Hinchliff et al., 2015). All 24 fish species were represented in this overall phylogeny. Following this the potential effect of fish species was accounted for in all models by 1) adding random intercepts for fish species identity, and 2) including the phylogenetic relatedness between fish species as part of the correlation structure (Smith & Brown, 2018). Furthermore, each effect size was nested within the corresponding study to incorporate the hierarchical error structure of multiple effects coming from the same study.
Meta-analytic path analysis via residual regressions
To explore the trophic interactions between plants and various groups of pests, we established a new data subset: 1) comprising paired trophic observations (i.e., rice performance vs. herbivore performance, rice performance vs. rice disease performance, rice performance vs. weed performance, and rice performance vs. pest performance) ; 2) encompassing the paired observations of rice performance vs. pest performance in different climate zones, rice field types and study types; 3) and comprising the paired observations of rice performance vs. herbivore performance vs. natural enemy performance . The lnRR effect sizes were derived from these pairwise data sets.
To analyze the direct effect of pests on rice performance in path analysis, residual regression was applied. In this analysis we used the lnRR as the response variable. Specifically, the direct effects of increasing fish species on pest were estimated using linear mixed-effects model:
Following Emmenegger & Bühlmann (2023). The direct effect of fish species on rice performance were estimated by regressing the Pearson residuals of the linear mixed-effects model:
on fish species. The direct effects of pests on rice were estimated by regressing the Pearson residuals extracted from the linear mixed-effects model:
on the Pearson residuals extracted from the linear mixed-effects model:
where β is the effect size, ri is the random effect represents the heterogeneity between studies or the
phylogenetic relatedness of the fish species within ith study, εi is the random error term with variance equal to the lnRR’s sample variance estimate.
We did above path analysis to understand the relative effect of i) increasing fish species on the performance of the overall animal community, ii) the performance of the pest and iii) the indirect effect on rice performance resulting from the response by the pest to increasing fish species. Specifically, the following models were constructed as input to residual regressions, and the β0 and β1 below were the effect of fish species and the effect of added fish species respectively:
In summary, for our path analyses, linear mixed-effects models were implemented within the R statistical environment using the function “lme” of the package “nlme”, with random intercepts for study IDs. Heteroscedasticity was accounted for by providing fixed variances based on lnRRs and setting sigma to 1 in the lme call.
We applied Student’s t-tests to assess statistical significance within the path analysis for each of the above equations. A significant positive intercept shows that fish species increases the mean value of the pest performance or rice performance for a trophic group and a negative intercept shows the reverse. We extracted the z-values of corresponding coefficients to test the effects of fish species on rice and different groups of pests (i.e., herbivores vs. rice, rice disease vs. rice, weeds vs. rice, and natural enemy vs. plants) . Likewise, we extracted the z-values of corresponding coefficients to test the effects of fish species on the interactions between performance values of rice and the aggregate indicator (i.e., pests) in different climatic zones (temperate and tropical zones - excluding greenhouse or indoor observations), study types (plot and pot experiments), and rice field type (organic and non-organic rice field), respectively. We explored the effects of number of added plant species relative to the control by constructing various models to explore the direct and indirect effects of added fish species on interactions. Likewise, we conducted the direct and indirect effects of the number of added fish species on rice and pests distinguishing between climatic zones, rice field type and in study type . The estimations and test statistics were extracted using the R function “coef()”. The relative goodness-of-fit analyses for path analyses of predictor variables (fish species), were conducted by extracting AIC (Akaike information criterion), AICc (corrected Akaike information criterion), BIC (Bayesian information criterion) and log-likelihood from the fitted models, using R functions “AIC()”, “AICc()”, “BIC()” and “logLik()” from the rsq package (version 2.5) (Zhang, 2022).
In the above path analysis for top-down effects in bi-trophic , we used mediation analysis to test whether the effect of an independent variable RF (i.e., rice-fish coculture vs. rice monoculture) on a dependent variable Y' (i.e., herbivore performance, including herbivore damage and reproduction) (RF→Y') is at least partly explained via the inclusion of a mediator variable X' (i.e., pest performance, including herbivore damage, herbivore reproduction, disease damage, predator and parasitoid abundance) (RF→X'→Y'). In our top-down analysis, the 3 causal paths E1, E2 and E3 correspond to the effect of rice-fish coculture on pest performance, the effect of pest performance on rice performance and the effect of rice-fish coculture on rice performance accounting for pest performance respectively. The three causal paths correspond to parameters from two regression models, one in which pest performance is the outcome and rice-fish coculture is the predictor, on in which rice performance is the outcome and rice-fish coculture and pest performance are the simultaneous predictors, one in which rice performance is the outcome and rice-fish coculture is the predictors. From these parameters, we can calculate the mediation effect (the product E1×E2) and total effect (E1×E2+ E3) of rice-fish coculture on rice performance. In tri-trophic, we test whether the effect of an independent variable RF (i.e., rice-fish coculture vs. rice monoculture) on a dependent variable Y (i.e., herbivore performance, including herbivore damage and reproduction) (RF→Y) is at least partly explained via the inclusion of a mediator variable X (i.e., natural enemy performance, including predator and parasitoid abundance) (RF→X→Y). The effect of an independent variable RF on a dependent variable Z (i.e., rice performance, including rice growth, reproduction and quality) is at least partly explained via the inclusion of a mediator variable Y (RF→Y→Z). In our top-down analysis, the five causal paths e1, e2, e3, e4, and e5 correspond to the effect of rice-fish coculture on natural enemy performance, the effect of natural enemy performance on herbivore performance, the effect of rice-fish coculture on herbivore performance accounting for natural enemy performance, the effect of herbivore performance on rice performance and the effect of rice-fish coculture on rice performance accounting for herbivore performance respectively. The five causal paths correspond to parameters from three regression models, one in which natural enemy performance is the outcome and rice-fish coculture is the predictor, on in which herbivore performance is the outcome and rice-fish coculture and natural enemy performance are the simultaneous predictors, one in which rice performance is the outcome and herbivore performance is the predictors and one in which rice performance is the outcome and rice-fish coculture and herbivore performance are the simultaneous predictors. From these parameters, we can calculate the mediation effect (the product e2×e3 and e3×e4) and total effect (e2×e3+ e1) of rice-fish coculture on herbivore performance and total effect (e3×e4+ e5) of rice-fish coculture on rice performance.
Publication bias test
We assessed publication bias using regression tests (Egger & Davey, 1997; Nakagawa et al 2022), which employ a partial slope test of association between effect size and the sample size. Here, a significant relationship (p<0.05) suggests publication bias. The trim-and-fill method was not employed as this is inappropriate for models with moderators (Duval & Tweedie, 2000; Nakagawa & Santos, 2012). Instead, we adopted the method suggested by Nakagawa et al. (2022), which uses a multilevel version of Egger regression to assess publication bias in mixed-effects meta-regression analysis. Following this approach, we fitted the following model to assess publication bias (where “sample size” represents the effective sample size that is calculated using):
Here, we considered lnRR as a response variable. Different categories of trophic groups, climatic zones, study type, and rice field types were considered as predictors, respectively, and sampling sizes were considered as an additional moderator in the mixed-effect model, and the test statistics for the coefficients of sampling sizes were used to test for publication bias. As a result, publication bias was assessed for all models by testing H0: β1=0 . The Egger’s regression test showed that there was publication bias for the 791 observations of the 113 cited articles in this paper (value of regression test = -0.96739, P = 0.334). We also derived the fail-safe N test to calculate the number of null-result studies that would have to be added to reduce the combined significance level (p-value) to a 0.05 alpha level. This fail-safe number can then be compared to a rule-of-thumb threshold, defined for this study to be 3965 (based on 5n + 10, where n is the number of observations in the analyses considered to be 791). This derived the fail-safe number of 17139 for the whole dataset of 113 cited articles. This means that at least 17139 studies with negative effects would be needed to be incorporated in our dataset to result in a non-significant overall effect. This provides strong evidence that the results in this paper are stable for publication bias.
Sensitivity tests to meta-analysis
We applied the trim and fill method to test the sensitivity. The trim and fill method is a nonparametric (rank-based) data augmentation technique proposed by Duval & Tweedie (2000). The method can be used to evaluate the number of studies missing from a meta-analysis based on the suppression of the most extreme results seen on one side of the funnel plot. The method then augments the observed data so that the funnel plot is symmetric and recomputes the summary estimate based on this derived data. Statistical values for this sensitivity analysis are presented for the different taxonomic groups of pest performance, invertebrate herbivore performance, rice disease performance, weed performance, natural enemy performance and rice performance from the overall 113 studies.
Field experiments
Field experiment design and data analysis
Field experiment was conducted at Situan Town, Fengxian District, Shanghai of China (30.94°N, 121.71°E). This site belongs to alluvial plains of the Yangtze River Delta in the east Asian monsoon zone with four distinctive seasons. The whole area of rice fields was about 46.7 ha. Field experiment was a randomized complete block design with three replicates. Each block consisted of one treatment plot (rice-fish Carassius auratus coculture; 90m×150m) surrounded on three sides by three ditches (0.7-0.9m in depth and 0.9m in width), and a control plot (rice monoculture; 90m×150m) without ditches or fish (Supplementary Figure 2). Blocks were spaced >50m apart to prevent spill-over effects among blocks. Two plots within each block were separated by >50m of buffer zone. Plots were established in late May 2020 when ditches of co-culture plots were filled with irrigation water. Field experiment was performed from 2021 to 2024.
In each year in late May, a late maturing variety of rice (Huruan 1212) was mechanical dibbled (cluster spacing was 16 cm × 25 cm) and seeded at 60 kg·ha-1. Each year plots received 15.0 ton·ha-1 of organic fertilizers as the base fertilizer in middle-to-late May, 225 kg·ha-1 of compound fertilizer (N: P2O5: K2O=15: 15: 15) at the seedling stage (early-to-middle June), 225 kg·ha-1 of compound fertilizer (N: P2O5: K2O=15: 15: 15) at the tillering stage (late June), and 180 kg·ha-1 of compound fertilizer (N: P2O5: K2O=15: 15: 15) at the booting stage (middle August).
During this experiment period from 2021 to 2024, to control rice disease and insect pests, each plot in each year received the same amount of the 24% Validamycin A Water Suspension (0.600kg·ha-1), 25% clothianidin water dispersible granules (0.300kg·ha-1), and 16% methiocarb-indoxacarb suspension concentrate (0.225 kg·ha-1). Weeds were removed by hands. Note that the application of these products was implemented to achieve consistency between treatments, but normal threshold based practices for pesticide application are in operation in this system under non-experimental conditions (See Wan et al., 2022b).
For each rice-fish coculture plot in late May 2020 when coculture plots were established, 1000 Carassius auratus (0.10 to 0.15 kg each) were released into the three surrounding ditches. Fish were not fed supplemental food to promote fish feeding on insect pests. To create ideal conditions for rice and fish growth, irrigation depth was maintained at 10–15 cm above the base of rice plants in coculture plots during the rice seedling stage; this avoided negative effects of flooding during the rice seedling stage and allowed small fish to move freely in rice fields. At the beginning of the tillering stage as abundance of insect pests and damage increased, irrigation depth was increased to 15–20 cm above the base of rice plants in co-culture plots, to provide a more suitable environment (i.e. sufficient water, enough activity space, reduced water temperature) for fish to prey on insect pests. In rice monoculture plots, depth of irrigation water was maintained at 3–5 cm. Irrigation was stopped 1 week before harvest for both plot types.
Harvest operations were conducted in middle-to-late October in each year from 2021 to 2024. In each rice plot, we collected the rice ears using a Z-style sampling method with 10 dots (20 rice ears in each dot were randomly selected). Ears were dried in the sun for 5 days, threshed to separate the ears into individual grains, and 1,000 grains were weighed to determine thousand kernel weight. Then a harvester was used to obtain the grain yield in each plot.
Insect pests and invertebrate predators were collected from 100 randomly selected rice clusters in each plot. Only three groups of rice insect pests (i.e., planthoppers, stemborers and leafrollers) which represent the main insect pests on rice in eastern China (Wan, et al., 2018), were sampled as target insect pests. As predatory spiders account for more than 90% of natural enemies and play an important role in pest control in rice fields (Tanwar et al., 2011), we sampled only spiders as target predators to test our hypothesis.
Insect pests and invertebrate predators were sampled at 2–4 week intervals from the seedling to maturity stages (i.e., from early to mid-July to early to mid-September), for a total of four sampling times per year. To sample, a white stainless-steel plate was placed at the base of rice plants (Wan, et al., 2018) then rice plants were shaken so that insect pests and invertebrate predators fell onto the plate and were counted. To further investigate the abundance of stemborers and leafrollers, we checked for infested rice stems and foliage for each sampled plant and counted the stemborers and leafrollers. Thus, abundance of rice stemborers or leafrollers at each sampling was the sum from the plate and from the infested stems or leaves (Wan et al., 2018).
Normality and homoscedasticity of all data were assessed by Kolmogorov–Smirnov and Levene’s tests. ANOVAs were used to test fixed effects of rice field type (rice-fish coculture and rice monoculture), experimental year (2021, 2022, 2023 and 2024), and their interaction, sampling period (1 to 4 within a year) and the interaction of treatment × year × period and random effect plot (treatment) on: abundance of herbivores (i.e., planthoppers, stemborers and leafrollers) and predatory spiders. Experimental year (2021, 2022, 2023 and 2024) and experiment sample period (1 to 4 within a year) were considered independent factors. Models without terms related to sampling period were used for rice thousand kernel weight and rice yield. Adjusted means contrasts were used to compare the differences among means for predictors with more than two combinations or levels. If the values of herbivore or spider abundance were significantly affected by the years, we considered the four years’ data separately to compare the differences among means of rice monoculture and rice-fish coculture with Tukey’s Honestly Significant Difference (HSD) test at 0.05 level. We use “aov()” function to conduct two-way ANOVA analyse, with term “abundance ~ year * treatment” in R.
In terms of path analysis, we using structure equation model (SEM) to do this. Our SEM path analyses were fitted with “psem” function in library “piecewiseSEM” (Lefcheck, 2016) using series of linear mixed-effects models with random intercepts experiment year. Specifically, we extracted the t-values and its corresponding coefficient to test the effects of fish species on the tri-trophic interactions of rice, herbivore, natural enemy performances (Supplementary Data 5). The estimations and test statistics were extracted using the “summary()” function in R on piecewise equation SEM models
Behaviour experiment of fish preying on brown planthoppers and their predatory spiders
To further explain above mechanisms, we conducted a behaviour experiment to test the predatory effect of fish on spiders and rice planthoppers, by using records of a video camera in 2025. We used two fish species: one is an indigenous, red, soft-scaled, common species of carp (Oujiang color common carp, Cypinius carpia color var.) from Qingtian which is a “globally important agricultural heritage system” (GIAHS) as designated by the United Nations Development Program and the Global Environment Facility in 2005 (Xie et al., 2011), and the other is a native species (crucian carp, Carassius auratus) from Situan, Shanghai, China (see above field experiment). We used the wolf spider (Pardosa pseudoannulata), collected from Jinggangshan Agricultural Science and Technology Park, Xingqiao Town, Jizhou District, Ji'an City, Jiangxi Province, China (27.14°N, 114.88°E) (provided by the research team of Professor Yonghong Xiao, Jinggangshan University, China)as this species is the most important dominant spider species in China's rice fields (Huang et al., 2018). High-quality fragrant japonica rice variety “Huxiangjing 106” was used, as this variety has relatively weak resistance to brown planthopper (Nilaparvata lugens), provided by Professor Jinglan Liu at Yangzhou University of China. Oujiang color common carps and crucian carps (body weight: 250–300g) were acclimated in an indoor recirculating aquaculture system for 7 days prior to the experiment and subjected to a 5-day fasting period to standardize their hunger levels and stimulate natural foraging motivation.
Rice seedlings were planted in pots. The potting container is a cylindrical plastic pot with a top diameter of 44 cm, a bottom diameter of 37 cm, and a height of 52 cm. Each pot is filled with 25 kg of paddy field soil that has been air-dried, crushed and sieved. The base fertilizer is compound fertilizer (N: P2O5: K2O=15: 15: 15) at 8.33g·kg-1 of soil. When seedlings are at the three-leaf and one-heart stages, we transplant them into above pots. We plant a single clump of rice seedlings (2 seedlings per clump) in the center of each pot. During the growth period of rice, regular water management should be carried out and the use of any pesticides should be avoided.
We select potted rice plants at the heading stage that are vigorous and have uniform growth. To ensure the water is clear and facilitate subsequent observation and video recording, evenly cover the surface of each pot of soil with a 2 cm thick layer of aquatic plant mud (Nitley, particle size 1-4mm). This material can effectively adsorb suspended particles, prevent fish activities from muddying the water body, and is non-toxic and harmless to aquatic organisms. We slowly fill the basin with water until the water depth remains at 10 cm, and then install the video camera properly. In each experimental pot, we used 10, 20, or 30 third-instar brown planthoppers, or 3 female spiders, or 3 male spiders, respectively (Wu et al., 2002). The five treatments were: T1 = Fish (1 x common carp + 1x crucian carp) + 10 planthoppers; T2 = Fish + 20 planthoppers; T3 = Fish + 30 planthoppers T4 = Fish + 3 female spiders; T5 = Fish + 3 male spiders. Each of the five treatment was replicated four times for each of the two fish species. Each replicate observation was based on five separate individual behavioral observations.
After the prey (planthoppers and spiders) are evenly distributed in the pots, we carefully place one Oujiang color common carp or one crucian carp that has been starved for continuous five days in one pot. Then the camera was immediately activated to start recording the behavior of fish preying on rice planthoppers and spiders. If rice plants are planted in clusters within the pots, the spiders and rice planthoppers introduced will affect the camera’s view due to their presence on the rice stems. Therefore, the leaves of the rice plants in the pots should be removed and only the main stems should be retained to record the fish predation activity. The video recording lasts for 20 minutes, and each processing is repeated 5 times. No one needs to be present during the recording to avoid interference to fish activity. After the finish of recording, we carefully collect and count the number of all surviving brown planthoppers and spiders, and then calculate the predation percentage of fish on planthoppers or spiders. Predation percentage is the number of planthoppers or spiders eaten by fish in a pot divided by the ushered number of planthoppers or spiders in a pot × 100%.
We employed a two-way ANOVA under a general linear model framework to assess predation rates, examining the interactive effects of fish species (common carp vs. crucian carp) and planthopper density (10, 20, or 30 third-instar brown planthoppers). We use “aov()” function to conduct two-way ANOVA analyse, with term “predation rate ~ year * treatment” in R. Post hoc comparisons of means were performed using Tukey’s HSD test (α = 0.05), with Bonferroni correction for multiple testing. Prior to ANOVA, predation percentage data underwent arcsine-square root transformation to meet parametric assumptions.
We used R version 4.3.1 (R Core Team, 2023) to conduct all statistical analyses, and used R package “metafor” 3.8-1 to perform meta-regression and publication bias assessment (Viechtbauer, 2010). In addition, we used R packages “nlme” (Pinheiro, 2023) and “piecewiseSEM” (Lefcheck, 2016) to residual regression in path analyses. A significance level of 0.05 was used for all tests.
Protocol references
Duval, S. & Tweedie, R. Trim and fill: A simple funnel-plot-based method of testing and adjusting for publication bias in meta-analysis. Biometrics 56, 455–463 (2000).
Emmenegger, C., & Bühlmann, P. Plug-in machine learning for partially linear mixed-effects models with repeated measurements. Scand. J. Statist. 50, 1553–1567 (2023).
Egger, M., Davey, S. G., Schneider, M. & Minder C. Bias in meta-analysis detected by a simple, graphical test. Brit. Med. J. 315, 629–634 (1997).
Fox, J. & Weisberg, S. An R Companion to Applied Regression, 3rd edition. Sage, Thousand Oaks CA.
https://socialsciences.mcmaster.ca/jfox/Books/Companion/index.html
(2019).
Hinchliff, C. E. et al. Synthesis of phylogeny and taxonomy into a comprehensive tree of life. Proc.Nat. Acad. Sci.U.S.A. 112, 12764–12769 (2015).
Huang, L. X., Wang, Z. Y., Yu, N., Li, J. J. & Liu, Z. W. Toxin diversity revealed by the venom gland transcriptome of Pardosa pseudoannulata, a natural enemy of several insect pests. Comp. Biochem. Phys. D28, 172–178 (2018).
Jin, Y. & Qian, H. V.PhyloMaker: an R package that can generate very large phylogenies for vascular plants. Ecography 42, 1353–1359 (2019).
Jin, Y. & Qian, H. V.PhyloMaker2: An updated and enlarged R package that can generate very large phylogenies for vascular plants. Plant Diversity 44, 335–339 (2022).
Lefcheck, J. S. piecewiseSEM: Piecewise structural equation modelling in R for ecology, evolution, and systematics. Methods Ecol. Evol. 7, 573–579 (2016).
Nakagawa, S. et al. Methods for testing publication bias in ecological and evolutionary meta-analyses. Methods Ecol. Evol. 13, 4–21 (2022).
Nakagawa, S., Noble, D. W. A., Lagisz, M., Spake, R., Viechtbauer, W. & Senior, A. M. A robust and readily implementable method for the meta-analysis of response ratios with and without missing standard deviations. Ecol. Lett. 26, 232–244 (2023a).
Nakagawa, S., Yang, Y., Macartney, E. L., Spake, R. & Lagisz, M. Quantitative ev dence synthesis: a practical guide on meta-analysis, meta-regression, and publication bias tests for environmental sciences. Environ.
Evid. 12, 8 (2023b).
Nakagawa, S. et al. orchaRd 2.0: An R package for visualising meta-analyses with orchard plots. Methods Ecol. Evol. 14, 2003–2010 (2023c).
Nakagawa, S. & Santos, E.S.A. Methodological issues and advances in biological meta-analysis. Evol. Ecol. 26, 1253–1274 (2012). R Core Team. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria (2023).
Pinheiro J, Bates D, R Core Team (2023). nlme: Linear and Nonlinear Mixed Effects Models. R package version 3.1-163, https://CRAN.R-project.org/package=nlme.).
R Core Team. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria (2023).
Ruppar, T. Meta-analysis: How to quantify and explain heterogeneity? Eur. J. Cardiovas. Nur. 19, 646–652 (2020).
Smith, S. A. & Brown, J. W. Constructing a broadly inclusive seed plant pylogeny. Am. J. Bot. 105, 302–314 (2018).
The Ministry of Agriculture of the People's Republic of China. Standardized Rice Quality Determination Method from Ministry of Agriculture of the People's Republic of China (NY 147-88). Standards Press of China, Beijing (2010) (in Chinese).
Tanwar, R. K., Garg, D. K., Singh, S. P., Dhandapani, A., & Bambawale, O. M. Enhancement of spiders' population through habitat management in rice (Oryza sativa) fields. Ind. J. Agri. Sci. 81, 462–464
(2011).
van Houwelingen, H. C., Arends, L. R. & Stijnen, T. Advanced methods in meta-analysis: multivariate approach and meta-regression. Statist. Med. 21, 589–624 (2002).
Viechtbauer, W. Conducting Meta-Analyses in R with the metafor Package. J. Stat. Soft. 36, 1–48 (2010).
Lefcheck, J. S. piecewiseSEM: Piecewise structural equation modelling in R for ecology, evolution, and systematics. Methods Ecol. Evol.7, 573–579 (2016).
Wan, N. F. et al. Increasing plant diversity with border crops reduces insecticide use and increases crop yield in urban agriculture. eLife7, e35103 (2018).
Wan, N. F. et al. Global synthesis of effects of plant species diversity on trophic groups and interactions. Nat. Plants 6, 503–510 (2020).
Wan, N. F. et al. Plant genetic diversity affects multiple trophic levels and trophic interactions. Nat. Commu. 13, 7312 (2022a).
Wan, N. F. et al. Spatial aggregation of herbivores and predators enhances tri-trophic cascades in paddy fields: Rice monoculture versus rice-fish co-culture. J. Appl. Ecol. 59, 2036–2045 (2022b
Wang, Y. Q. et al. Understanding biodiversity effects on trophic interactions with a robust approach to path analysis. Cell Rep. Sustain. 2, 100362 (2025).
Wu, J. et al. Comparison of daily predation of rice plathoppers by spiders with their predatory effect in laboratory and caged-rice. Acta Ecol. Sin. 22, 1260–1268 (2002).
Zhang, D. rsq: R-Squared and Related Measures_. R package version 2.5, <https://CRAN.R-project.org/package=rsq> (2022).