An ecological network approach for detecting and validating influential organisms for rice growth

  1. Masayuki Ushio  Is a corresponding author
  2. Hiroki Saito
  3. Motoaki Tojo
  4. Atsushi J Nagano
  1. Hakubi Center, Kyoto University, Japan
  2. Center for Ecological Research, Kyoto University, Japan
  3. Department of Ocean Science, The Hong Kong University of Science and Technology, Clear Water Bay, Kowloon, China
  4. Tropical Agriculture Research Front, Japan International Research Center for Agricultural Sciences, Japan
  5. Graduate School of Agriculture, Osaka Metropolitan University, Japan
  6. Faculty of Agriculture, Ryukoku University, Japan
  7. Institute for Advanced Biosciences, Keio University, Japan

Abstract

How to achieve sustainable food production while reducing environmental impacts is a major concern in agricultural science, and advanced breeding techniques are promising for achieving such goals. However, rice is usually grown under field conditions and influenced by surrounding ecological community members. How ecological communities influence the rice performance in the field has been underexplored despite the potential of ecological communities to establish an environment-friendly agricultural system. In the present study, we demonstrate an ecological-network-based approach to detect potentially influential, previously overlooked organisms for rice (Oryza sativa). First, we established small experimental rice plots, and measured rice growth and monitored ecological community dynamics intensively and extensively using quantitative environmental DNA metabarcoding in 2017 in Japan. We detected more than 1000 species (including microbes and macrobes such as insects) in the rice plots, and nonlinear time series analysis detected 52 potentially influential organisms with lower-level taxonomic information. The results of the time series analysis were validated under field conditions in 2019 by field manipulation experiments. In 2019, we focused on two species, Globisporangium nunn and Chironomus kiiensis, whose abundance was manipulated in artificial rice plots. The responses of rice, namely, the growth rate and gene expression patterns, were measured before and after the manipulation. We confirmed that, especially in the G. nunn-added treatment, rice growth rate and gene expression pattern were changed. In the present study, we demonstrated that intensive monitoring of an agricultural system and the application of nonlinear time series analysis were helpful to identify influential organisms under field conditions. Although the effects of the manipulations were relatively small, the research framework presented here has future potential to harness the ecological complexity and utilize it in agriculture. Our proof-of-concept study would be an important basis for the further development of field-basis system management.

eLife assessment

There is a tremendous need to increase agricultural productivity with means that are both practical and efficient. Drawing on data from variable field environments, this important study provides a theoretical framework for the identification of new factors with presumed relevance for crop growth. This framework can be applied in the context of both agricultural and ecological studies. There is solid evidence for several of the authors' claims, but the impact of the study is limited due to missing functional validation of candidate species in the field. Plant biologists and ecologists working in agricultural and natural environments will find the work interesting.

https://doi.org/10.7554/eLife.87202.3.sa0

Introduction

Global food production supports energy requirements for human beings (Godfray et al., 2010), while the food production is a major driver of greenhouse gas emission and other environmental loads (Aleksandrowicz et al., 2016). How to achieve sustainable food production while reducing environmental impacts is thus a major concern in agricultural science (Lenaerts et al., 2019). Rice (Oryza sativa), one of the world’s major staple crops, is an essential component of the diets and livelihood of over 3.5 billion people (Wing et al., 2018). A large number of studies have investigated how the performance of rice can be improved, and genomics-based breeding, combined with other advanced technologies such as high-throughput field-based phenotyping, is one of the most promising approaches to improve the performance (Wing et al., 2018).

While the advanced breeding techniques are promising, rice is usually grown under field conditions and inevitably influenced by the surrounding biotic and abiotic environment. Previous studies have investigated how meteorological and endogenous (e.g. plant age and genotype) variables influence the gene expression patterns (transcriptome dynamics) of rice under fluctuating field conditions (Kashima et al., 2021; Nagano et al., 2012), and have suggested that transcriptome dynamics can be predominantly explained by endogenous factors and abiotic variables such as ambient air temperature and solar radiation. Nonetheless, biotic variables such as insect herbivory and microbial mutualists/pathogens also play an important role in determining the transcriptome dynamics and productivity of crops (Cohen and Leach, 2019; Savary et al., 2019). However, the dynamics of biotic variables (i.e. ecological community members) under field conditions are difficult to predict because they often show more complex, nonlinear dynamics than abiotic variables (Hsieh et al., 2005). Thus, understanding whether and how biotic variables influence the rice performance (i.e. interspecific interactions) has been underexplored despite its importance for sustainable agriculture (Toju et al., 2018).

In ecological studies, biotic interactions have been studied theoretically and empirically for decades (Allesina and Tang, 2012; May, 1972; Ushio et al., 2018a). Ecologists have traditionally studied interspecific interactions based on observations and manipulations. For example, Paine, 1966 measured the strength of interspecific interactions by removing predator in a rocky tidal system, and showed that the relatively less abundant predator species affected the community diversity disproportionally strongly through interspecific interactions (i.e. such a species is a so-called ‘keystone’ species). Reynolds and Bruno, 2013 also showed that multiple predator species altered population growth of an estuarine food web through complex interactions among species. Thus far, interspecific interactions and their influences on the system dynamics have been studied as in Wootton and Emmerson, 2005. However, despite the substantial contributions to ecology, these observation- and manipulation-based approaches have critical limitations: the identification of multitaxa species and the quantification of their abundance under field conditions are challenging, and the quantification of their interactions is even more difficult (but see Ushio, 2022). Overcoming these difficulties and understanding how ecological community members influence the rice performance under field conditions will provide insights into how we can improve the rice performance and how rice responds to the ongoing and future anthropogenic impacts.

One of the promising approaches for overcoming the current limitations is to monitor the system frequently and detect interspecific interactions using time series data. Recent advances in empirical and statistical methods provide a practical way to achieve this goal. First, use of environmental DNA (eDNA) enables researchers to efficiently detect ecological community members under field conditions (Taberlet et al., 2018). Previous studies have shown that eDNA metabarcoding, an approach to comprehensively amplify and sequence DNAs belonging to a target taxa in an environmental sample, is a cost- and time-effective means to detect a large number of species (e.g. Miya et al., 2015), and the eDNA-based community data is especially informative when it is obtained quantitatively (e.g. sequencing with internal spike-in DNAs: Ushio, 2022; Ushio et al., 2018b). For example, quantitative eDNA metabarcoding enabled effective evaluation of intraspecific genetic diversity (Tsuji et al., 2020) and frequent and comprehensive monitoring of community dynamics (Ushio, 2022; Ushio et al., 2023). Second, nonlinear time series analytical tools enable researchers to reconstruct complex interaction networks (Chang et al., 2021; Deyle et al., 2016; Osada et al., 2023; Schreiber, 2000; Sugihara et al., 2012), and previous studies adopting the approach have detected causality among many variables (e.g., reviewed in Runge et al., 2019). These methods detected and quantified biological interactions in complex systems such as microbiomes (Chang et al., 2021; Fujita et al., 2023), and contributed to understanding and forecasting complex dynamics driven by the interactions.

In the present study, we reconstructed the interaction network surrounding rice and detected potentially influential organisms for rice under field conditions. First, in 2017, we established small experimental rice plots, and monitored rice performance (i.e. growth rates) and ecological community dynamics intensively and extensively. Rice performance was quantified by measuring the growth rates (cm/day), and ecological community members were monitored by a quantitative eDNA metabarcoding approach (Ushio, 2022; Ushio et al., 2018b; Figure 1). We selected growth rates as a monitoring variable for rice because frequent and inexpensive monitoring is possible, and because it integrates various physiological states. The monitoring was performed daily from 23 May 2017 to 22 September 2017 (122 consecutive days). Second, we analyzed the generated extensive time series containing 1197 species and rice growth rates and produced a list of 52 potentially influential species using a time-series-based causality analysis (Figure 2). Third, in 2019, we empirically tested the effects of two species that had been identified as potentially influential in 2017 using manipulative experiments (Figure 3). During the growing season in 2019, an Oomycetes species, Globisporangium nunn (syn. Pythium nunn), was added, and a midge species, Chironomus kiiensis, was removed from small rice plots. The rice responses (the growth rate and gene expression patterns) were measured before and after the manipulation. We confirmed that the two species, especially G. nunn, indeed had statistically clear effects on the rice performance, which demonstrated that integration of the eDNA-based monitoring and time series analysis would be effective for detecting previously overlooked influential organisms in agricultural systems.

Figure 1 with 2 supplements see all
Rice plot, growth rate, air temperature, and ecological community dynamics.

(a) 90 cm × 90 cm rice plot and Wagner pot alignments. Three rice individuals were grown in each pot. Heights and SAPD of the four red individuals in each plot were measured every day during the monitoring period, and the average values of the four individuals were regarded as representative values for each plot. (b) Rice growth rate (cm/day). (c) Daily mean air temperature measured at each rice plot. Upper and lower dotted lines indicate daily maximum and minimum air temperature. (d) Ecological community compositions and average DNA copy numbers per ml water (reported in Ushio, 2022). (e) The number of amplicon sequence variants (ASVs) from each water sample (reported in Ushio, 2022). For (b), (c), and (e), different colors indicate data from different rice plots.

Information transfer between rice growth and ecological community members.

(a) An example of the results of the unified information-theoretic causality (UIC) analysis. Information transfer between air temperature and rice growth rates was quantified. Much higher information transfer was detected from air temperature to rice growth (left panel) compared with the opposite direction (right panel). (b) Strength of causal influence from ecological community members to rice growth. Transfer entropy (TE) quantified by the UIC method was used as an index of causal influence. Colors indicate taxa assigned to ASVs. y-axis indicates ASV ID. Note that the prefix (e.g. ‘Fungi_’) of the IDs corresponds a major target group of the primer and does not necessarily indicate a taxonomic group assigned to the ASV (see Supplementary file 1). (c) eDNA dynamics of putative Globisporangium nunn (Fungi_Taxa00402 in Supplementary file 1). (d) eDNA dynamics of Chironomus kiiensis (total DNA copy numbers of five midge ASVs). For (c) and (d), different colors indicate data from different rice plots.

Figure 3 with 2 supplements see all
The manipulation experiment performed in 2019 and ecological community compositions before and after the manipulation.

(a) Setting of the manipulation experiment in 2019. The number and characters next to each plot indicate the plot number and treatment. The left-top inset shows daily mean air temperature (thick line) and daily maximum and minimum air temperature (dashed lines). The right-bottom inset shows three individuals (red and green points) in each Wagner pot, and the number in each red individual indicates the pot location number. Heights and SPAD of the red individuals were measured. Total eDNA copy numbers of (b) putative Globisporangium spp. and (c) midge (Chironomus kiiensis) in the rice plots. ‘before’ and ‘after’ indicate ‘from 18 June to 24 June’ and ‘from 25 June to 12 July’, respectively. (d) Overall community compositions after the manipulation. Gray, red, and blue indicate CT (control), GN (Globisporangium nunn added), and CK (Chironomus kiiensis removed) treatments, respectively.

Results

Field monitoring of rice growth and ecological communities in 2017

In five rice plots established in 2017 in an experimental field at Kyoto University, Japan, daily rice growth rate (cm/day in height) was monitored during the growing season by measuring rice leaf height of target individuals (red points in the right panel of Figure 1a) every day using a ruler, which showed consistent patterns among the plots (Figure 1b). Daily mean air temperature also showed consistent patterns among the five plots (Figure 1c). The daily growth rate reached the maximum during late June to early July, and the height of rice individuals did not increase after the middle of August (first headings appeared on 12 or 13 August in the five plots). During the monitoring period, we occasionally observed decreases in the rice heights due to mechanical damage or insect herbivores, but it is unlikely that they affected our causal inferences because the changes in the rice heights due to the damages and herbivories were smaller and less frequent than those due to growth. See Figure 1—figure supplement 1 and a video (https://doi.org/10.6084/m9.figshare.19029650.v1) for additional information regarding the monitoring.

Daily ecological community dynamics was monitored using quantitative eDNA metabarcoding by analyzing eDNA in water samples taken from the five plots with four universal primer sets (16 S rRNA, 18 S rRNA, ITS, and COI regions targeting prokaryotes, eukaryotes, fungi, and animals, respectively; see Materials and methods), and detailed patterns of the ecological community dynamics and their implications were reported in Figure 1—figure supplement 2, a video (https://doi.org/10.6084/m9.figshare.23514150.v1), and Ushio, 2022. Briefly, the total eDNA copy number increased late in the sampling period (Figure 1d). In contrast, ASV diversity (a surrogate of species diversity) was highest in August and then decreased in September (Figure 1e). Prokaryotes largely accounted for this pattern (Figure 1d). In the previous study, a large ecological interaction network was reconstructed (Figure 1—figure supplement 1c), and possible mechanisms of the ecological dynamics were discussed in detail (Ushio, 2022). Importantly, the time series-based causality analysis used in the next section requires quantitative time series, and our quantitative eDNA time series are suitable for this purpose.

Detection of potential causal relationships between rice growth and ecological community dynamics

Potential causal relationships among the daily rice growth rates, eDNA time series, and climate data were detected using unified information-theoretic causality (UIC) analysis, which quantifies information flow in terms of transfer entropy (TE; Figure 2; for UIC, see ‘Detection of potential causal relationships between rice growth and ecological community members’ in Materials and methods and Osada et al., 2023). The information flow between eDNA time series and rice growth rates was regarded as a sign of interactions between rice and the ecological community member. Figure 2a is a demonstration of the performance of the UIC analysis, showing the influence of air temperature on the rice growth. The left panel in Figure 2a indicates that there were strong effects of air temperature on the rice growth (x-axis indicates the time-lag of the effects), while the right panel indicates that there were virtually no effects of the rice growth on air temperature, which is biologically convincing. In total, 718 ASVs were identified as potentially causal (i.e. statistically clear information flow; p<0.05; we use the term ‘statistical clarity’ instead of ‘statistical significance’, according to the recommendation by Dushoff et al., 2019), among which 52 ASVs were identified with lower-level taxonomic information (Supplementary file 1). The potentially causal ASVs belong to various taxa, including fungi and insects, with varying degree of the effects, measured as the information flow (Figure 2b). Except for one species (Penicillium citrinum), we could not find information regarding the effects of the detected ASVs on rice performance, at least at the time of this analysis.

Among the 52 ASVs, we particularly focused on two ASVs, Globisporangium nunn (Oomycete; previously known as Pythium nunn) (Uzuhashi et al., 2010) and Chironomus kiiensis (midge species; see the eDNA dynamics of putative G. nunn and C. kiiensis in Figure 2c and d), in the field manipulation experiments in 2019. The reasons why we focused on these two species in 2019 are multifold. First, they can be relatively easily manipulated; G. nunn had already been isolated and cultivated (Kobayashi et al., 2010; Tojo et al., 1993) and C. kiiensis is an insect species and relatively large (several mm to cm). Second, the two species had not been reported to be pathogens; using pathogen species under field conditions is not recommended in the experimental field because they might have adverse effects on nearby agricultural systems. G. nunn is a potential biocontrol agent and is reported to be beneficial for vegetables (Kobayashi et al., 2010). C. kiiensis is a common midge species around the study region (including in rice paddy fields). The larvae are usually about 10 mm and adults are about 5 mm. Additional information for the two species is described in Materials and methods and Supplementary file 2 (see also ‘Notes on Globisporangium nunn’).

The ecological community composition after the manipulation experiments in 2019

In 2019, we established nine experimental plots to monitor rice growth and to perform field manipulation experiments in the same experimental field. The plots were similar to, but smaller than, those used in 2017 so that the number of replicates could be increased (Figure 3a). We prepared three treatments, namely, the control treatment (CT), G. nunn-added treatment (GN), and C. kiiensis-removed treatment (CK), and there were three replicates for each treatment (3 treatments ×3 replicates = 9 plots). We manipulated the abundance of G. nunn and C. kiiensis by adding G. nunn and removing C. kiiensis, respectively (Figure 3—figure supplement 1), and the manipulation experiments were performed three times (on 24, 26, and 28 June 2019). We quantified the eDNA dynamics of the ecological communities weekly or biweekly except during the intensive, daily monitoring period before and after the field manipulation experiments (from 18 June to 12 July; see the sampling design in Supplementary file 3 and overall eDNA dynamics in Figure 3—figure supplement 2). As for the GN treatment, we detected 19 putative Globisporangium sequences belonging to the family Pythiaceae from the four marker genes (i.e. 16 S, 18 S, ITS and COI), most of which we could not identify with lower-level taxonomic information partially because of the conservative assignment algorithm of the sequence analysis pipeline, Claident (Tanabe and Toju, 2013). We described them as ‘putative Globisporangium spp.’ and regarded them as G. nunn that we added because the eDNA concentrations of these Pythiaceae sequences were increased after we added G. nunn isolates.

The increase in the eDNA concentrations of putative Globisporangium spp. in the GN treatment after the manipulation was statistically clear (p<2.0 × 10–16; Figure 3b). In addition, there was a slight, but statistically clear increase in the eDNA concentrations of putative Globisporangium spp. in the CT treatment after the manipulation (p=0.032; Figure 3b). On the other hand, the eDNA concentrations of putative Globisporangium spp. in the CK treatment was not changed after the manipulation (p>0.05). The eDNA concentrations of C. kiiensis were statistically clearly changed after the manipulation in the CT (slightly decreased; p=0.045; Figure 3c) and GN (increased; p=6.6 × 10–6; Figure 3c) treatments. In the CK treatment, the eDNA concentration of C. kiiensis was not changed (p>0.05; Figure 3c). These results were further visualized and confirmed by the t-distributed stochastic neighbor embedding (t-SNE) (Van der Maaten and Hinton, 2008) analysis of the overall community composition. Ecological communities of the GN treatment were generally separated from those of the CT treatment, while those of the CK treatment almost overlapped with those of the CT treatment (Figure 3d). Alternative statistical modeling that included the treatments (the control versus GN or CK treatments) and manipulation timing (i.e. before or after the manipulation), which simultaneously took the temporal changes of all the treatments into account, also showed qualitatively similar results (Supplementary file 4), further supporting the results.

Rice performance after the field manipulation experiments in 2019

We measured the rice growth as in 2017 and compared those of the GN and CK treatments with that of the CT treatment before and after the manipulation experiments. In general, we did not find statistically clear differences among the treatment in the daily rice growth rate before and after manipulation (p>0.05; Figure 4a; see Figure 4—figure supplement 1 for the rice growth trajectory during the growing season). However, we found that the cumulative rice growth (during 10 days after the manipulation) increased in the GN and CK treatments after the manipulation, and the effect was statistically clear (p=0.016 and p=0.021 for the GN and CK treatments, respectively; Figure 4b). Alternative statistical modeling that included the treatments (the control versus GN or CK treatments) and manipulation timing (i.e. before or after the manipulation), which simultaneously took the temporal changes of all the treatments into account, also showed qualitatively similar results (Materials and methods and Supplementary file 5). We did not find any statistically clear difference in the rice yields (e.g. the number of rice grains per head, the grain weights, or the proportion of fertile grains) among the treatments (see Supplementary file 6 for the raw rice yield data).

Figure 4 with 1 supplement see all
Rice growth rate and cumulative growth before and after the manipulation experiment in 2019.

(a) Growth rates and (b) cumulative growth of the rice individuals in the three treatments (CT = control; GN = Globisporangium nunn added; CK = Chironomus kiiensis removed) before and after the manipulation (cumulative growth was calculated by summing up height growth before the third manipulation or during 10 days after the third manipulation).

Changes in rice gene expression pattern in 2019

We further explored changes in the rice performance before and after the manipulation experiments using RNA-seq. We collected rice leaf samples four times (1 day before the first manipulation and 1, 14, and 38 days after the third manipulation). As a result of the RNA-seq analysis, 875,534,703 reads were generated (Supplementary file 7), and 8,105,902 reads per sample were used for the sequence analysis. DESeq2 analysis was performed, and here we focus on the results of the second sampling event (i.e. 1 day after the third manipulation). The analysis suggested that, when all three rice pot locations within each plot were merged (see Figure 3a for the pot locations), there were differentially expressed genes (DEGs) in the GN treatment (Figure 5a; see Supplementary file 7 for nuclear DEGs) while there were no differentially expressed genes in the CK treatment (Figure 5b). However, when the data was separately analyzed for each Wagner pot location, we found that there were several differentially expressed genes even in the CK treatment (Figure 5—figure supplement 1 and Supplementary file 8), suggesting that there was a pot location-specific effect and that the CK treatment also had small effects on the rice gene expression. We found almost no DEGs for leaf samples taken one day before and 14 and 38 days after the third manipulation (the leaf sampling event 1, 3, and 4), suggesting that the influences of the treatments on the gene expression patterns were transient. In Figure 6, we show examples of differentially expressed genes. In the GN treatment, the expressions of three genes, Os12g0504050, Os11g0184900, and Os01g0678500 (p<0.0001, negative binomial GLMM; Figure 6a–c), were increased, while those of three other genes (Os01g0642200, Os08g0162800, and Os03g0285700) were decreased (p<0.0001, negative binomial GLMM; Figure 6d–f). These DEGs are mostly related to the photosynthetic system and stress responses (Supplementary file 7).

Figure 5 with 1 supplement see all
Differential expression genes analysis.

(a)Globisporangium nunn-added and (b)Chironomus kiiensis-removed treatment. Red and blue points indicate significant up- and down-regulated genes, respectively. Upper and lower dashed lines indicate lo2(1.5) and −log2(1.5), respectively.

Examples of differentially expressed genes after the manipulation experiment.

Results of (a) Os12g0504050, (b) Os01g0642200, (c) Os08g0162800, (d) Os11g0184900, (e) Os03g0285700, and (f) Os01g0678500 are presented. y-axis represents DESeq2-normalized read counts. Gray, red, and blue indicate CT (control), GN (Globisporangium nunn added), and CK (Chironomus kiiensis removed) treatments, respectively. The gene expressions of the GN treatment in all six genes are statistically clearly different from those of the other two treatments (p<0.0001) except for the panel GN vs. CK in c (p=0.00014) and GN vs. CK in (d) (p=0.0087).

Discussion

In the present study, we demonstrated a novel research framework to integrate the ecological network concept and agricultural practices to screen potentially influential organisms. We performed intensive field monitoring in 2017, obtained intensive and extensive time series of rice and ecological communities, detected potentially influential organisms using time series-based causality analysis, and based on the results, we performed field manipulation experiments and evaluated the effect of the manipulation. Throughout the study, we detected statistically clear effects of the manipulation on the rice growth and gene expression patterns especially for the GN treatment, proving that our proposed framework might work as a novel approach to comprehensively detect potentially influential, previously overlooked ecological community members for an agricultural system (in our case, rice). Below, we discuss the results, advantages, disadvantages, potential limitations, and future perspectives of our approach.

Monitoring of ecological community using eDNA and rice growth

Our framework started from monitoring ecological community members using quantitative eDNA metabarcoding (Ushio, 2022). We successfully detected more than 1000 ASVs and described detailed temporal dynamics of the ASVs (Figure 1d and e; Ushio, 2022). We included spike-in standard DNAs whose concentrations were known a priori, and our previous study had shown that, in almost all cases, there were clear positive linear correlations between the sequence reads and the copy numbers of the spike-in standard DNAs (e.g. Supplementary figures in Ushio, 2022). This suggested that our quantitative eDNA metabarcoding reasonably quantified the concentrations of eDNA, which can be a ‘rough’ index of species abundance in the artificial rice plots. Some time-series-based causal inference methods required quantitative data (i.e. relative abundance data that is typical for amplicon sequencing is not suitable; Osada et al., 2023; Schreiber, 2000; Sugihara et al., 2012), and our method justifies the use of the advanced time series method to infer the causal relationship between rice growth and ecological community members.

Rice growth was monitored by a manual measurement, for example, rice height measurements made using rulers (Figure 1b). Although the measurements successfully captured the seasonal dynamics of the rice growth, one might speculate whether only a single variable is sufficient to capture the rice performance. Indeed, it is true that assessing more variables, for example, gene expression (Nagano et al., 2019; Nagano et al., 2012) and multispectral image monitoring during the whole growth season (Candiago et al., 2015), would be potentially useful to fully capture the trajectory of the rice growth. Nonetheless, we suggest that very frequent monitoring of a single variable is sufficient to reconstruct and predict the dynamics of the whole system dynamics (Takens, 1981). Thus, we suggest that our monitoring approach can resolve the dynamics at least as a proof-of-concept study, but taking more variables at fine time resolution is certainly beneficial for more accurate delineation of the system dynamics.

Detection of potentially influential organisms using nonlinear time series analysis

By using UIC (Osada et al., 2023), we detected 718 potentially influential ASVs. Short read amplicon sequencing and conservative taxa assignment algorithm of Claident (Tanabe and Toju, 2013) did not allow to assign the species name to most of the ASVs, and we could assign the species name only for 52 ASVs. Among them, we detected a previously known detrimental fungus, Penicillium citrinum, that may colonize rice and produce a toxin called citrinin that could cause human health problems (Ali, 2018; Ostry et al., 2013), which partially supported the validity of our time series analysis. However, at least at the time of the analysis, it is unknown whether the species in Supplementary file 1 really include influential organisms.

In 2019, we focused on two species, G. nunn and C. kiiensis (Figure 3—figure supplement 1). At the time of that analysis, these choices were reasonable: G. nunn was previously reported to have some antagonistic activity against pathogens to vegetables such as Pythium ultimum (Kobayashi et al., 2010) and did not have any clear detrimental effects on rice (which was necessary in order not to spread potential pathogens to nearby farms), and C. kiiensis (midge) is a common freshwater insect species in the study region. C. kiiensis may have effects on water quality (e.g. the concentration of phosphorus) through their feeding behavior (Kawai et al., 2003). However, they have never been focused on in terms of rice growth management, and thus, they provided an opportunity to test the validity of our approach in the present study.

Field manipulation experiments in 2019

In 2019, we performed the manipulation experiments. Although the CK treatment did not change the community composition statistically clearly, the GN treatment changed the community compositions (Figure 3b–d). In the GN treatment, the abundance of putative Globisporangium spp. clearly increased (Figure 3b), suggesting that the introduction of G. nunn was successful. However, we also observed increases in the putative Globisporangium spp. abundance in some of the CT and CK samples. The rice plots were located close to each other (50–60 cm apart; Figure 3a) and thus immigrations might have occurred between the GN and the other plots. As for C. kiiensis, we did not observe a statistically clear decrease in C. kiiensis eDNA in the CK treatments (Figure 3c), suggesting that C. kiiensis individuals or their eDNA persisted even after the removal treatment. This might suggest that the treatment was not successful, or that quantitative eDNA metabarcoding failed to capture the actual decrease in their abundance. On the other hand, we observed statistically clear increases in C. kiiensis eDNA in the CT and GN treatments. This response was unexpected, but the effects of the increase in G. nunn might be transmitted to other species through the interaction networks which were reconstructed in the previous study (Ushio, 2022).

As for the rice growth, although the manipulations were mild (the GN treatment) or might have been relatively unsuccessful (the CK treatment) compared with other conventional manipulations such as the addition of fertilizers and insecticides, we surprisingly found statistically clear effects of the manipulations under field conditions (Figure 4b). The addition of G. nunn increased the cumulative rice growth and changed the gene expression pattern after the manipulations statistically clearly. The removal of C. kiiensis also affected the rice growth (Figure 4b), and the gene expression patterns were changed when the pot locations were taken into account (Figure 5—figure supplement 1c and f and Supplementary file 8), suggesting that the effects of ecological community members might be situation dependent (e.g., microclimate changes depending on the pot locations). The unclear effects of the CK treatment relative to those of the GN treatment could be due to the relatively unstable removal method (i.e. C. kiiensis larvae were manually removed by a hand net) or incomplete removal of the larvae (some larvae might have remained after the removal treatment). Unfortunately, detailed mechanisms of these effects cannot be elucidated in the present study, but some inferences could be made from the gene expression patterns. For example, nuclear DEGs included genes related to photosynthesis and stress responses (Supplementary files 7 and 8). These changes were observed relatively quickly and transient (only found in leaf samples collected one day after the third manipulation; Supplementary file 3) and thus they could be direct effects from G. nunn. Up-regulations of photosynthesis-related genes (e.g. Os08g0433350, Os04g0473150, Os12g0504050, and Os01g0791033 in Supplementary file 7) might have contributed to the increase in the cumulative rice heights. Stress response genes such as Os11g0184900 might also have contributed to the changes in the rice performance, but further studies are required to elucidate the mechanisms by which ecological community members affect the rice performance under field conditions.

In the present study, we could not test the effect of other species on rice growth due to the limitations of time, labor, and money, but the species listed in Supplementary file 1 could be candidates of novel biological agents that influence the rice performance. Thus, validating the effects of these species on the rice performance is important, and targeting ASVs with a large TE value would be a promising direction to find previously overlooked beneficial/detrimental organisms for rice under field conditions. In addition, we could not test the effects of species that had no clear effects on the rice performance according to the time series analysis. It will also be important to validate that such species do not have clear influence on the rice performance in the field, which would further validate the effectiveness of our research framework.

Potential limitations and future perspectives

Although the present study has shown that our framework is potentially useful for detecting potentially influential, previously overlooked organisms for rice growth under field conditions, there are potential limitations that should be acknowledged.

One limitation is the quantitative capacity and species identification accuracy of the eDNA metabarcoding. As for the first issue, the quantitative eDNA metabarcoding approach used in this study provides reasonable estimates of eDNA concentrations (see the following studies for the comparison between quantitative eDNA metabarcoding and the fluorometric method and quantitative PCR; Ushio, 2022; Ushio et al., 2018b), and the eDNA concentrations usually positively correlate with the number or biomass of individuals. However, there are species- and environment-specific eDNA dynamics such as degradation and release processes (Dejean et al., 2011; Maruyama et al., 2014), which may bias the estimation of absolute abundance. Although our causality detection might not be strongly influenced by the species-specific bias (i.e. time series were standardized before the analysis), accurate estimations of the absolute abundance will improve the accuracy of analyses. Identifications and enumerations of individuals using eDNA would be a fascinating direction to improve the accuracy of the species abundance estimations. As for the second issue, short-read sequencing has dominated current eDNA studies, but it is often not sufficient for lower-level taxonomic identification (i.e. species or subspecies identification). Using long-read sequencing techniques (e.g. Oxford Nanopore MinION) for eDNA studies is a promising approach to overcome the second issue (Baloğlu et al., 2021).

Another limitation is that the monitoring of rice and the introduced organisms in this study was not comprehensive. While intensive monitoring was conducted during the rice growing season, only a limited number of variables were measured for a limited number of rice individuals and the detailed dynamics of the two introduced species was unclear (i.e. the fate of the introduced species). This is particularly important for understanding how the introduced organisms affected rice performance. At the molecular level, acquiring larger scale ‘omics’ data for rice and ecological community members, for example long-term gene expression monitoring data (Nagano et al., 2019), will facilitate understanding of detailed mechanisms. In addition, at the field level, future studies should consider applying a more advanced monitoring approach, such as unamended aerial vehicle (Fujiwara et al., 2022), multispectral camera (Candiago et al., 2015), and/or neural network-based image analysis (Toda and Okura, 2019).

Lastly, while we used an advanced nonlinear time series analysis to detect potential causal organisms for rice growth (Osada et al., 2023), there is still room for improvement. Nonlinear time series analysis is also a rapidly developing field, and many new techniques are emerging for causal inference (Runge et al., 2019; Suzuki et al., 2022), quantification of interaction strength (Chang et al., 2021) and near future forecasting (Chen et al., 2020). Applying these techniques to agricultural data may help to identify previously overlooked beneficial organisms under field conditions.

Conclusions

In this study, we demonstrated that intensive monitoring of an agricultural system and the application of nonlinear time series analysis are helpful for identifying potentially influential organisms under field conditions. Although the effects of the two species were relatively small and how the two species influenced the rice performance (i.e. the fate of the introduced species) was still unclear, the research framework presented here has future potential. The intensive monitoring, nonlinear time series analysis, and in situ validation of the statistical results may be applicable for other agricultural and aquatic systems. For example, the application of our framework to fishery aquaculture systems may be helpful for detecting potential pathogens of fish. Also, if the continuous monitoring of soil ecological communities becomes possible, our approach can be extended to other crops and vegetables. An automated, real-time monitoring system combined with advanced machine learning methods would provide a promising tool for increasing the accuracy and applicability of our approach. In conclusion, we proposed a novel framework to integrate the ecological network concept and agricultural practices to explore and detect potentially influential organisms. Our proof-of-concept study would be an important basis for the further development of field-based system management.

Materials and methods

Field experimental setting and rice monitoring in 2017

Request a detailed protocol

Five artificial rice plots were established using small plastic containers (90×90 × 34.5 cm; 216 L total volume; Risu Kogyo, Kagamigahara, Japan) in an experimental field at the Center for Ecological Research, Kyoto University, in Otsu, Japan (34° 58′ 18″ N, 135° 57′ 33″ E) (Figure 1a and Figure 1—figure supplement 1a and b). Sixteen Wagner pots (φ174.6–160.4 mm [top–bottom] × 197.5 mm; AsOne, Osaka, Japan) were filled with commercial soil, and three rice seedlings (var. Hinohikari) were planted in each pot on 23 May 2017 and then harvested on 22 September 2017 (122 days). The containers (‘plots’) were filled with well water. This experimental system was previously reported in Ushio, 2022, which investigated the mechanism of ecological community dynamics. In the present study, we focused on the relationship between rice performance and ecological community members.

Daily rice growth was monitored by measuring rice leaf height of target individuals every day using a ruler (the highest leaf heights were measured). We selected four rice individuals at the center of each plot as the target individuals (indicated by the four red points in the right panel of Figure 1a). The average height of the four individuals were used as a representative rice height of each plot (a time-lapse movie that shows rice growth is available here: https://doi.org/10.6084/m9.figshare.19029650.v1). Leaf SPAD was also measured every day using a SPAD meter (SPAD-502Plus, KONICA-MINOLTA, Inc, Tokyo, Japan) for the same leaf whose height was measured, but not analyzed intensively in the present study. Climate variables (temperature, light intensity, and humidity) were monitored using a portable, automated climate logger (Logbee, Chitose Industries, Osaka, Japan). The Logbee data was corrected and refined using the climate station data provided by the Center for Ecological Research, Kyoto University.

Water sampling for the eDNA-based monitoring of ecological communities in 2017

Request a detailed protocol

The complete information about eDNA-based ecological community monitoring is available in Ushio, 2022 (https://ndownloader.figstatic.com/files/34067324 and https://ndownloader.figstatic.com/files/34067327). Approximately 200 ml of water in each rice plot was collected from each of the four corners of the plot using a 500 ml plastic bottle and taken to the laboratory within 30 min. The water was filtered using Sterivex filter cartridges (Merck Millipore, Darmstadt, Germany). Two types of filter cartridges were used to filter water samples: to detect microorganisms, φ0.22-µm Sterivex (SVGV010RS) filter cartridges that included zirconia beads inside were used (Ushio, 2019), and to detect macroorganisms, φ0.45-µm Sterivex (SVHV010RS) filter cartridges were used. After filtration, 2 ml of RNAlater solution (ThermoFisher Scientific, Waltham, Massachusetts, USA) were added to each filter cartridge to prevent DNA degradation during storage. In total, 1220 water samples (122 days × 2 filter types × 5 plots) were collected during the census term. In addition, 30 field-level negative controls, 32 PCR-level negative controls with or without the internal standard DNAs, and 10 positive controls to monitor the potential DNA cross-contamination and degradation during the sample storage, transport, DNA extraction and library preparations were used.

Quantitative analysis of environmental DNA: Library preparations, sequencing, and sequence data processing

Request a detailed protocol

Detailed protocols of the quantitative eDNA metabarcoding are described in Ushio, 2022, and the protocols below were also adopted for the samples collected in 2019. Briefly, DNA was extracted and purified using a DNeasy Blood & Tissue kit (Qiagen, Hilden, Germany). After the purification, DNA was eluted using 100 µl of the elution buffer and stored at −20 °C until further processing.

A two-step PCR approach was adopted for the library preparation for quantitative MiSeq sequencing. The first-round PCR (first PCR) was carried out with the internal standard DNAs to amplify metabarcoding regions using primers specific to prokaryotes (515 F and 806 R; Bates et al., 2011; Caporaso et al., 2011), eukaryotes (Euk_1391 f and EukBr; Stoeck et al., 2010), fungi (ITS1-F-KYO1 and ITS2-KYO2; Toju et al., 2012) and animals (mlCOIintF and HCO2198’ Folmer et al., 1994; Leray et al., 2013). The second-round PCR (second PCR) was carried out to append indices for different samples for sequencing with MiSeq. The DNA library was sequenced on the MiSeq (Illumina, San Diego, CA, USA).

The raw MiSeq data were converted into FASTQ files using the bcl2fastq program provided by Illumina (bcl2fastq v2.18). The FASTQ files were then demultiplexed using the command implemented in Claident v0.2.2019.05.10 (http://www.claident.org; Tanabe and Toju, 2013). Demultiplexed FASTQ files were then analyzed using the Amplicon Sequence Variant (ASV) method implemented in ‘dada2’ package (Callahan et al., 2016) of R (R Development Core Team, 2022). Taxonomic identification was performed using Claident (Tanabe and Toju, 2013). We detected negligible sequences from the negative control samples.

Estimations of DNA copy numbers and validation of the quantitative capability of the MiSeq sequencing with internal standard DNAs

Request a detailed protocol

For all analyses after this subsection, the free statistical environment R was used (R Development Core Team, 2022). The procedure used to estimate DNA copy numbers consisted of two parts, following previous studies (Ushio, 2019; Ushio et al., 2018b). Briefly, we performed (i) linear regression analysis to examine the relationship between sequence reads and the copy numbers of the internal standard DNAs for each sample, and (ii) conversion of sequence reads of non-standard DNAs to estimate the copy numbers using the result of the linear regression for each sample. The regression equation was: MiSeq sequence reads = sample-specific regression slope ×the number of standard DNA copies [/µl]. Then, the estimated copy numbers per µl extracted DNA (copies/µl) were converted to DNA copy numbers per ml water in the rice plot (copies/ml water). The quantitative capacity of this method was thoroughly evaluated by comparing with quantitative PCR, fluorescence-based DNA measurement, and shotgun metagenomic analysis (Ushio, 2022; Ushio et al., 2018b), and the method was shown to have reasonable capacity to quantify DNA. The method has already been shown to be effective for eDNA-based quantitative assessment of ecological community (Sato et al., 2021; Tsuji et al., 2020). The estimation of eDNA concentrations was necessary to apply the time series-based causality analysis explained in the next section.

Detection of potential causal relationships between rice growth and ecological community members

Request a detailed protocol

We detected potentially influential organisms for rice growth analyzing the quantitative, 1197-species eDNA time-series (Ushio, 2022) and the daily rice growth rate obtained in 2017 using nonlinear time series analysis. We quantified information flow from eDNA time series to rice growth rate (i.e. a proxy of interaction between an organism and rice growth) by the ‘unified information-theoretic causality (UIC)’ method (Osada et al., 2023) implemented in the ‘rUIC’ package (Osada and Ushio, 2021) of R. UIC tests the statistical clarity of information flow between variables in terms of transfer entropy (TE; Frenzel and Pompe, 2007; Schreiber, 2000) computed by nearest neighbor regression based on time-embedded explanatory variables (i.e. cross-mapping: Sugihara et al., 2012). In contrast to the standard method used to measure TE, UIC quantifies information flow as follows:

(1) TE=1Tt=1Tlog(p(yt+tpxt,xtτ,,xt(E1)τ,zt)p(yt+tpxtτ,xt2τ,,xt(E1)τ,zt)),

where x, y, and z represent a potential causal variable, effect variable, and conditional variable (in our case, z is/are climate variable(s) if they are causal factors), respectively. p(AB) represents conditional probability: the probability of A conditioned on B. t, tp, τ, and E represent the time index, time step, a unit of time lag, and the optimal embedding dimension, respectively. T is the number of times to evaluate TE. For example, if tp = –1 in Eqn. (1), UIC tests the causal effect from yt–1 to xt. Optimal E was selected by measuring TE as follows:

(2) TE=1Tt=1Tlog(p(xt+tpyt,xt,xtτ,,xt(E1)τ,zt)p(xt+tpyt,xt,xtτ,,xt(ER1)τ,zt)),

where ERE is the optimal embedding dimension of lower dimensional models. Eqn. (2) is a TE version of simplex projection (Sugihara and May, 1990), and we set tp = 1 in Eqn. (2) to determines the optimal E (i.e. the optimal E based on one-step forward prediction). Statistical clarity was tested by bootstrapping data after embedding (a threshold was set to 0.05). TE measured according to Eqn. (1) gains the advantage of previous causality tests, that is, standard TE methods (Runge et al., 2012; Schreiber, 2000) and convergent cross mapping (CCM) (Sugihara et al., 2012), as explained in Osada et al., 2023.

By using UIC, we quantified TE from eDNA time series to rice growth. We standardized the eDNA time series (copies/ml water) and rice growth rates (cm/day) to have zero means and a unit of variance before the analysis. We tested time-lag up to 14 days (i.e. tp in Eqn. (1) was from 0 to –14; effects from up to 14 days ago were considered). We interpret a finding that TE from eDNA time series to rice growth rate is statistically greater than zero to mean that the ecological community member statistically clearly influences rice growth rate.

Field experimental setting in 2019

Request a detailed protocol

Based on the results in 2017, we performed a field manipulative experiment. We focused on two potentially influential organisms as explained in the next section. Nine artificial rice plots were established using smaller plastic containers than those used in 2017 (42.2 × 32.0 × 30.0 cm; 40.5 L total volume; Sanko, Tokyo, Japan) in the same experimental field at Kyoto University (i.e. the identical location with Plot 3 in 2017) (Figure 3a). Three Wagner pots were filled with commercial soil, and three rice seedlings (var. Hinohikari) were planted in each pot on 20 May 2019. Three treatments were prepared (see the next section), including the control treatment, and there were three replicates for each treatment (3 treatments × 3 replicates × 3 Wanger pots × 3 rice seedlings = 81 rice seedlings in nine rice plots). The rice individuals were harvested on 20 September 2019 (124 day growing period).

Field manipulation experiments in 2019

Request a detailed protocol

Field manipulation experiments were performed using two potentially influential species, Globisporangium nunn (Oomycete; previously known as Pythium nunn; Uzuhashi et al., 2010) and Chironomus kiiensis (midge species; Figure 3—figure supplement 1). The reasons why we focused on these two species are described in the Results. As for G. nunn, we referred to the ASV found by the eDNA metabarcoding in 2017 as ‘putative Globisporangium’ because there are differences in the ASV sequence and the G. nunn type strain sequence. Detailed explanations about the difference in the sequences are explained in ‘Notes on Globiosporangium nunn’ in Supplementary file 2.

We performed the rice plot manipulations three times at around noon on 24, 26, and 28 June 2019. There were three treatments, namely, ‘Globisporangium nunn-added’ (GN), ‘C. kiiensis-removed’ (CK), and control (CT) treatments with three replicates for each treatment. For the GN treatment, isolated and incubated G. nunn was mixed in vermiculite (Kobayashi et al., 2010; Tojo et al., 1993) and approximately 200 ml of vermiculite was added to each Wagner pot (Figure 3—figure supplement 1a and b). Vermiculite without G. nunn was added to Wagner pots in the two other treatments. For the CK treatment, midge larvae were removed instead of adding them because isolation and cultivation of this insect species are technically difficult and time- and effort-consuming. Midge larvae were removed using a φ1.0 mm net (Figure 3—figure supplement 1c). On average, 283 midge larvae were removed from each plot and the DNA of about 100 removed larvae was sequenced using Sanger sequencing, and all were confirmed as C. kiiensis. For the other treatments, the net was also inserted into plots, the water in each plot was mixed by the net for a while, but midge larvae were not removed.

Rice growth and ecological community monitoring, and rice leaf sampling during the manipulation experiments

Request a detailed protocol

Rice growth was monitored by measuring rice leaf height of target individuals using a ruler, as in 2017. We selected one rice individual in the middle of each of the Wagner pots as the target individuals (indicated by the three red points in the inlet panel of Figure 3a; labelled as locations ‘1’, ‘2’, and ‘3’ from east to west). Leaf SPAD was also measured using a SPAD meter as in 2017. The ecological communities in the rice plots were monitored by analyzing eDNA as in 2017. After we harvested rice, we counted the number of fertile and sterile grains and quantified rice yields.

In 2019, the rice growth was monitored weekly or biweekly except during the period of the manipulation experiments. The manipulations were performed three times (i.e. 24, 26, and 28 June 2019), and the daily monitoring of the rice plots was performed before and after the manipulations (Supplementary file 3). Rice performance (growth rates and SPAD) was monitored every day from 11 June to 12 July (32 days). Water samples for eDNA-based ecological community monitoring were collected every day from 18 June to 12 July (25 days). Rice leaf samples for RNA-seq were collected four times before and after the field manipulations (Supplementary file 3). The leaf samples were immediately frozen by dipping into liquid nitrogen under field conditions, and kept frozen until stored in a freezer at –20 °C. In total, 108 leaf samples were collected for RNA-seq analysis (Supplementary file 7).

Effects of the field manipulation experiment on rice growth and ecological communities

Request a detailed protocol

The differences in rice performance among the three treatments and total eDNA concentrations of ecological community members (copies/ml water) before and after the manipulation experiments were tested using linear mixed model (LMM) and general linear mixed model (GLMM), respectively, using ‘lme4’ package of R (Bates et al., 2015). For rice growth rates, the effect of the treatment was tested with the random effects of rice individual and rice plot (in R, this is lme4::lmer(growth_rate ~treatment + (1|ind/plot), data = data)), and the effects were separately analyzed for before/after the manipulation. Rice cumulative growth (i.e. cumulative growth before the third manipulation or 10 days after the third manipulation) was analyzed similarly, but the random effect of rice individuals was not included because there was only one cumulative value for each rice individual. Also, we tested alternative models including the timing and manipulation treatments (in R, these are lme4::lmer(growth_rate ~before_or_after_manipulation*treatment + (1|ind/plot), data = data) and lme4::lmer(cumulative_growth ~before_or_after_manipulation*treatment + (1|plot), data = data)), and obtained general agreement with the results of the above-mentioned analysis. For ecological community, we analyzed the eDNA data taken during the intensive monitoring period before and after the field manipulation experiments (18 June to 12 July; Supplementary file 3), and the effects of the manipulation (i.e. the addition of G. nunn or the removal of C. kiiensis) were separately tested for each treatment with the random effect of rice plot assuming a gamma error distribution (in R, this is lme4::glmer(DNA_conc ~ before_or_after_manipulation + (1|plot), data = data, family = Gamma(link="log"))). Also, as in the rice growth analysis, we tested alternative models including the timing and manipulation treatments (in R, these are lme4::glmer(DNA_conc ~ before_or_after_manipulation*treatment + (1|plot), data = data, family = Gamma(link="log"))). Differences in the ecological community compositions were visualized using t-distributed stochastic neighbor embedding (t-SNE) (Van der Maaten and Hinton, 2008). In the present study, we used the term ‘statistical clarity’ instead of ‘statistical significance’ to avoid misinterpretations, according to the recommendations by Dushoff et al., 2019, when p<0.05.

Quantifications of RNA expressions of rice leaves

Request a detailed protocol

The leaf samples were ground under cryogenic conditions using a Multi-Beads Shocker (Yasui Kikai, Osaka, Japan). Total RNA was extracted using the Maxwell 16 LEV Plant RNA Kit (Promega, Madison, WI, USA). RNA concentration was measured using the broad-range Quant-iT RNA Assay Kit (Thermo Fisher Scientific, Waltham, MA, USA). The RNA concentrations were adjusted and 500 ng of RNA was used as the input of each sample for library preparation. Library preparation for RNA-sequencing was conducted using Lasy-Seq (Kamitani et al., 2019) version 1.1 (https://sites.google.com/view/lasy-seq/). The library was sequenced using HiSeq X (Illumina, San Diego, CA, USA) with paired-end sequencing lengths of 150 bp. On average, 8,105,902 reads per sample (±5,494,529 S.D.) were generated (total reads = 875,534,703 reads; Supplementary file 7).

All obtained reads were trimmed using Trimmomatic version 0.33 (Bolger et al., 2014) using the following parameters: TOPHRED33, ILLUMINACLIP:TruSeq3-SE.fa:2:30:10, LEADING:19, TRAILING:19, SLIDINGWINDOW:30:20, AVGQUAL:20, MINLEN:40, indicating that reads with more than 39 nucleotides and average quality scores over 19 were reported. Then, the trimmed reads were mapped onto the reference sequences of the IRGSP-1.0_transcript (Kawahara et al., 2013) and the virus reference sequences, which were composed of complete genome sequences of 7,457 viruses obtained from NCBI GenBank (Kashima et al., 2021) using RSEM version 1.3.0 (Li and Dewey, 2011) and Bowtie version 1.1.2 (Langmead et al., 2009) with default parameters. The output of the analysis, ‘expected_counts’, was used as inputs to the analysis of differentially expressed genes.

Detection of differentially expressed genes (DEGs)

Request a detailed protocol

The expected counts were imported as a phyloseq object using ‘phyloseq’ package (McMurdie and Holmes, 2013) of R (R Development Core Team, 2022). Then, the object was converted to a DESeq2 object using phyloseq::phyloseq_to_deseq2() function, and differentially expressed genes (DEGs) were detected using ‘DESeq2’” package (Love et al., 2014) of R. DESeq2 provides statistical frameworks for determining differential expression using a model based on the negative binomial distribution. Briefly, size factors and dispersions were estimated, generalized linear models (GLMs) assuming a negative binomial error distribution were fitted, and the differences were tested. These analyses were performed using DESeq2::DESeq() function implemented in the package. In the DESeq2 analysis, false discovery rate was set as 0.05, and genes with an adjusted p<0.05 found by DESeq2 were assigned as differentially expressed. We measured rice gene expression four times and we detected almost no DEGs at sampling events 1, 3, and 4 except for two DEGs at sampling event 4. Therefore, we report only the results of sampling event 2 in the present study (one day after the third field manipulation). In Supplementary files 8 and 9, we report a list of nuclear DEGs.

Data availability

Computer codes and data used in the study are available in Github (https://github.com/ong8181/rice-ecolnet-2019, copy archived at Ushio, 2023) and archived in Zenodo (https://doi.org/10.5281/zenodo.8124471). Sequence data were deposited in DDBJ Sequence Read Archives (DRA). The accession numbers are as follows: DRA009658, DRA009659, DRA009660 and DRA009661 for eDNA data of ecological communities in 2017 (Ushio, 2022), DRA015682, DRA015683, DRA015685, and DRA015686 for eDNA for eDNA data of ecological communities in 2019, and DRA015706 for rice RNA expression data.

The following data sets were generated
    1. Ushio M
    (2023) DDBJ Sequence Read Archive
    ID DRA009658. eDNA data of ecological communities in 2017.
    1. Ushio M
    (2023) DDBJ Sequence Read Archive
    ID DRA009659. eDNA data of ecological communities in 2017.
    1. Ushio M
    (2023) DDBJ Sequence Read Archive
    ID DRA009660. eDNA data of ecological communities in 2017.
    1. Ushio M
    (2023) DDBJ Sequence Read Archive
    ID DRA009661. eDNA data of ecological communities in 2017.
    1. Ushio M
    (2023) DDBJ Sequence Read Archive
    ID DRA015682. eDNA data of ecological communities in 2019.
    1. Ushio M
    (2023) DDBJ Sequence Read Archive
    ID DRA015683. eDNA data of ecological communities in 2019.
    1. Ushio M
    (2023) DDBJ Sequence Read Archive
    ID DRA015685. eDNA data of ecological communities in 2019.
    1. Ushio M
    (2023) DDBJ Sequence Read Archive
    ID DRA015686. eDNA data of ecological communities in 2019.
    1. Ushio M
    (2023) Zenodo
    ong8181/rice-ecolnet-2019: v0.9.0.
    https://doi.org/10.5281/zenodo.8124472

References

    1. Folmer O
    2. Black M
    3. Hoeh W
    4. Lutz R
    5. Vrijenhoek R
    (1994)
    DNA primers for amplification of mitochondrial cytochrome c oxidase subunit I from diverse metazoan invertebrates
    Molecular Marine Biology and Biotechnology 3:294–299.
  1. Software
    1. R Development Core Team
    (2022) R: A language and environment for statistical computing
    R Foundation for Statistical Computing, Vienna, Austria.
  2. Book
    1. Takens F.
    (1981) Detecting strange attractors in turbulence In
    In: Rand D, Young L-S, editors. Dynamical Systems and Turbulence, Lecture Notes in Mathematics. Berlin, Heidelberg: Springer. pp. 366–381.
    https://doi.org/10.1007/BFb0091924
    1. Van der Maaten L
    2. Hinton G
    (2008)
    Visualizing data using t-SNE
    Journal of Machine Learning Research: JMLR 9:2579–2605.

Peer review

Reviewer #1 (Public Review):

This manuscript describes identification of influential organisms on rice growth and an attempt of validation. The analysis of eDNA on rice pot and mimic field provides rice growth promoting organisms. This approach is novel for plant ecology field. However current results did not fully support whether eDNA analysis-based detection of influencing organism.

The strength of this manuscript is to attempt application of eDNA analysis-based plant growth differentiation. The weakness is too preliminary data and experimental set-up to make any conclusion. The trials of authors experiments are ideal. However, the process of data analysis did not meet certain level. For example, eDNA analysis of different time points on rice growth stages resulted two influential organisms for rice growth. Then they cultivate two species and applied rice seedlings. Without understanding of fitness and robustness, how we can know the effect of the two species on rice growth.

The authors did not check the fate of two species after introducing into rice. If this is true, it is difficult to link between the rice gene expression after treatments and the effectiveness of two species. I think the validation experiment in 2019 needs to be re-conducted.

As authorized gave answered, no strong rationale to select the two species was found. However, I insist that the method has enough novelty to present to general audiences.

https://doi.org/10.7554/eLife.87202.3.sa1

Reviewer #2 (Public Review):

Most farming is done by subtracting or adding what people want based in nature. However, in nature, crops interact with various objects, and mostly we are unaware of their effects. In order to increase agricultural productivity, finding useful objects is very important. However, in an uncontrolled environment, it coexists with so many biological objects that it is very inefficient to verify them all experimentally. It is therefore necessary to develop an effective screening method to identify external environmental factors that can increase crop productivity. This study identified factors presumed to be important to crop growth based on metabarcoding analysis, field sampling, and non-linear analysis/information theory, and conducted a mesocosm experiment to verify them experimentally. In conclusion, the object proposed by the author did not increase rice yield, but rather rice growth rate.

The authors responded to my general concerns and all of my specific comments. The manuscript has significantly improved. The flow of aims and approaches is more understandable. Extra supplementary material -especially the visual ones, is useful.

I agree with the other reviewers that the study needs more data and evidence. However, this study aims to introduce ecological concepts and advanced statistical methods to the field. Also, most time series analyses require absolute abundance data, but the manuscript provides solutions for the sequencing data.

https://doi.org/10.7554/eLife.87202.3.sa2

Author response

The following is the authors’ response to the original reviews.

eLife assessment

It is very important to find practical and efficient means in order to increase agricultural productivity. Drawing on data from variable field environments, this study provides a useful theoretical framework to identify new factors that could increase agricultural production. There is solid evidence to support the authors' claims, though following the fate of candidate species after introduction into rice fields would have strengthened the study. Plant biologists and ecologists working in nature and fields will find the work interesting.

Thank you so much for your careful evaluation of our manuscript. We are very pleased to hear that you found our framework useful. We have revised our manuscript according to the "Recommendations for the Authors" to improve our manuscript.

Public Review

Reviewer #1 (Public Review):

This manuscript describes the identification of influential organisms on rice growth and an attempt of validation. The analysis of eDNA on rice pot and mimic field provides rice growth promoting organisms. This approach is novel for plant ecology field. However current results did not fully support whether eDNA analysis-based detection of influencing organism.

Thank you so much for evaluating our manuscript. We have carefully read and responded to your comments. We hope our responses resolve your concerns on our study.

The strength of this manuscript is to attempt application of eDNA analysis-based plant growth differentiation. The weakness is too preliminary data and experimental set-up to make any conclusion. The trials of authors experiments are ideal. However, the process of data analysis did not meet certain levels. For example, eDNA analysis of different time points on rice growth stages resulted in two influential organisms for rice growth. Then they cultivate two species and applied rice seedlings. Without understanding of fitness and robustness, how we can know the effect of the two species on rice growth.

We agree with your comments that we did not have the fitness data of the two species and/or rice seedlings. Thus, it is still difficult to obtain deep understanding of the mechanisms of our findings that the species introduced in the system would influence rice growth. Nonetheless, our study demonstrated the effectiveness of our research framework as we found evidence that the species that were discovered by the eDNA monitoring and time series analysis indeed cause changes in the system. We believe that the first step is to show that the framework is workable and that detailed understanding of the mechanisms or genetic pathway was not a focus of our study. To avoid misunderstanding, we have added several explanations regarding this point in L426–431 and L447. For example, in L426, we have added the following statement: "... the detailed dynamics of the two introduced species was unclear (i.e., the fate of the introduced species). This is particularly important for understanding how the introduced organisms affected rice performance...".

The authors did not check the fate of two species after introducing into rice. If this is true, it is difficult to link between the rice gene expression after treatments and the effectiveness of two species. I think the validation experiment in 2019 needs to be re-conducted.

We did not check the fate of the two species (except measuring the eDNA concentrations of the species), and it is true that we cannot show evidence of "how" these two species influence the rice gene expression. Understanding molecular mechanisms of the phenomenon that we found is important (especially from the viewpoint of molecular biology), but our primary objective was to demonstrate that our "eDNA x time series analysis" framework is feasible for detecting previously overlooked but influential organisms. To this end, we believe that we achieved our objective and repeating the validation experiment should be for a different purpose (i.e., for understanding molecular mechanisms). We have clarified these points in L426–431 and L447 as explained above.

Reviewer #2 (Public Review):

The manuscript "Detecting and validating influential organisms for rice growth: An ecological network approach" explores the influence of biotic and abiotic entities that are often neglected on rice growth. The study has a straightforward experimental design, and well thought hypothesis for explorations. Monitoring data is collected to infer relationships between species and the environment empirically. It is analyzed with an up-to-date statistical method. This allowed the manuscript to hypothesize and test the effects most influential entities in a controlled experiment.

Thank you so much for your careful evaluations. We are pleased to see that you evaluated our manuscript positively. We have further revised our manuscript according to your comments and hope the revision has resolved your concerns.

The manuscript is interesting and sets up a nice framework for future studies. In general, the manuscript can be improved significantly, when this workflow is smoothly connected and communicated how they follow each other more than the sequence and dates provided. It is valuable philosophical thinking, and the research community can benefit from this framework.

Thank you for your suggestions. In order to improve the logic flow and readability of our manuscript, we have revised the descriptions of workflow and clarified how the experimental and statistical steps were connected to each other. To do so, we have added brief explanations about what/how we did at the first sentence of Results subsections (some of these explanations were only in Materials and Methods in the original manuscript). Also, we have moved all of the Supplementary Materials and Methods to the main text. We have thoroughly revised the manuscript, and we hope that all the parts of our manuscript have been connected more smoothly than in the original manuscript.

I understand the length and format of the manuscript make it difficult to add more details, but I am sure it can refer to/clear some concepts/methods that might be new for the audience. How/why variables are selected as important parts of the system, a tiny bit of information about the nonlinear time series analysis in the early manuscript, and the biological reasoning behind these statistically driven decisions are some examples.

We have explained how/why variables are selected (in L125), added more information about the nonlinear time series analysis (in L129 and L175) , and added the biological reasoning behind the statistical decisions (L195).

Reviewer #3 (Public Review):

Most farming is done by subtracting or adding what people want based in nature. However, in nature, crops interact with various objects, and mostly we are unaware of their effects. In order to increase agricultural productivity, finding useful objects is very important. However, in an uncontrolled environment, it coexists with so many biological objects that it is very inefficient to verify them all experimentally. It is therefore necessary to develop an effective screening method to identify external environmental factors that can increase crop productivity. This study identified factors presumed to be important to crop growth based on metabarcoding analysis, field sampling, and non-linear analysis/information theory, and conducted a mesocosm experiment to verify them experimentally. In conclusion, the object proposed by the author did not increase rice yield, but rather rice growth rate.

Thank you so much for your evaluation of our manuscript. We have revised our manuscript based on your comments, and hope it has been improved compared with the original version.

Strength

In actual field data, since many variables are involved in a specific phenomenon, it is necessary to effectively eliminate false positives. Based on the metabarcoding technique, various variables that may affect rice growth were quantitatively measured, although not perfectly, and the causal relationship between these variables and rice growth was analyzed by using information transfer analysis. Using this method, two new players capable of manipulating rice growth were verified, despite their unknown functions until now. I found this process to be very logical, and I think it will be valuable in subsequent ecological studies.

We are very pleased to see that you found our framework is very logical and potentially beneficial for future ecological studies.

Weaknesses

CK treatment's effectiveness remains questionable. Rice's growth was clearly altered by CK treatment. The validation of the CK treatment itself is not clear compared to the GN treatment, and the transcriptome data analysis results do not show that DEG is not present. The possibility of a side effect caused by a variable that the author cannot control remains a possibility in this case. Even though this part is mentioned in Discussion, it is necessary to discuss various possibilities in more detail.

We agree that the effectiveness of the CK treatment was questionable. We have added some more discussion about this point in L376: "The unclear effects of the CK treatment relative to those of the GN treatment could be due to the relatively unstable removal method (i.e., C. kiiensis larvae were manually removed by a hand net) or incomplete removal of the larvae (some larvae might have remained after the removal treatment)."

Reviewer #1 (Recommendations For The Authors):

Comment #1-1 This manuscript describes identification of influential organisms on rice growth and an attempt of validation. The analysis of eDNA on rice pot and mimic field provides rice growth promoting organisms. This approach is novel for plant ecology field. However current results did not fully support whether eDNA analysis-based detection of influencing organism.

Thank you for your careful evaluations of our manuscript. We are pleased to see you found that our approach is novel. We have revised our manuscript in accordance with your comments, and we hope that the revision and responses resolved your concerns.

Comment #1-2 1. Experimental setting: Authors made up small scale pot system in 2017 and then expanded manipulative experiment. I do not understand how two influencing organism sequences were identified from the single treatment depending on different time points. How they can be convince the two organisms affect the rice growth rather than other biological and environmental factors.

In 2017, we performed an intensive monitoring of the experimental rice plots and obtained large time series data (122-day consecutive monitoring x 5 plots = 610 data points). The time series data were analyzed using the information-theoretic causal analysis. The analysis is critically different from correlational analyses and designed to identify causal relationships among variables. Although we understand that field manipulation experiments are a common and straightforward approach to identify causal relationships among organisms, we chose the "fieldmonitoring + time-series-based causal analysis" approach. This is because, as explained in the main text, there are numerous factors that could influence rice performance, and it is practically impossible to perform manipulative experiments for all the potential factors that could influence rice growth. On the other hand, our "field-monitoring + timeseries-based causal analysis" approach has a potential to identify multiple factors under field conditions, even by the single experimental treatment.

Nonetheless, we must admit that our time-series-based approach still has a chance to misidentify causal factors. Our framework relies on statistics, so the chance of false-positive detection of causality cannot be zero. This was exactly the reason why we performed the "validation" experiment in 2019. To complement the statistical results of the 2017 experiments, we performed another experiment in 2019.

Comment #1-3 2. eDNA technology: The eDNA analysis based on four universal primers 16s rRNA, 18s rRNA, ITS, and COI regions must not be enough to identify specific species. The resolution of species classification may not meet to confirm exact species. Thus, the accuracy of two species that they selected for further experiment is difficult to be confirmed. Authors also referred to "putative Globisporangium".

Your point is correct. The DNA barcoding regions we selected are short and it is often difficult to identify species. However, this limitation could not have been overcome even if we had chosen a different genetic marker. The long-read sequencing technology could partially solve the issue, but the number of sequence reads generated by the long-read technique is less than that by the short-read sequencing technology, and comprehensive detection of all species in an ecological community was still challenging. Our approach struck a balance among the identification resolution, comprehensiveness of the analysis, and sequencing costs. In addition, even though we could not identify most ASVs at the species level, some ASVs could be identified at the species level (52 ASVs among the 718 ASVs which had causal influences on rice growth), and we selected the two species (G. nunn and C. kiiensis) from the 52 species.

Further, the taxa assign algorithm we used here (i.e., Claident; Tanabe & Toju 2012 PLoS ONE 10.1371/journal.pone.0076910) adopted conservative criteria for species identification and has a low falsepositive probability.

More importantly, this is also the reason why we performed the "validation" experiment in 2019. The species identified in the 2017 experiment are still "potential" organisms that influence rice growth (i.e., the hypothesis-generating phase), and we tested the hypothesis in 2019.

Nonetheless, we must admit that clear description of potential limitations is important. Thus, we have discussed this in L418: "As for the second issue, short-read sequencing has dominated current eDNA studies, but it is often not sufficient for lower-level taxonomic identification. Using long-read sequencing techniques (e.g., Oxford Nanopore MinION) for eDNA studies is a promising approach to overcome the second issue".

Comment #1-4 3. Biological relevance 1: Authors identify two organisms as influencing organism for rice growth. As conducting the first experiment in 2017, the 2019 experiment was different from natural condition. The two experiments in 2017 and 2019 were conducted under different conditions. How do they compare the experiments? At least, the eDNA analyses in 2017 and 2019 should be very similar. I cannot find such data.

The experimental conditions were different between 2017 and 2019 because they were conducted in different years. Theoretically, it is ideal if the experimental conditions in 2019 are covered by the range of experimental conditions in 2017 (e.g,. rice variety, air temperature, rainfall, and solar radiation). If this condition were satisfied, the attractor (i.e., rice growth trajectory delineated in the state space) in 2019 would be within that in 2017, and our model prediction in 2017 would be used to predict dynamics in 2019 accurately. To fulfill the conditions, we made as much effort as possible: we used the same rice variety and soils in 2019 as those used in 2017, and started our experiment at the same timing in 2019 as that in 2017.

Although natural ecological dynamics cannot be precisely controlled, our monitoring revealed that the ecological dynamics in 2019 was qualitatively similar to that in 2017. To demonstrate that the experimental conditions and eDNA community data were similar between the two experiments, we have presented the climate and eDNA data in an inset figure in Figure 3a, Figure 1–figure supplement 2, Figure 3–figure supplement 2. We must admit that these dynamics are not identical, but we hope that this resolves your concern.

Comment #1-5 4. Lack of detail description: In the Materials and Methods, there are many parts which lack on detail description. For instance, authors must described the two species cultivation, application concentrations, and application methods.

We have moved Supplementary Materials and Methods to the main text and added more detailed descriptions in Materials and Methods. Also, to improve the logical flow and readability of our manuscript, we have added brief explanations about what/how we did at the first sentence of Results subsections (some of these explanations were only in Materials and Methods in the original manuscript). We have added the reference for how to cultivate G. nunn in L608 (Kobayashi et al., 2010; Tojo et al., 1993) (C. kiiensis was not cultivated but removed from the system as in Materials and Methods), and application concentrations. Application methods were described in Materials and Methods, the section Field manipulation experiments in 2019 in L596.

Comment #1-6 5. Validation: Application of one species clearly resulted to promote rice growth. They must include appropriate control treatment. If they pick same genus but different species that identified no specific effect on rice growth through eDNA analysis, no effect on growth can be provided. Generally application of large population of certain non-harmful organism confer plant growth promotion. It is not surprising result. Authors need to prove effectiveness of eDNA analysis. In addition, the field experiments required at least two years of consistent data for publication because environmental factors are so dynamic.

Thank you for pointing this out. We agree with your comment that species that were predicted to have no effect should not promote rice growth in a validation experiment. It was also one of our inititial experimental plans to include such species in our manipulation experiment in 2019, but we could not include them because of the limitation of time, labor, and money. More extensive validation of the statistical results of the 2017 data, including multi-year experiments, would further validate the effectiveness of our approach, which should be done as future studies. To clarify this point, we have added statements in the paragraph starting at L396.

Comment #1-7 In conclusion, I suggest that authors need more large data analysis and validate with more accurate and meaningful protocol.

As we explained in the revised manuscript and the Response to Comments #1-2 to #1-7, our study demonstrated a novel research framework to detect previously overlooked influential organisms under field conditions. We agree that larger data analysis would be ideal to further validate our approach, but whether and how to collect larger data is constrained by time, money, and labor. We believe that our study was designed carefully and could provide meaningful avenues for developing an ecological-network based, novel, and environment-friendly agriculture solutions.

Reviewer #2 (Recommendations For The Authors):

Comment #2-1 Lines 97-110: This is so cool. Modeling with empirical data is very powerful. But a rice field is an open system consisting of metacommunity dynamics. Maybe a tiny bit of biological and biogeochemical background here would be good.

Thank you for your comments. We have added a few examples of how and in which systems these methods were used to evaluate community dynamics and detect biological interactions in L109-L118.

Comment #2-2 Lines 111-126: I like the summary of the study here. I think the influential species concept can be a little more elevated. Paine's famous keystone species work has been cited but a couple more pieces of literature can help to enhance the ecological importance of this work.

We have explained the work by Paine (1966) a bit more and added one more paper that showed the effect of multiple predator species on the system dynamics at L88. We have also added a relevant sentence at L137 to emphasize the ecological/agricultural significance of our work.

Comment #2-3 Experimental design/Figure 1:

Is there any rationale behind choosing red individuals to measure the growth?

Is there any competition between the individuals in the pots?

Figure 1e: It is nice to show the ASVs in time. I wonder how the plot would look like when normalized by biomass/DNA content/coverage/rarefaction because of the seasonality.

As for the first question, we chose the four individuals to minimize the edge effects (i.e., effects of microclimates and neighboring rice would be different between the four rice individuals and those planted in the edge regions). We have mentioned this in the legend of Figure 1.

As for the second question, there might be competition among the individuals in the pot. However, we did not measure the effect of competition (e.g., by comparing the growth with/without other rice individuals).

As for the third question, we published detailed dynamics of ecological community in the Supplementary Figures in Ushio (2022) Proceedings Bhttps://doi.org/10.6084/m9.figshare.c.5842766.v1. In addition, we haveuploaded a video showing the temporal dynamics of some top ( = most abundant) ASVs in https://doi.org/10.6084/m9.figshare.23514150.v2.

We have mentioned the supporting information in L153.

Comment #2-4 Line 146-147: Is this damage influence the inferences? Maybe it is better to justify.

While we occasionally observed physical damages, it is unlikely that they affected our causal inference because the changes in the rice heights due to the damages were smaller and less frequent than those due to growth. We have noted this at L151.

Comment #2-5 Line 161-162: Maybe refer readers to the methods section where you explain UIC analysis. It'd be easier to interpret the figures.

Mentioned.

Comment #2-6 Line 175-176: I believe very brief information in the intro about the organisms might help explain the hypothesis and interpret the results better.

We have included brief information of the two species at L197.

Comment #2-7 Figure 2: Species interaction strength: Are these proxies to the Jacobians? Is there a threshold for the influence we can consider strong/weak? For example, influential species compared to diagonal elements of the Jacobians (intraspecies interactions) could be shown as a mean vertical line in Figure 2b.

"Influences to rice growth" in Figure 2b is transfer entropy (TE) from a target ASV to rice growth. They are not proxies of the Jacobians, but they might positively correlate with the absolute value of the Jacobians. We have clarified this point in the legend (L953). More direct estimations of the Jacobian can be done using the MDR S-map method (Chang et al. 2021 DOI:10.1111/ele.13897), but we did not perform the MDR S-map in the present manuscript (see Ushio et al. 2023https://doi.org/10.7554/eLife.85795 for the application of the MDR S-map). As for TE, there is no clear threshold to distinguish strong/weak interactions.

Comment #2-8 Figure 2: Looking at panels c and d, it looks like there is a negative frequency selection between two influential species. Is it a reasonable observation?

This is an interesting point. In this manuscript, we have not carefully examined the interspecific relationship between these two particular species. However, the interspecific interactions were examined in detail and reported in Ushio (2022) Proceedings of the Royal Society B (DOI:10.1098/rspb.2021.2690). We re-checked the result in Ushio (2022); although there is a negative correlation between them, we did not find any (statistical) causal relationship between them.

Comment #2-9 Line 209: What is t-SNE analysis? Because of the manuscript's format, maybe methods should be shortly referred to in the relevant section or explained in brackets.

We have spelled out t-SNE.

Comment #2-10 Line 212-214: Maybe briefly explain what the hypotheses are for the alternative analysis, and what is the contribution of the results to the study.

We have added a brief explanation at L241: "Alternative statistical modeling that included the treatments (the control versus GN or CK treatments) and manipulation timing (i.e., before or after the manipulation), which simultaneously took the temporal changes of all the treatments into account, also showed qualitatively similar results (Supplementary file 4), further supporting the results."

Comment #2-11 Figure 3b/c: Maybe species names as panel titles could be helpful. d: Treatment names with initials in the legend could be also helpful to read the plots.

We have added species name as panel titles of Figure 3b,c. Treatment names were included in the legend of Figure 3.

Comment #2-12 Line 233: Maybe mention why the manuscript uses the word "clear".

We have mentioned this in L185.

Comment #2-13 Line 234-236: I think that these alternative tests should be explained somewhere.

We have revised the sentence so that it includes some explanations (L241). Also, we have referred to Materials and Methods.

Comment #2-14 Figure 4: The title says ecological community compositions, and panels show the growth rates and cumulative growth.

Thank you for pointing this out. This was a typo and we have corrected it.

Comment #2-15 Lines 246-269: Can these expression patterns be transient and relevant to the time point that the sample is taken?

Yes, these expression patterns were transient. We collected rice leaf samples for RNA-seq 1 day before the first manipulation and 1, 14, and 38 days after the third manipulation (see Supplementary file 3 for the sampling design). When we merged the pot locations, we observed no difference in the gene expression for samples 1 day before the first manipulation and 14 and 38 days after the third manipulation (except for two genes in samples 38 days after the manipulation), and thus, we consider the DEGs that appeared only in the short period after the manipulation. We have mentioned this in L278 and L383: "We found almost no DEGs for leaf samples taken one day before and 14 and 38 days after the third manipulation (the leaf sampling event 1, 3, and 4), suggesting that the influences of the treatments on the gene expression patterns were transient." (L278) and "These changes were observed relatively quickly and transient." (L383)

Comment #2-16 I wonder if a conceptual framework figure would help to generalize the workflow that can be used for other studies.

Thank you for your suggestion. Although we agree with your comment that such a figure would be helpful to generalize the workflow, we believe that our framework is clear and decided not to include it in the present manuscript. We might consider including such a figure (like Figure 1a in Ushio 2022) if we have an opportunity to write a review paper regarding this topic.

Comment #2-17 Lines 329-335: I feel this information is unclear in the early manuscript.Maybe it's necessary to clearly communicate in the beginning.

We have explained that we could not find any relevant information at least at the time we detected the ASVs in L189.

Comment #2-18 Lines 336-337: Can these species be identified in the previous data set from the ASV sequences?

Yes, these species were identified in the DNA data set obtained in 2017.

Comment #2-19 Lines 387-397: Are there any measurements such as total biomass, and statistical methods to help with the eDNA bias and data compositionality?

We have confirmed that our quantitative eDNA metabarcoding generates comparable results with the fluorescence-based method and quantitative PCR (e.g., see Supplementary Figures in Ushio 2022) (mentioned in L310 in the revised manuscript). However, at least in this study, we could not perform a direct comparison of the eDNA data with species abundance and/or biomass. This is partly because the number of our target species was too large (> 1,000 species). The accurate estimation of species abundance and/or biomass is one of our next goals.

Comment #2-20 Line 472: Maybe mention transfer entropy somewhere in the early manuscript.

We have mentioned this in L175.

Comment #2-21 Lines 494-503: Maybe a summary of this reasoning should be mentioned somewhere in the early manuscript too.

We have described a brief summary of the reasoning in L195.

Comment #2-22 Lines 29-33 If this sentence is simplified it might be easier to follow.

The sentence has been divided into two sentences in L28. Also, each sentence has been simplified.

Comment #2-23 Line 38 Maybe "macrobes" can be explicitly mentioned. Fungi, protozoa, etc.

Mentioned.

Comment #2-24 Line 139: I am not sure if the date should be in the title.

Similar monitoring was done in 2017 and 2019. Thus, we think the date is necessary in the section title.

Comment #2-25 Figure 1: There are 4 red individuals in the design but 5 measurements in the plots.

Heights and SPAD of the four individuals were measured for each plot and the averaged values were used as representative values for each plot. Therefore, 20 measurements ( = 4 rice individuals 5 plots) were done every day, but each plot has one rice height for each day. We have clarified this in the legend of Figure 1: "the average values of the four individuals were regarded as representative values for each plot."

Comment #2-26 Figure 1b: Maybe use the same axis length for the temperature as the other plots?

Corrected.

Comment #2-27 Lines 259-261: Are there the names of the genes in databases?

Yes, these are gene names used in the rice databases (e.g., The Rice Annotation Project Database; https://rapdb.dna.affrc.go.jp/inde x.html).

Reviewer #3 (Recommendations For The Authors):

Comment #3-1 Additionally, RGR is not statistically significant, but statistical significance is observed only in cumulative growth because data presentation does not reflect plant characteristics. RGR changes according to the developmental stage of the plant. Therefore, if RGR data are shown separately according to the rice growing season, the cumulative growth pattern and the pattern will appear similar.

RGRs were calculated daily (i.e., cm/day) and they changed depending on the developmental stage of the rice (Figure 1 and Figure 4–figure supplement 1). Therefore, we might find similar RGR patterns if we focus on a specific period of the growing season. However, unfortunately, we performed the intensive (i.e., daily) monitoring in 2019 only during the field manipulation period (middle June to middle July 2019), and we cannot investigate the changes in cumulative growth throughout the growing season (this depends on how many days we add up RGR to calculate the cumulative growth, though). We agree that, if we had investigated the detailed pattern of RGR throughout the growing season in 2019, we could have found similar pattens between RGR and cumulative growth rate at a certain period in the growing season. In Figure 4, the cumulative growths were calculated based on the RGRs before the third manipulation or during 10 days after the third manipulation. We clarified this in the legend of Figure 4.

https://doi.org/10.7554/eLife.87202.3.sa3

Article and author information

Author details

  1. Masayuki Ushio

    1. Hakubi Center, Kyoto University, Kyoto, Japan
    2. Center for Ecological Research, Kyoto University, Otsu, Japan
    3. Department of Ocean Science, The Hong Kong University of Science and Technology, Clear Water Bay, Kowloon, Hong Kong SAR, China
    Contribution
    Conceptualization, Data curation, Formal analysis, Supervision, Funding acquisition, Validation, Investigation, Visualization, Methodology, Writing - original draft, Project administration, Writing – review and editing
    For correspondence
    ong8181@gmail.com
    Competing interests
    The authors declare that there is no conflict of interest
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-4831-7181
  2. Hiroki Saito

    Tropical Agriculture Research Front, Japan International Research Center for Agricultural Sciences, Okinawa, Japan
    Contribution
    Investigation, Methodology, Writing – review and editing
    Competing interests
    No competing interests declared
  3. Motoaki Tojo

    Graduate School of Agriculture, Osaka Metropolitan University, Osaka, Japan
    Contribution
    Resources, Investigation, Methodology, Writing – review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-9660-4359
  4. Atsushi J Nagano

    1. Faculty of Agriculture, Ryukoku University, Otsu, Japan
    2. Institute for Advanced Biosciences, Keio University, Tsuruoka, Japan
    Contribution
    Software, Investigation, Methodology, Writing – review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-7891-5049

Funding

Japan Science and Technology Agency (JPMJPR16O2)

  • Masayuki Ushio

Japan Society for the Promotion of Science (KAKENHI (B) 20H03323)

  • Masayuki Ushio

The Hong Kong University of Science and Technology (Startup Fund)

  • Masayuki Ushio

Kyoto University (Hakubi Project)

  • Masayuki Ushio

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We thank Asako Kawai and Mutsumi Kato for assistance in field monitoring and DNA library preparations, Yutaka Osada for advice on the data analysis, Akira Matsumoto and Satomi Yoshinami for maintenance of the experimental field and assistance in the field monitoring, Eiso Inoue for comments on the midge species, and Natsumi Fujii for assistance in the field experiment. This research was supported by PRESTO (JPMJPR16O2) from the Japan Science and Technology Agency (JST), KAKENHI (B) 20H03323, the Hakubi Project in Kyoto University, and The Hong Kong University of Science and Technology Startup Fund to MU.

Senior Editor

  1. Detlef Weigel, Max Planck Institute for Biology Tübingen, Germany

Reviewing Editor

  1. Youngsung Joo, Seoul National University, Republic of Korea

Version history

  1. Preprint posted: February 21, 2023 (view preprint)
  2. Sent for peer review: March 3, 2023
  3. Preprint posted: May 5, 2023 (view preprint)
  4. Preprint posted: August 25, 2023 (view preprint)
  5. Version of Record published: September 13, 2023 (version 1)

Cite all versions

You can cite all versions using the DOI https://doi.org/10.7554/eLife.87202. This DOI represents all versions, and will always resolve to the latest one.

Copyright

© 2023, Ushio et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 533
    Page views
  • 80
    Downloads
  • 0
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Masayuki Ushio
  2. Hiroki Saito
  3. Motoaki Tojo
  4. Atsushi J Nagano
(2023)
An ecological network approach for detecting and validating influential organisms for rice growth
eLife 12:RP87202.
https://doi.org/10.7554/eLife.87202.3

Share this article

https://doi.org/10.7554/eLife.87202

Further reading

    1. Ecology
    2. Plant Biology
    Daniel Fuks, Yoel Melamed ... Ehud Weiss
    Research Article Updated

    Global agro-biodiversity has resulted from processes of plant migration and agricultural adoption. Although critically affecting current diversity, crop diffusion from Classical antiquity to the Middle Ages is poorly researched, overshadowed by studies on that of prehistoric periods. A new archaeobotanical dataset from three Negev Highland desert sites demonstrates the first millennium CE’s significance for long-term agricultural change in Southwest Asia. This enables evaluation of the ‘Islamic Green Revolution (IGR)’ thesis compared to ‘Roman Agricultural Diffusion (RAD)’, and both versus crop diffusion during and since the Neolithic. Among the findings, some of the earliest aubergine (Solanum melongena) seeds in the Levant represent the proposed IGR. Several other identified economic plants, including two unprecedented in Levantine archaeobotany—jujube (Ziziphus jujuba/mauritiana) and white lupine (Lupinus albus)—implicate RAD as the greater force for crop migrations. Altogether the evidence supports a gradualist model for Holocene-wide crop diffusion, within which the first millennium CE contributed more to global agricultural diversity than any earlier period.

    1. Ecology
    2. Physics of Living Systems
    William R Shoemaker, Jacopo Grilli
    Research Advance

    The structure and diversity of microbial communities are intrinsically hierarchical due to the shared evolutionary history of their constituents. This history is typically captured through taxonomic assignment and phylogenetic reconstruction, sources of information that are frequently used to group microbes into higher levels of organization in experimental and natural communities. Connecting community diversity to the joint ecological dynamics of the abundances of these groups is a central problem of community ecology. However, how microbial diversity depends on the scale of observation at which groups are defined has never been systematically examined. Here, we used a macroecological approach to quantitatively characterize the structure and diversity of microbial communities among disparate environments across taxonomic and phylogenetic scales. We found that measures of biodiversity at a given scale can be consistently predicted using a minimal model of ecology, the Stochastic Logistic Model of growth (SLM). This result suggests that the SLM is a more appropriate null-model for microbial biodiversity than alternatives such as the Unified Neutral Theory of Biodiversity. Extending these within-scale results, we examined the relationship between measures of biodiversity calculated at different scales (e.g. genus vs. family), an empirical pattern previously evaluated in the context of the Diversity Begets Diversity (DBD) hypothesis (Madi et al., 2020). We found that the relationship between richness estimates at different scales can be quantitatively predicted assuming independence among community members, demonstrating that the DBD can be sufficiently explained using the SLM as a null model of ecology. Contrastingly, only by including correlations between the abundances of community members (e.g. as the consequence of interactions) can we predict the relationship between estimates of diversity at different scales. The results of this study characterize novel microbial patterns across scales of organization and establish a sharp demarcation between recently proposed macroecological patterns that are not and are affected by ecological interactions.