Ecology

Collection of marine organismsMarine invertebrates such as Corallium rubrum are ideal organisms to perform controlled experiments and to gather useful information on a variety of environmental conditions74. This species, whose diet is based on small zooplankton captured with the polyp tentacles, has been already used in long-term experiments74,75,76. Coral specimens were collected in March 2017 at ca. 35-m depth in the Marine Protected Area of Portofino (Punta del Faro, 44°17′41.02 N; 9°13′31.30 E) in the Ligurian Sea (North-Western Mediterranean Sea) by scuba divers (using TRIMIX blending).Experimental designAfter recovery, the coral specimens were brought to the laboratory and maintained in a tank (30 L) at in situ temperature (13 ± 0.8 °C) and subjected to the continuous flux of natural seawater filtered onto 0.7-µm pore-size membranes (micro-glass fibre paper, Munktell) by using two submersible pumps (Euronatale, 203 V, 50 Hz, 4 Watt).Sixty coral branches obtained from different colonies, with similar morphology, and a surface of ~2 cm2 each, were distributed among 12 experimental tanks, in order to have 5 coral branches per tank (12 L glass tanks, containing on average, 274 ± 26.4 coral polyps each). The corals were acclimatised for 20 days in a temperature-controlled room, and dim light conditions, before starting experiments. Each tank, filled with natural seawater, was equipped with a prefiltered (0.2 µm) channelled aeration system combined with motor-driven paddles in order to create convective currents, which allowed the resuspension of the microplastic mixture, thus ensuring as much as possible a homogeneous distribution of the polymers. This experimental system was designed and set up according to Sutherland et al.77. To assess the potential effects of increasing microplastic  microparticles L−1 (here defined as low, medium and high concentrations of microplastic particles). We also quantified the exact amount of particles actually interacting with the corals, by discounting the fractions loss due to experimental manipulations (see details in Supplementary Methods). According to the results reported in the Supplementary Results, the systems were responsible for the loss of ca. 40% of the microplastic particles, thus the corals in experimental systems were actually exposed to 60, 300 and 600 microplastic particles per litre (to which we referred the low, medium and high concentrations).The highest concentrations of microplastics (up to 600 microplastic particles per litre) can reflect future contamination on the basis of estimates obtained by numerical models9, whereas the low and medium concentrations have been selected to represent highly-contaminated marine habitats, including the areas where the corals were collected (Ligurian Sea)78,79. In particular, for the Ligurian Sea, we estimated an average value of 94 microplastic particles L−1, based on the concentrations of microplastic particles ( >200 µm) determined by Fossi et al.78,79, and the most cautionary correction factor (105) calculated by Brandon et al.10 for the unaccounted smaller fraction of microplastics (25–75% of the fragments falling approximately in the 20–100 µm dimensional class with median range: 59–116).Microplastic mixtures were also prepared considering the concentration and composition of dominant polymers in different coastal marine environments, especially in hot spots of microplastic contamination5,7,8.The microplastic mixture added to the tanks was composed of 76.6% polyethylene, 10.9% polypropylene, 7.3% polystyrene, 3.3% polyvinylchloride and 1.8% polyethylene terephthalate particles. These particles were obtained by milling plastic objects from everyday life (i.e., containers, bottles, cups, pipes) according to Paul-Pont et al.8 (Supplementary Table 5). Plastic milling was carried out under a laminar flow hood in chilled sterilized and 0.02 µm prefiltered milliQ water. All the tools used for handling plastics were pre-treated with 1% sodium hypochlorite in water and rinsed 10 times with sterilized and 0.02 µm prefiltered milliQ water, and then dried under laminar flow hood. Details on the preparation of microplastic mixtures are reported in the Supplementary Methods.The low, medium and high concentrations of microplastics were added in triplicate tanks (n = 3 for each concentration). Additional systems containing seawater and coral branches without microplastics (n = 3, here defined Controls), and seawater added with microplastics (at the highest concentration) without red corals (n = 3, here define CTRL MPs) were used as controls. Overall, the experimental setup comprised 15 tanks.The experiments for assessing the impact of microplastics on red corals started immediately after the microplastic mixture addition (time 0). During the experiment, seawater temperature (range: 13.10 ± 0.01–13.13 ± 0.05 °C), salinity (range: 38.35 ± 0.18–38.65 ± 0.18) and oxygen levels (7.10 ± 0.08–7.36 ± 0.2 mg L−1) were monitored daily in all tanks using a probe (YSI Professional Plus, USA) and corals were fed three times a week with 103 Artemia salina nauplii L−1.After ten days the condition of corals that were exposed to microplastics was deteriorating, so we collected one coral branch from each tank for molecular analyses (i.e., associated microbiome, gene expression and DNA damage). After 14 days, the experiment was stopped because the coral branches that were exposed to the medium and high microplastic particle concentrations were completely wrapped in mucus, with a large portion of damaged tissue and without polyp activity, therefore corals were defined dead (overall 12 branches, see ‘Results’ section for details). Coral branches in the controls showed no visible signs of necrosis or other macroscopic stress. The tissue remained intact, and the colour unchanged until the end of the experiment.Effects of microplastic ingestionCoral feeding activityTo assess the impact of microplastics on feeding activity of C. rubrum, analyses based on the use of Artemia salina were performed after 2 and 10 days from the start of the experiment (t0) in replicate systems (n = 3 for each treatment, n = 3 for the controls) according to standard international protocols80. The nauplii of Artemia salina were reared in the laboratory, incubating 0.5 g of cysts (Ocean Nutrition) in 1 L of seawater filtered onto 0.2-µm filter in a separatory funnel, 2 days before the analysis of feeding rate. At hatching, nauplii were counted and maintained in vials to obtain the concentration of 1000 nauplii L−1. To avoid stress, corals (one branch for each tank) were transferred underwater to beakers along with 1 L seawater of each tank. After addition of live A. salina nauplii (1000 nauplii L−1) to the 1 L beakers containing the coral branches and to the controls, three aliquots of 10 ml seawater were collected after ~30 s from the start of the experiment (t0) and after 2 and 4 h. The remaining nauplii present in each seawater aliquot were counted under a stereomicroscope at ×3.2 magnification (Zeiss Stami 2000). Mean ingestion rates (nauplii removed h−1) were determined by linear regression analysis.Accumulation of microplastics by C. rubrum
To investigate the accumulation of plastic polymers by C. rubrum polyps, the number of microplastic particles ingested by coral polyps was evaluated after 14 days of exposure to microplastic mixture, by dissolving polyps and skeleton of the corals (one for each tank at the concentration of 1000 microplastic particles L−1) using an acid/base digestion protocol36 with some modifications.To exclude biases on the estimate of the number of microplastic particles actually accumulated within the polyps, coral branches were accurately rinsed with milliQ water and checked under stereomicroscope (at ×50 magnification) for the potential presence of microplastic particles adherent to the coral tissue. Coral branches were then soaked in 5 ml of 4.5% sodium hypochlorite (NaClO) for 24 h and dissolved in 5 ml of 37% HCl for 30 min. Particulate material was retained on a 0.2-μm filter in a vacuum filtration system, and microplastic particles were counted under a stereomicroscope at ×50 magnification. The chemical composition of the polymers ingested by corals was confirmed by FT-IR analyses (Perkin Elmer, software Packages Spectrum 5.3.1). To evaluate possible damage to plastic polymers due to the use of acid/base solutions, we exposed polypropylene, polyethylene, polystyrene, polyvinylchloride and polyethylene terephthalate at the same volume and concentration of NaClO and HCl during digestion of the coral.Potential transfer of microplastics by zooplanktonWhile testing the exposure of the red corals to microplastics, we also determined the rate of microplastic ingestion by A. salina nauplii used to feed the red corals, in order to assess their role as potential vectors of microplastics. To do this, additional tanks (n = 3) were added with 0.2 µm prefiltered 12 L natural seawater, 1000 nauplii L−1 of A. salina and the same microplastic mixture used for the experiment on the red corals (at the highest concentration). Three other tanks were used as controls containing 0.2 µm prefiltered 12 L natural seawater and 1000 nauplii L−1 of A. salina.Microplastic ingestion by A. salina was determined after 2 and 10 days of experiment following the enzymatic digestion protocol previously developed81 with some modifications. Such a procedure degrades biological tissues without affecting shape, colour and composition of plastic fragments. Gut contents of 100 individuals of A. salina (n = 5) were assessed under a stereomicroscope (Leica MZ125) and light microscope (Zeiss Axiovert 200) and photographed with a Zeiss Axiocam digital camera. Afterwards, nauplii were processed immediately according to the modified enzymatic digestion protocol. Nauplii were dried in an oven for 3 h at 60 °C, transferred to glass jars containing a buffer homogenizing solution (400 mM Tris-HCl pH 8, 60 mM EDTA pH 8, 5 M NaCl, SDS 1%) incubated at 50 °C for 15 min and exposed to Proteinase K (1 mg ml−1). Then, samples were dried for 2 h at 50 °C, homogenized and re-incubated at 60 °C for 20 min and sonicated on ice (1–2 min) three times. After digestion, the microplastic-containing suspensions were placed in Utermöhl chambers and the microplastics were examined at the inverted light microscope (Leica DMI3000-Bat ×200 magnification) and counted. Microplastics obtained from nauplii digested after 10 days of incubation were also measured and categorized by colours and shape to evaluate the numbers and the size spectra of microplastics ingested by A. salina during the experiment.Physical impact on coral coenenchymaScanning electron microscopy (SEM) analysesTo investigate the physical damage of the microplastic mixture on the coral tissues, samples (one branch from each tank including the control) were collected before the start of the experiment (t0), after 7 days and at the end of the experiment and prepared for SEM analyses according to standard protocols82 with some modifications. Coral branches were stored in 0.7 µm prefiltered seawater with 4% buffered formalin. After 24 h, samples were washed with 0.7 µm prefiltered seawater and dehydrated for 3 h in 20% ethanol. After 3 h they were washed in the same way and dehydrated in ethanol 50%. After 3 h, samples were stored in 70% ethanol. Samples were stored at +4 °C. We dehydrated samples using different gradients of ethanol solutions (70–80%, 80–90%, 90–95%, 95–99% in 2 days)82. Then, samples were dried using HMDS (Hexamethyldisilazane, Aldrich 440191)83. Dried samples were mounted on aluminium stubs using Leit-C glue (conductive carbon cement, Neubauer Chemikalien) and sputter-coated with gold. Samples were examined with a Scanning Electron Microscope (Zeiss SUPRA 40). In addition, the tissue damage percentage was assessed on SEM micrographs at ×200 of magnification by using PhotoQuad v1.4 software84. Such a software for advanced image processing of 2D photographic quadrat samples, dedicated to ecological applications, was used for the analysis of three randomly selected areas from the apex to the base of each coral rotating it on three sides (n = 9). Additional analyses through random SEM observations (n = 20) at 3.00KX to 17.00KX of magnification were carried out to determine prokaryotic cell abundances around lesions of corals (n = 3) exposed to high concentrations of microplastic particles. Data were standardised to the coral surface analysed.Mucus release and trapped microplastics and prokaryotic cellsTo evaluate the first symptoms of coral stress, a photographic report was conducted daily. The abundance of microplastic particles trapped in coral mucus was estimated using an enzymatic digestion protocol81 with some modifications. Mucus produced by corals exposed to higher microplastics concentrations was dried in oven at 60 °C for 12 h. After 12 h, five ml of homogenizing solution was added to the samples and incubated at 50 °C for 15 min. Proteinase K (1 mg mL−1) was added to the samples, which subsequently were incubated at 50 °C for 2 h. Then, samples were homogenized and incubated again at 60 °C for 20 min, after that samples were sonicated three times (three 1-min treatments using a Branson Sonifier 2200; 60 W). After digestion, microplastics-containing suspension was filtered on 0.2-μm filters in a vacuum filtration system (Whatman, Nuclepore). Filters were analysed at stereomicroscope at ×50 magnification (Zeiss Stemi 2000).Stress signals at the molecular levelRNA extraction, cDNA synthesis and gene expression level by qPCRTo assess potential changes in the gene expression pattern of C. rubrum due to microplastics, total RNA was extracted from ca. 20 mg of tissue (wet weight) from one coral branch randomly collected from each treatment (n = 3) and control (n = 3) after 10 days of experiment by using Quick-RNA™ MiniPrep (Zymo Research, Freiburg, Germany) according to the manufacturer’s instructions. Total RNA was also extracted from additional samples of coral branches collected randomly at the beginning of the experiment. Once scraped by surgical disposable scalpels (Braun), coral tissues were placed in new 2 ml sterile tubes and washed three times with phosphate-buffered saline (PBS 1×). Samples were centrifuged at 1800 rpm for 10 min in an Eppendorf® 5810r refrigerated centrifuge using a swing-out rotor at 4 °C and, after removing the supernatant, were homogenized for 5 min with a RNase-free sterile glass stick in RNA lysis buffer. Contaminating DNA was degraded by treating each sample with DNase dissolved in RNase-free water included in the kit. For each sample, 250 ng of total RNA extracted was retrotranscribed with an iScript™ cDNA Synthesis kit (Bio-Rad, Milan, Italy), following the manufacturer’s instructions. The reaction was performed on the Veriti™ 96-Well Thermal Cycler (Applied Biosystem, Monza, Italy). To evaluate the efficiency of cDNA synthesis, a PCR was performed with primers of the reference gene, cytochrome oxidase I (COI, Supplementary Table 6). The reaction was carried out using MyTaq™ HS DNA Polymerase (Bioline, Luckenwalde, Germany) on the Veriti™ 96-Well Thermal Cycler (Applied Biosystem, Monza, Italy). The PCR programme consisted of a denaturation step at 95 °C for 1 min, 35 cycles at 95 °C for 45 s, 60 °C for 45 s, and 72 °C for 45 s and a final extension step at 72 °C for 10 min.The expression levels of the six genes of hsp70, hsp60, MnSOD, mtMutS, EF1 and cytb, involved in a broad range of functional responses, such as stress, detoxification processes, and DNA repair, were followed by real-time qPCR to identify potential stress of corals exposed to microplastics61. For the cytb, target-specific primer pairs were designed with the Primer 3 software (http://primer3.ut.ee85) using nucleotide sequences retrieved from the GenBank database for C. rubrum as template (https://www.ncbi.nlm.nih.gov/genbank/; Supplementary Table 6). SensiFAST™ SYBR® & Fluorescein mix (Bioline, Luckenwalde, Germany) were used for measuring the levels of mRNAs on CFX Connect™ Real-Time PCR detection system (Biorad, Milan, Italy). Fluorescence was measured using CFX Manager™ software (Biorad, Milan, Italy). All genes tested by qPCR in this study were amplified with primers purchased from Life Technologies/Thermo Fisher Scientific (Milan, Italy). The fold change in target gene mRNA expression of corals exposed to microplastics compared with the control was calculated using the comparative CT method using the 2−ΔΔCt equation86. COI was used as reference gene for normalising the gene expression analyses.DNA oxidative damageFor evaluating oxidative DNA damage potentially due to microplastic exposure on C. rubrum, the content of 8-hydroxydeoxyguanosine (8-OHdG) was analysed. DNA was extracted from 20 mg (wet weight) of tissue randomly collected from one coral branch for each treatment (n = 3) and control (n = 3) after 10 days of experiment using DNeasy Blood & Tissue Kits (Qiagen, Valencia, CA) and following the manufacturer’s protocol. Finally, samples were kept at −20 °C before subsequent analyses. Nucleic acids extracted (2 μg) were transferred into new 2-ml tubes and incubated for 5 min at 95 °C, then rapidly chilled on ice. Samples were digested to nucleosides by incubating the denatured DNA in sodium acetate 20 mM, pH 5.2 with 2 μl of nuclease P1 (6 U/μl; Merck KGaA, Darmstadt, Germany) for 2 h at 37 °C. Each sample was then incubated with 5 μl alkaline phosphatase (1 U/μl; Roche, Mannheim, Germany) in Tris-HCl 100 mM, pH 7.5 for 1 h at 37 °C. The reaction mixtures were then centrifuged for 5 min at 6000 × g and the supernatants tested for DNA oxidation with an OxiSelect™ Oxidative DNA Damage ELISA Kit (8-OHdG Quantitation; Cell Biolabs, CA, USA). As positive control, Escherichia coli genomic DNA (2 μg) was incubated in a final concentration of 50 and 100 mM H2O2 overnight at 37 °C, and subsequently tested.Prokaryotic abundance in coral mucus and surrounding seawaterTo highlight possible effects in terms of prokaryotic contamination associated with the exposure of the corals to microplastics, we determined prokaryotic abundances in the mucus released by C. rubrum and the surrounding seawater.Prokaryotic abundances in the coral mucus collected from each tank (except for the control where coral mucus was not released) after 14 days of the experiment, were analysed by epifluorescence microscopy. The extraction of prokaryotic cells from the mucus (ca. 1 mL for each tank) was performed using pyrophosphate (final concentration, 5 mM) and ultrasound treatment (three 1-min treatments using a Branson Sonifier 2200; 60 W)87. Then, samples were diluted from 50- to 100-fold with sterile water filtered onto 0.2-μm pore-size filters (Anodisc filters; black-stained polycarbonate). The filters were stained using SYBR Green I (10,000× in anhydrous dimethyl sulfoxide, Molecular Probes-Invitrogen) diluted 1:20 in prefiltered TE buffer (pH 7.5) and incubated in the dark for 20 min; a drop (20 µl) of antifade solution (composed of 50% 6.7 mmol L−1 phosphate buffer at pH 7.8 and 50% glycerol with the addition of 0.5% ascorbic acid) was laid both on a glass slide and on the filter mounted on it. Prokaryotic counts were performed under epifluorescence microscopy (magnification, ×1000; Zeiss filter set #09, 488009-9901-000, excitation BP 450–490 nm, beam splitter FT 515, emission LP 520), by examining at least 20 fields per slide and counting at least 400 cells per filter.For the determination of prokaryote abundance in seawater surrounding corals, three replicates of 10 ml of seawater were collected from each tank. Total prokaryotic abundance was determined according to Danovaro87. Samples were filtered onto 0.2-μm pore-size filters (Anodisc black-stained polycarbonate filters, Whatman) into a funnel with vacuum pressure no greater than 20 kPa (or 150 mmHg) to avoid cell damage. When the sample had passed through, filters were stained with 20 µl of SYBR Green I (10,000× in anhydrous dimethyl sulfoxide, Molecular Probes-Invitrogen) diluted 1:20 in prefiltered TE buffer (pH 7.5) and incubated in the dark for 20 min. Then, to remove the excess stain, filters were washed three times using 3 ml of Milli-Q water; a drop (20 µl) of antifading solution (composed of 50% 6.7 mmol L−1 phosphate buffer at pH 7.8 and 50% glycerol with the addition of 0.5% ascorbic acid) was laid both on a glass slide and on the filter mounted on it. Prokaryotic counts were carried out as described above.Microbiome of corals exposed to microplasticsThe coral microbiome was analysed immediately before the start of the experiment (before the addition of microplastics) and after 10 days of the experiment, both in replicated coral branches exposed to microplastics and in unexposed corals (Control t10). For the analysis of the microbiome, ca. 20 mg of tissue (wet weight) from one coral branch randomly collected from two tanks of each treatment and control was scraped from the skeleton by using surgical disposable scalpels (Braun) and DNA extraction was performed using the QIAGEN DNeasy Blood & Tissue Kit. Briefly, samples were digested with proteinase K at 56 °C overnight or until the tissue was completely lysed, then samples were processed following the manufacturer’s protocol. Finally, samples were held at −20 °C before PCR amplification and sequencing. The molecular size of the DNA extracts was analysed by agarose gel electrophoresis (1%) and the amount and purity of DNA were determined by Nanodrop spectrophotometer (ND-1000). For PCR amplification of the 16S V3 region, the Bacteria-specific primer pair 805R/341F was chosen with Illumina-specific adapters and barcodes. Sequencing was performed on an Illumina MiSeq platform by LGC Genomics GmbH (Berlin, Germany).Raw sequencing paired-end reads were first joined using the bbmerge tool from the BBMap suite88 in a two-step process: reads that did not merge in a first step were quality-trimmed to remove low-quality bases (Q  More

DataPathosystemWMV is widespread in cucurbit crops, mostly in temperate and Mediterranean climatic regions of the world16. WMV has a wide host range including some legumes, orchids and many weeds that can be alternative hosts16. Like other potyviruses, it is non-persistently transmitted by at least 30 aphid species16. In temperate regions, WMV causes summer epidemics on cucurbit crops, and it can overwinter in several common non-cucurbit weeds when no crops are present16,34. WMV has been common in France for more than 40 years, causing mosaics on leaves and fruits in melon, but mostly mild symptoms on zucchini squash. Since 2000, new symptoms were observed in southeastern France on zucchini squash: leaf deformations and mosaics, as well as fruit discoloration and deformations that made them unmarketable. This new agronomic problem was correlated to the introduction of new molecular groups of WMV strains. At least four new groups have emerged since 2000 and they have rapidly replaced the native “classical” strains, causing important problems for the producers35. These new groups, hereafter “emerging strains” (ES) are significantly more related molecularly to worldwide strains than to any other isolates from the French populations36. As emphasised in35, this supports that the new group of emerging strains has arisen through introductions, mostly from Southeastern Asia, rather than through local differentiation.In this study, we focus on the pathosystem corresponding to a classical strain (CS) and four emerging strains (ESk, (k = 1, ldots ,4)) of WMV and their cucurbit hosts.Study area and samplingThe study area, located in Southeastern France, is included in a rectangle of about 25,000 km2 and is bounded on the South by the Mediterranean Sea. Between 2004 and 2008, the presence of WMV had been monitored in collaboration with farmers, farm advisers and seed companies. Each year, cultivated host plants were collected in different fields and at different dates between May 1st and September 30th. In total, more than two thousand plant samples were collected over the entire study area. All plant samples were analyzed in the INRAE Plant Pathology Unit to confirm the presence of WMV and determine the molecular type of the virus strain causing the infection (see35 for detail on field and laboratory protocols). All infected host plants were cucurbits, mostly melon and different squashes (e.g., zucchini, pumpkins).Observations In the absence of individual geographic coordinates, all infected host plants were attributed to the centroid of the municipality (French administrative unit, median size about 10 km2) where they have been collected. Then for one date, one observation corresponded to a municipality in which at least one infected host plant was sampled. Table 1 summarizes for each year, the number of observations (i.e. number of municipalities), the number of infected plants sampled and the proportion of each WMV strain (CS, and ES1 to ES4) found in the infected host plants. Errors in assignment of virus samples to the CS or ES strains was negligible because of the large genetic distance separating them: 5 to 10% nucleotide divergence both in the fragment used in the study and in complete genomes35, also precluding the possibility of local jumps between groups by accumulation of mutations.Table 1 Number of observations and corresponding proportions of classical and emerging strains.Full size tableLandscapeTo approximate the density of WMV host plants over the study area, we used 2006 land use data (i.e. BD Ocsol 2006 PACA and LR) produced by the CRIGE PACA (http://www.crige-paca.org/) and the Association SIG-LR (http://www.siglr.org/lassociation/la-structure.html). Based on satellite images, land use is determined at a spatial resolution of 1/50,000 using an improved three-level hierarchical typology derived from the European Corine Land Cover nomenclature. Here we used the third hierarchical level of the BD Ocsol typology (i.e. 42 land use classes) to classify the entire study area in three habitats: (1) WMV-susceptible crops, (2) habitats unfavorable to WMV host plants (e.g. forests, industrial and commercial units…) and, (3) non-terrestrial habitat (i.e. water). The proportion of WMV-susceptible crops was then computed within all cells of a raster covering the entire study area, with a spatial resolution of (1.4 times 1.4) km2. These proportions were used to approximate host plant density (zleft( {varvec{x}} right)), which was normalized, so that (zleft( {varvec{x}} right) = 0) corresponds to an absence of host plants and (zleft( {varvec{x}} right) = 1) to the maximum density of host plants (Fig. 1).Figure 1Approximated density (zleft( x right)) of the host plants in the study area. The density was normalized, so that (zleft( x right) = zleft( {x_{1} ,x_{2} } right) = 0) corresponds to an absence of cucurbit plants and (zleft( x right) = 1) to the maximum density. The axes (x_{1}) and (x_{2}) correspond to Lambert93 coordinates (in km). The white regions are non-terrestrial habitats (water). Land use data were not available in the gray regions; the host plant density was then computed by interpolation.Full size imageMechanistic-statistical modelThe general modeling strategy is based on a mechanistic-statistical approach12,22,33. This type of approach combines a mechanistic model describing the dynamics under investigation with a probabilistic model conditional on the dynamics, describing how the measurements have been collected. This method that has already proved its theoretical effectiveness in determining dispersal parameters using simulated genetic data12 aims at bridging the gap between the data and the model for the determination of virus dynamics.Here, the mechanistic part of the model describes the spatio-temporal dynamics of the virus strains, given the model parameters (demographic parameters, introduction dates/sites). This allows us to compute the expected proportions of the five types of virus strains (CS and ES1 to ES4) at each date and site of observation. The probabilistic part of the mechanistic-statistical model describes the conditional distribution of the observed proportions of the virus strains, given the expected proportions. Using this approach, it is then possible to derive a numerically tractable formula for the likelihood function associated with the model parameters.Population dynamicsThe model is segmented into two stages: (1) the intra-annual stage describes the dispersal and growth of the five virus strains during the summer epidemics on cucurbit crops, and the competition between them, during a period ranging from May 1st (noted (t = 0)) to September 30 (noted (t = t_{f}), (t_{f} = 153) days); (2) the inter-annual stage describes the winter decay of the different strains when no crops are present and the virus overwinters in weeds. We denote by (c^{n} left( {t,{varvec{x}}} right)) and (e_{k}^{n} left( {t,{varvec{x}}} right)) the densities of classical strain (CS) and emerging strains (ESk, (k = 1, ldots ,4)), at position ({varvec{x}}) and at time (t) of year (n.)Dynamics of the classical strain before the first introduction eventsBefore the introduction of the first emerging strain, the intra-annual dynamics of the population CS is described by a standard diffusion model with logistic growth: (partial_{t} c^{n} = D{Delta }c^{n} + rc^{n} left( {zleft( {varvec{x}} right) – c^{n} } right)). Here, ({Delta }) is the Laplace 2D diffusion operator (sum of the second derivatives with respect to coordinate). This operator describes uncorrelated random walk movements of the viruses, with the coefficient (D) measuring the mobility of the viruses (e.g.,26,37). The term (r zleft( {varvec{x}} right)) is the intrinsic growth rate (i.e., growth rate in the absence of competition) of the population, which depends on the density of host plants (zleft( {varvec{x}} right)) and on a coefficient (r) (intrinsic growth rate at maximum host density). Under these assumptions, the carrying capacity at a position ({varvec{x}}) is equal to (zleft( {varvec{x}} right)), which means that the population densities are expressed in units of the maximum host population density. The model is initialized by setting (c^{1980} left( {0,{varvec{x}}} right) = (1 – m_{c} ) zleft( {varvec{x}} right)), where (m_{c}) is the winter decay rate of the CS (see the description of the inter-annual stage below). In other terms, we assume that the CS density is at the carrying capacity in 1979, i.e., 5 years after its first detection and 20 years before the first detections of ESs38.Introduction eventsThe ESs are introduced during years noted (n_{k} ge 1981), at the beginning of the intra-annual stage (other dates of introduction within the intra-annual stage would lead—at most—to a one-year lag in the dynamics). Their densities are (0) before introduction: (e_{k}^{n} = 0) for (n < n_{k}). Once introduced, the initial density of any ES is assumed to be 1/10th of the carrying capacity at the introduction point (other values have been tested without much effect, see Supplementary Fig. S1), with a decreasing density as the distance from this point increases:$$e_{k}^{{n_{k} }} left( {0,x} right) = frac{{zleft( {varvec{x}} right)}}{10}exp left( { - frac{|{{varvec{x}} - {varvec{X}}_{{varvec{k}}}|^{2} }}{{2sigma^{2} }}} right),$$where ({varvec{X}}_{{varvec{k}}}) is the location of introduction of the strain (k.) In our computations, we took (sigma = 5) km for the standard deviation.Intra-annual dynamics after the first introduction eventIntra-annual dynamics were described by a neutral competition model with diffusion (properties of this model have been analyzed in [54]):$$left{ {begin{array}{*{20}c} {partial_{t} c^{n} left( {t,x} right) = DDelta c^{n} + rc^{n} left( {zleft( {varvec{x}} right) - c^{n} - mathop sum limits_{i = 1}^{4} e_{i}^{n} left( {t,{varvec{x}}} right)} right)} \ {partial_{t} e_{k}^{n} left( {t,x} right) = DDelta e_{k}^{n} + re_{k}^{n} left( {zleft( {varvec{x}} right) - c^{n} - mathop sum limits_{i = 1}^{4} e_{i}^{n} left( {t,{varvec{x}}} right)} right)} \ end{array} } right.,$$for (t = 0 ldots t_{f}) and for all introduced emerging strains, i.e. all (k) such that (n ge n_{k} .) We assume reflecting boundary conditions, meaning that the population flows vanish at the boundary of the study area, due to truly reflecting boundaries (e.g., sea coast in the southern part of the site) or symmetric inward and outward fluxes26. In addition, in order to limit the number of unknown parameters, and in the absence of precise knowledge on the differences between the strains, we assume here that the diffusion, competition and growth coefficients are common to all the strains during the intra-annual stage (see the discussion for more details on this assumption).Inter-annual dynamicsThe population densities at time (t = 0) of year (n) are connected with those of year (n - 1,) at time (t = t_{f} ,) through the following formulas:$$left{ {begin{array}{*{20}c} {c^{n} left( {0,{varvec{x}}} right) = left( {1 - m_{c} } right)c^{n - 1} left( {t_{f} ,{varvec{x}}} right) hbox{ for } n ge 1981} \ {e_{k}^{n} left( {0,{varvec{x}}} right) = left( {1 - m_{e} } right)e_{k}^{n - 1} left( {t_{f} ,{varvec{x}}} right) hbox{ for }n ge n_{k} + 1} \ end{array} } right.$$with (m_{c}) and (m_{e}) the winter decay rates of the CS and ESs strains (note that (m_{e}) is common to all of the ESs). Estimation of CS and ES decay rates provides an assessment of the competitive advantage of one type of strain vs the other.Numerical computationsThe intra-annual dynamics were solved using Comsol Multiphysics time-dependent solver, which is based on a finite element method (FEM). The triangular mesh which was used for our computations is available as Supplementary Fig. S2.Probabilistic model for the observation processDuring the years (n = 2004, ldots ,2008), (I_{n}) observations were made (see Section Observations above and Table 1). They consist in counting data, that we denote by (C_{i}) and (E_{k,i}) for (k = 1, ldots ,4) and (i = 1, ldots ,I_{n}), corresponding to the number of samples infected by the CS and ESs strains, respectively, at each date of observation and location (left( {t_{i} ,{varvec{x}}_{i} } right)). Note that these dates and locations depend on the year of observation (n). More generally, the above quantities should be noted (C_{i}^{n} , E_{k,i}^{n} , t_{i}^{n} , {varvec{x}}_{i}^{n} ;) for simplicity, the index (n) is omitted in the sequel, unless necessary.We denote by (V_{i} = C_{i} + mathop sum nolimits_{k = 1}^{4} E_{k,i}) the total number of infected samples observed at (left( {t_{i} ,{varvec{x}}_{i} } right)). The conditional distribution of the vector (left( {C_{i} ,E_{1,i} ,E_{2,i} ,E_{3,i} ,E_{4,i} } right)), given (V_{i}) can be described by a multinomial distribution ({mathcal{M}}left( {V_{i} ,{varvec{p}}_{i} } right)) with ({varvec{p}}_{i} = left( {p_{i}^{c} ,p_{i}^{{e_{1} }} ,p_{i}^{{e_{2} }} ,p_{i}^{{e_{3} }} ,p_{i}^{{e_{4} }} } right)) the vector of the respective proportions of each strain in the population at (left( {t_{i} ,{varvec{x}}_{i} } right)). We chose to work conditionally to (V_{i}) because the sample sizes are not related to the density of WMV.Statistical inferenceUnknown parametersWe denote by ({{varvec{Theta}}}) the vector of unknown parameters: the diffusion coefficient (D,) the intrinsic growth rate at maximum host density (r), the winter decay rates ((m_{c} , m_{e} )) and the locations ((x_{k} in {mathbb{R}}^{2})) and years ((n_{k})) of introduction, for (k = 1, ldots ,4.) Thus ({{varvec{Theta}}} in {mathbb{R}}^{16} .)Computation of a likelihoodGiven the set of parameters ({{varvec{Theta}}}), the densities (c^{n} left( {t,{varvec{x}}|{{varvec{Theta}}}} right)) and (e_{k}^{n} left( {t,{varvec{x}}|{{varvec{Theta}}}} right)) can be computed explicitly with the mechanistic model for population dynamics. Thus, at a given year (n), at (left( {t_{i} ,x_{i} } right)), the parameter ({varvec{p}}_{i}) of the multinomial distribution ({mathcal{M}}left( {V_{i} ,{varvec{p}}_{i} } right)) writes:$$p_{i}^{c} left( {{varvec{Theta}}} right) = frac{{c^{n} left( {t_{i} ,{varvec{x}}_{i} |{{varvec{Theta}}}} right)}}{{c^{n} left( {t_{i} ,{varvec{x}}_{i} |{{varvec{Theta}}}} right) + mathop sum nolimits_{i = 1}^{4} e_{i}^{n} left( {t_{i} ,{varvec{x}}_{i} {|}{{varvec{Theta}}}} right)}}, p_{i}^{{e_{k} }} left( {{varvec{Theta}}} right) = frac{{e_{k}^{n} left( {t_{i} ,{varvec{x}}_{i} |{{varvec{Theta}}}} right)}}{{c^{n} left( {t_{i} ,{varvec{x}}_{i} |{{varvec{Theta}}}} right) + mathop sum nolimits_{i = 1}^{4} e_{i}^{n} (t_{i} ,{varvec{x}}_{i} |{{varvec{Theta}}})}}.$$The probability (P(C_{i} ,E_{1,i} ,E_{2,i} ,E_{3,i} ,E_{4,i} |{{varvec{Theta}}},{text{V}}_{{text{i}}} )) of the observed outcome (C_{i} ,E_{1,i} ,E_{2,i} ,E_{3,i} ,E_{4,i}) is then$$Pleft( {C_{i} ,E_{1,i} ,E_{2,i} ,E_{3,i} ,E_{4,i} {|}{{varvec{Theta}}},{text{V}}_{{text{i}}} } right) = frac{{left( {V_{i} } right)!}}{{C_{i} ! mathop prod nolimits_{k = 1}^{4} E_{k,i} !}}left( {p_{i}^{c} left( {{varvec{Theta}}} right)} right)^{{C_{i} }} mathop prod limits_{k = 1}^{4} (p_{i}^{{e_{k} }} left( {{varvec{Theta}}} right))^{{E_{k,i} }} .$$Assuming that the observations during each year and at each date/location are independent from each other conditionally on the virus strain proportions, we get the following formula for the likelihood:$${mathcal{L}}left( {{varvec{Theta}}} right) = mathop prod limits_{n = 2004}^{2008} mathop prod limits_{{i = 1, ldots , I_{n} }} Pleft( {C_{i} ,E_{1,i} ,E_{2,i} ,E_{3,i} ,E_{4,i} {|}{{varvec{Theta}}},{text{V}}_{{text{i}}} } right).$$A priori constraints on the parameters By definition and for biological reasons, the parameter vector ({{varvec{Theta}}}) satisfies some constraints. First, (D in left( {10^{ - 4} ,10} right){text{ km}}^{2} /{text{day}}), (r in left( {0.1,1} right) {text{day}}^{ - 1} ,) and (m_{c} , m_{e} in left{ {0,0.1,0.2, ldots ,0.9} right},) (see Supplementary Note S7 for a biological interpretation of these values). Second, we assumed that the locations of introductions ({varvec{X}}_{{varvec{k}}}) belong to the study area. To facilitate the estimation procedure, the points ({varvec{X}}_{{varvec{k}}}) were searched in a regular grid with 20 points (see Supplementary Fig. S3), and the dates of introduction (n_{k}) were searched in (left{ {1985,1990,1995,2000} right}.)Inference procedureDue to the important computation time (4 min in average for one simulation of the model on an Intel(R) Core(R) CPU i7-4790 @ 3.60 GHz), we were not able to compute an a posteriori distribution of the parameters in a Bayesian framework. Rather, we used a simulated annealing algorithm to compute the maximum likelihood estimate (MLE), i.e., the parameter ({{varvec{Theta}}}^{*}) which leads to the highest log-likelihood. This is an iterative algorithm, which constructs a sequence (({{varvec{Theta}}}_{j} )_{j ge 1}) converging in probability towards a MLE. It is based on an acceptance-rejection procedure, where the acceptance rate depends on the current iteration (j) through a "cooling rate" ((alpha )). Empirically, a good trade-off between quality of optimization and time required for computation (number of iterations) is obtained with exponential cooling rates of the type (T_{0} alpha^{j}) with (0 < alpha < 1) and some constant (T_{0} gg 1) (this cooling schedule was first proposed in= 39 = 39). Too rapid cooling ((alpha ll 1)) results in a system frozen into a state far from the optimal one, whereas too slow cooling ((alpha approx 1)) leads to important computation times due to very slow convergence. Here, we ran (6) parallel sequences with cooling rates (alpha in left{ {0.995,0.999,0.9995} right}). For this type of algorithm, there are no general rules for the choice of the stopping criterion [HenJac03], which should be heuristically adapted to the considered optimization problem. Here, our stopping criterion was that ({{varvec{Theta}}}_{j}) remained unchanged during 500 iterations. The computations took about 100 days (CPU time).Confidence intervals and goodness-of-fitTo assess the model’s goodness-of-fit, 95% confidence regions were computed for the observations (left( {C_{i} ,E_{1,i} ,E_{2,i} ,E_{3,i} ,E_{4,i} } right)) at each date/location (left( {t_{i} ,{varvec{x}}_{i} } right),) and for each year of observation. The confidence regions were computed by assessing the probability of each possible outcome of the observation process, at each date/location, based on the computed proportions ({varvec{p}}_{i} = left( {p_{i}^{c} ,p_{i}^{{e_{1} }} ,p_{i}^{{e_{2} }} ,p_{i}^{{e_{3} }} ,p_{i}^{{e_{4} }} } right)), corresponding to the output of the mechanistic model using the MLE ({{varvec{Theta}}}^{user2{*}}) and given the total number of infected samples (V_{i}). Then, we checked if the observations belonged to the 95% most probable outcomes. More

1.Capece MC, Clark E, Saleh JK, Halford D, Heinl N, Hoskins S, et al. Polyextremophiles and the constraints for terrestrial habitability. In: Seckbach J, Oren A, Stan-Lotter H, editors. Polyextremophiles. Life under muliple forms of stress. Dordrecht, Neaderlands: Springer; 2013. p. 3–60.2.Harrison JP, Gheeraert N, Tsigelnitskiy D, Cockell CS. The limits for life under multiple extremes. Trends Microbiol. 2013;21:204–12.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
3.Oger PM, Jebbar M. The many ways of coping with pressure. Res Microbiol. 2010;161:799–809.PubMed 
Article 
PubMed Central 

Google Scholar 
4.Whitman WB, Coleman DC, Wiebe WJ. Prokaryotes: the unseen majority. Proc Natl Acad Sci USA. 1998;95:6578–83.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
5.Bartlett DH. Pressure effects on in vivo microbial processes. Biochim Biophys Acta. 2002;1595:367–81.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
6.Aertsen A, Meersman F, Hendrickx ME, Vogel RF, Michiels CW. Biotechnology under high pressure: applications and implications. Trends Biotechnol. 2009;27:434–41.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
7.Jannasch HW, Taylor CD. Deep-sea microbiology. Ann Rev Microbiol. 1984;38:487–514.CAS 
Article 

Google Scholar 
8.Fang J, Zhang L, Bazylinski DA. Deep-sea piezosphere and piezophiles: geomicrobiology and biogeochemistry. Trends Microbiol. 2010;18:413–22.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
9.Yayanos AA. Microbiology to 10,500 meters in the deep sea. Ann Rev Microbiol. 1995;49:777–805.CAS 
Article 

Google Scholar 
10.Eloe EA, Lauro FM, Vogel RF, Bartlett DH. The deep-sea bacterium Photobacterium profundum SS9 utilizes separate flagellar systems for swimming and swarming under high-pressure conditions. Appl Environ Microbiol. 2008;74:6298–305.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
11.Horikoshi K, Bull AT Prologue: Definition, categories, distribution, origin and evolution, pioneering studies, and emerging fields of extremophiles. In: Horikoshi K, editor. Extremophiles handbook. Tokyo, Japan: Springer; 2011. p. 3–18.12.Holt RD. Bringing the Hutchinsonian niche into the 21st century: ecological and evolutionary perspectives. Proc Natl Acad Sci USA. 2009;106:19659–65.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
13.Talley LD, Pickard GL, Emery WJ, Swift JH. Typical distributions of water characteristics. In: Descriptive physical oceanography, 6th ed. London, UK: Elsevier; 2011. p. 67–110.14.Jebbar M, Franzetti B, Girard E, Oger P. Microbial diversity and adaptation to high hydrostatic pressure in deep-sea hydrothermal vents prokaryotes. Extremophiles. 2015;19:721–40.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
15.Berhardt G, Jaenicke R, Ludemann H-D, Konig H, Stetter KO. High pressure enhances the growth rate of the thermophilic archaebacterium Methanococcus thermolithotrophicus without extending its temperature range. Appl Environ Microbiol. 1998;54:1258–61.Article 

Google Scholar 
16.Scoma A, Garrido-Amador P, Nielsen SD, Roy H, Kjeldsen KU. The polyextremophilic bacterium Clostridium paradoxum attains piezophilic traits by modulating its energy metabolism and cell membrane composition. Appl Environ Microbiol. 2019;85:e00802–19.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
17.Wiegel J. Temperature spans for growth: hypothesis and discussion. FEMS Microbiol Rev. 1990;75:155–70.Article 

Google Scholar 
18.Morita RY. Psychrophilic bacteria. Bacteriol Rev. 1975;39:144–67.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
19.Zeikus JG. Thermophilic Bacteria—Ecology. Physiol Technol Enz Microb Technol. 1979;1:243–52.CAS 
Article 

Google Scholar 
20.Yayanos AA. Evolutional and ecological implications of the properties of deep-sea barophilic bacteria. Proc Natl Acad Sci USA. 1986;83:9542–6.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
21.Jannasch HW, Wirsen CO. Variability of pressure adaptation in deep sea bacteria. Arch Microbiol. 1984;139:281–8.Article 

Google Scholar 
22.Yayanos AA, Chastain R. The influence of nutrition on the physiology of piezophilic bacteria. In: Bell CR, Brylinsky M, Johnson-Green P, Eds. Proceedings of the 8th International Symposium on Microbial Ecology. Halifax, NS, Canada: Atlantic Canada Society for Microbial Ecology; 6; 1999.23.Matsumura P, Keller DM, Marquis RE. Restricted pH ranges and reduced yields for bacterial growth under pressure. Microb Ecol. 1974;1:176–89.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
24.Oren A. Bioenergetic aspects of halophilism. Microbiol Mol Biol Rev. 1999;63:334–48.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
25.Yayanos AA, Dietz AS, Van, Boxtel R. Dependence of reproduction rate on pressure as a hallmark of deep-sea bacteria. Appl Environ Microbiol. 1982;44:1356–61.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
26.Deming JW, Hada H, Colwell RR, Luehrsen KR, Fox GE. The ribonucleotide sequence of 5S rRNA from two strains of deep-sea barophilic bacteria. J Gen Microbiol. 1984;130:1911–20.CAS 
PubMed 
PubMed Central 

Google Scholar 
27.Lauro FM, Chastain RA, Blankenship LE, Yayanos AA, Bartlett DH. The unique 16S rRNA genes of piezophiles reflect both phylogeny and adaptation. Appl Environ Microbiol. 2007;73:838–45.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
28.Marteinsson VT, Birrien J-L-, Reysenbach A-L, Vernet M, Marie D, Gambacorta A, et al. Thermococcus barophilus sp. nov., a new barophilic and hyperthermophilic archaeon isolated under high hydrostatic pressure from a deep-sea hydrothermal vent. Int J Syst Bacteriol. 1999;49:351–9.PubMed 
Article 
PubMed Central 

Google Scholar 
29.Alain K. Marinitoga piezophila sp. nov., a rod-shaped, thermo-piezophilic bacterium isolated under high hydrostatic pressure from a deep-sea hydrothermal vent. Int J Sys Evol Microbiol. 2002;52:1331–9.CAS 

Google Scholar 
30.Canganella F, Gonzalez JM, Yanagibayashi M, Kato C, Horikoshi K. Pressure and temperature effects on growth and viability of the hyperthermophilic archaeon Thermococcus peptonophilus. Arch Microbiol. 1997;168:1–7.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
31.Canganella F, Gambacorta A, Kato C, Horikoshi K. Effects of hydrostatic pressure and temperature on physiological traits of Thermococcus guaymasensis and Thermococcus aggregans growing on starch. Microbiol Res. 2000;154:297–306.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
32.Tamburini C, Boutrif M, Garel M, Colwell RR, Deming JW. Prokaryotic responses to hydrostatic pressure in the ocean-a review. Environ Microbiol. 2013;15:1262–74.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
33.Nogi Y, Masui N, Kato C. Photobacterium profundum sp. nov., a new, moderately barophilic bacterial species isolated from a deep-sea sediment. Extremophiles. 1998;2:1–7.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
34.Arakawa S, Nogi Y, Sato T, Yoshida Y, Usami R, Kato C. Diversity of piezophilic microorganisms in the closed ocean Japan Sea. Biosci Biotechnol Biochem. 2006;70:749–52.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
35.Xu Y, Nogi Y, Kato C, Liang Z, Ruger H-J, De Kegel D, et al. Moritella profunda sp. nov. and Moritella abyssi sp. nov., two psychropiezophilic organisms isolated from deep Atlantic sediments. Int J Syst Evol Microbiol. 2003;53:533–8.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
36.Sekiguchi T, Sato T, Enoki M, Kanehiro H, Kato C. Procedure for isolation of the plastic degrading piezophilic bacteria from deep-sea environments. J Jap Soc Extremophil. 2010a;9:25–30.Article 

Google Scholar 
37.Sekiguchi T, Sato T, Enoki M, Kanehiro H, Uematsu K, Kato C. Isolation and characterization of biodegradable plastic degrading bacteria from deep-sea environments. JAMSTEC Rep. Res Dev. 2010b;11:33–41.
Google Scholar 
38.Nogi Y, Kato C, Horikoshi K. Psychromonas kaikoae sp. nov., a novel piezophilic bacterium from the deepest cold-seep sediments in the Japan Trench. Int J Syst Evol Microbiol. 2002;52:1527–32.CAS 
PubMed 
PubMed Central 

Google Scholar 
39.Yayanos AA, Dietz AS, van Boxtel R. Isolation of a deep-sea barophilic bacterium and some of its growth characteristics. Science. 1979;205:808–10.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
40.Nogi Y, Hosoya S, Kato C, Horikoshi K. Colwellia piezophila sp. nov., a novel piezophilic species from deep-sea sediments of the Japan Trench. Int J Syst Evol Microbiol. 2004;54:1627–31.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
41.Kato C, Sato T, Horikoshi K. Isolation and properties of barophilic and barotolerant bacteria from deep-sea mud samples. Biodiv Cons. 1995;4:1–9.Article 

Google Scholar 
42.Kato C, Inoue A, Horikoshi K. Isolating and characterizing deep-sea marinemicroorganisms. Tibtech. 1996;14:6–12.CAS 
Article 

Google Scholar 
43.Nogi Y, Kato C. Taxonomic studies of extremely barophilic bacteria isolated from the Mariana Trench and description of Moritella yayanosii sp. nov., a new barophilic bacterial isolate. Extremophiles. 1999;3:71–7.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
44.Deming JW, Somers LK, Straube WL, Swartz DG, Macdonell MT. Isolation of an Obligately Barophilic Bacterium and Description of a New Genus, Colwellia gen. nov. Systematic and Applied Microbiology. 1988;10:152–60.Article 

Google Scholar 
45.Kusube M, Kyaw TS, Tanikawa K, Chastain RA, Hardy KM, Cameron J, et al. Colwellia marinimaniae sp. nov., a hyperpiezophilic species isolated from an amphipod within the Challenger Deep, Mariana Trench. Int J Syst Evol Microbiol. 2017;67:824–31.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
46.Cao J, Lai Q, Liu P, Wei Y, Wang L, Liu R, et al. Salinimonas sediminis sp. nov., a piezophilic bacterium isolated from a deep-sea sediment sample from the New Britain Trench. Int J Syst Evol Microbiol. 2018;68:3766–71.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
47.Liu P, Ding W, Lai Q, Liu R, Wei Y, Wang L, et al. Physiological and genomic features of Paraoceanicella profunda gen. nov., sp. nov., a novel piezophile isolated from deep seawater of the Mariana Trench. MicrobiologyOpen. 2019;00:e966.
Google Scholar 
48.Quéméneur M, Erauso G, Frouin E, Zeghal E, Vandecasteele C, Ollivier B, et al. Hydrostatic Pressure Helps to Cultivate an Original Anaerobic Bacterium From the Atlantis Massif Subseafloor (IODP Expedition 357): Petrocella atlantisensis gen. nov. sp. nov. Front Microbiol. 2019;10:1497.PubMed 
PubMed Central 
Article 

Google Scholar 
49.Xiao X, Wang P, Zeng X, Bartlett DH, Wang F. Shewanella psychrophila sp. nov. and Shewanella piezotolerans sp. nov., isolated from west Pacific deep-sea sediment. Int J Syst Evol Microbiol. 2007;57:60–5.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
50.Alazard D, Dukan S, Urios A, Verhe F, Bouabida N, Morel F, et al. Desulfovibrio hydrothermalis sp. nov., a novel sulfate-reducing bacterium isolated from hydrothermal vents. Int J Syst Evol Microbiol. 2003;53:173–8.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
51.Pathom-Aree W, Nogi Y, Sutcliffe IC, Ward AC, Horikoshi K, Bull AT, et al. Dermacoccus abyssi sp. nov., a piezotolerant actinomycete isolated from the Mariana Trench. Int J Syst Evol Microbiol. 2006;56:1233–7.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
52.Takai K, Miyazaki M, Hirayama H, Nakagawa S, Querellou J, Godfroy A. Isolation and physiological characterization of two novel, piezophilic, thermophilic chemolithoautotrophs from a deep-sea hydrothermal vent chimney. Environ Microbiol. 2009;11(8):1983–97.PubMed 
Article 
PubMed Central 

Google Scholar 
53.Erauso G, Reysenbach A-L, Godfroy A, Meunier J-R, Crump B, Partensky F, et al. Pyrococcus abyssi sp. nov., a new hyperthermophilic archaeon isolated from a deep-sea hydrothermal vent. Arch Microbiol. 1993;160:338–49.CAS 
Article 

Google Scholar 
54.Li Y, Mandelco L, Wiegel J. Isolation and Characterization of a Moderately Thermophilic Anaerobic Alkaliphile. Clostridium paradoxum sp. nov. Int J Sys Bacteriol. 1993;43:450–60.Article 

Google Scholar 
55.Zhao W, Zeng X, Xiao X. Thermococcus eurythermalis sp. nov., a conditional piezophilic, hyperthermophilic archaeon with a wide temperature range for growth, isolated from an oil-immersed chimney in the Guaymas Basin. Int J Sys Evol Microbiol. 2015;65:30–5.CAS 
Article 

Google Scholar 
56.Takai K, Nakamura K, Toki T, Tsunogai U, Miyazaki M, Miyazaki J, et al. Cell proliferation at 122 degrees C and isotopically heavy CH4 production by a hyperthermophilic methanogen under high-pressure cultivation. Proc Natl Acad Sci U S A. 2008;105:10949–54.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
57.González JM, Kato C, Horikoshi K. Thermococcus peptonophilus sp. nov., a fast-growing, extremely thermophilic archaebacterium isolated from deep-sea hydrothermal vents. Arch Microbiol. 1995;164:159–64.PubMed 
Article 
PubMed Central 

Google Scholar 
58.Jones WJ, Leigh JA, Mayer F, Woese CR, Wolfe RS. Methanococcusjannaschii sp. nov., an extremely thermophilic methanogen from a submarine hydrothermal vent. Arch Microbiol. 1983;136:254–61.CAS 
Article 

Google Scholar  More

1.Prosser, J. I. & Nicol, G. W. Relative contributions of archaea and bacteria to aerobic ammonia oxidation in the environment. Environ. Microbiol. 10, 2931–2941 (2008).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
2.Bernhard, A. E. & Bollmann, A. Estuarine nitrifiers: new players, patterns and processes. Estuar. Coast. Shelf Sci. 88, 1–11 (2010).CAS 
Article 

Google Scholar 
3.Martiny, J. B. H., Eisen, J., Penn, K., Allison, S. D. & Horner-Devine, M. C. Drivers of bacterial beta-diversity depend on spatial scale. Proc. Natl Acad. Sci. USA 108, 7850–7854 (2011).4.Nelson, M. B., Martiny, A. C. & Martiny, J. B. H. Global biogeography of microbial nitrogen-cycling traits in soil. Proc. Natl Acad. Sci. USA 113, 8033–8040 (2016).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
5.Daims, H. et al. Complete nitrification by Nitrospira bacteria. Nature 528, 504–509 (2015).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
6.Marton, J. M., Roberts, B. J., Bernhard, A. E. & Giblin, A. E. Spatial and temporal variability of nitrification potential and ammonia-oxidizer abundances in Louisiana salt marshes. Estuaries Coast. 38, 1824–1837 (2015).CAS 
Article 

Google Scholar 
7.Martens-Habbena, W., Berube, P. M., Urakawa, H., de la Torre, J. R. & Stahl, D. A. Ammonia oxidation kinetics determine niche separation of nitrifying archaea and bacteria. Nature 461, 976–981 (2009).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
8.Dimitri Kits, K. et al. Kinetic analysis of a complete nitrifier reveals an oligotrophic lifestyle. Nature 549, 269–272 (2017).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
9.Hink, L., Nicol, G. W. & Prosser, J. I. Archaea produce lower yields of N2O than bacteria during aerobic ammonia oxidation in soil. Environ. Microbiol. 19, 4829–4837 (2017).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
10.Bernhard, A. E., Donn, T., Giblin, A. E. & Stahl, D. A. Loss of diversity of ammonia-oxidizing bacteria correlates with increasing salinity in an estuary system. Environ. Microbiol. 7, 1289–1297 (2005).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
11.Moin, N. S., Nelson, K. A., Bush, A. & Bernhard, A. E. Distribution and diversity of archaeal and bacterial ammonia oxidizers in salt marsh sediments. Appl. Environ. Microbiol. 75, 7461–7468 (2009).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
12.Bernhard, A. E. et al. Abundance of ammonia-oxidizing archaea and bacteria along an estuarine salinity gradient in relation to potential nitrification rates. Appl. Environ. Microbiol. 76, 1285–1289 (2010).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
13.Francis, C. A., O’Mullan, G. D. & Ward, B. B. Diversity of ammonia monooxygenase (amoA) genes across environmental gradients in Chesapeake Bay sediments. Geobiology 1, 129–140 (2003).CAS 
Article 

Google Scholar 
14.Ward, B. B. et al. Ammonia-oxidizing bacterial community composition in estuarine and oceanic environments assessed using a functional gene microarray. Environ. Microbiol. 9, 2522–2538 (2007).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
15.Mills, H. J. et al. Characterization of nitrifying, denitrifying, and overall bacterial communities in permeable marine sediments of the northeastern Gulf of Mexico. Appl. Environ. Microbiol. 74, 4440–4453 (2008).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
16.Newell, S. E. et al. A shift in the archaeal nitrifier community in response to natural and anthropogenic disturbances in the northern Gulf of Mexico. Environ. Microbiol. Rep. 6, 106–112 (2014).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
17.Bernhard, A. E., Sheffer, R., Giblin, A. E., Marton, J. M. & Roberts, B. J. Population dynamics and community composition of ammonia oxidizers in salt marshes after the Deepwater Horizon oil spill. Front. Microbiol. 7, 854 (2016).PubMed 
PubMed Central 

Google Scholar 
18.Bernhard, A. E., Chelsky, A., Giblin, A. E. & Roberts, B. J. Influence of local and regional drivers on spatial and temporal variation of ammonia-oxidizing communities in Gulf of Mexico salt marshes. Environ. Microbiol. Rep. 11, 825–834 (2019).CAS 
PubMed 
PubMed Central 

Google Scholar 
19.Nelson, K. A., Moin, N. S. & Bernhard, A. E. Archaeal diversity and the prevalence of Crenarchaeota in salt marsh sediments. Appl. Environ. Microbiol. 75, 4211–4215 (2009).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
20.Peng, X. et al. Differential responses of ammonia-oxidizing archaea and bacteria to long-term fertilization in a New England salt marsh. Front. Microbiol. 3, 445 (2012).PubMed 
PubMed Central 

Google Scholar 
21.Bernhard, A. E., Marshall, D. & Yiannos, L. Increased variability of microbial communities in restored salt marshes nearly 30 years after tidal flow restoration. Estuaries Coast. 35, 1049–1059 (2012).CAS 
Article 

Google Scholar 
22.Marton, J. M. & Roberts, B. J. Spatial variability of phosphorus sorption dynamics in Louisiana salt marshes. J. Geophys. Res. Biogeosci. 119, 451–465 (2014).CAS 
Article 

Google Scholar 
23.Hill, T. D. & Roberts, B. J. Effects of seasonality and environmental gradients on Spartina alterniflora allometry and primary production. Ecol. Evol. 7, 9676–9688 (2017).PubMed 
PubMed Central 
Article 

Google Scholar 
24.Bernhard, A. E., Tucker, J., Giblin, A. E. & Stahl, D. A. Functionally distinct communities of ammonia-oxidizing bacteria along an estuarine salinity gradient. Environ. Microbiol. 9, 1439–1447 (2007).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
25.Schutte, C. A., Marton, J. M., Bernhard, A. E., Giblin, A. E. & Roberts, B. J. No evidence for long-term impacts of oil spill contamination on salt marsh soil nitrogen cycling processes. Estuaries Coast. 43, 865–879 (2020).Article 

Google Scholar 
26.Bernhard, A. E., Dwyer, C., Idrizi, A., Bender, G. & Zwick, R. Long-term impacts of disturbance on nitrogen-cycling bacteria in a New England salt marsh. Front. Microbiol. 6 https://doi.org/10.3389/fmicb.2015.00046 (2015).27.Pjevac, P. et al. AmoA-targeted polymerase chain reaction primers for the specific detection and quantification of comammox Nitrospira in the environment. Front. Microbiol. 8 https://doi.org/10.3389/fmicb.2017.01508 (2017).28.Francis, C. A., Roberts, K. J., Beman, J. M., Santoro, A. E. & Oakley, B. B. Ubiquity and diversity of ammonia-oxidizing archaea in water columns and sediments of the ocean. Proc. Natl Acad. Sci. USA 102, 14683–14688 (2005).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
29.Park, S.-J., Park, B.-J. & Rhee, S.-K. Comparative analysis of archaeal 16S rRNA and amoA genes to estimate the abundance and diversity of ammonia-oxidizing archaea in marine sediments. Extremophiles 12, 605–615 (2008).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
30.Ludwig, W. et al. ARB: a software environment for sequence data. Nucleic Acids Res. 32, 1363–1371 (2004).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
31.Kumar, S., Stecher, G. & Tamura, K. MEGA7: Molecular Evolutionary Genetics Analysis Version 7.0 for bigger datasets. Mol. Biol. Evol. 33, 1870–1874 (2016).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
32.Schloss, P. D. et al. Introducing mothur: open-source, platform-independent, community-supported software for describing and comparing microbial communities. Appl. Environ. Microbiol. 75, 7537–7541 (2009).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
33.R Core Team. R: A Language and Environment for Statistical Computing (R Foundation for Statistical Computing, 2017).34.Turner, R. E., Rabalais, N. N. & Justic, D. Predicting summer hypoxia in the northern Gulf of Mexico: riverine N, P, and Si loading. Mar. Pollut. Bull. 52, 139–148 (2006).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
35.Tian, H. et al. Long-term trajectory of nitrogen loading and delivery from Mississippi river basin to the Gulf of Mexico. Glob. Biogeochem. Cycles 34, 6475 (2020).Article 
CAS 

Google Scholar 
36.Dang, H. et al. Diversity, abundance, and spatial distribution of sedimet ammonia-oxidizing Betaproteobacteria in response to environmental gradients and coastal eutrophication in Jiaozhou Bay, China. Appl. Environ. Microbiol. 76, 4691–4702 (2010).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
37.Sims, A., Zhang, Y., Gajaraj, S., Brown, P. B. & Hu, Z. Toward the development of microbial indicators for wetland assessment. Water Res. 47, 1711–1725 (2013).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
38.Zhang, Q. -F. et al. Impacts of Spartina alterniflora invasion on abundance and composition of ammonia oxidizers in estuarine sediment. J. Soils Sediment. 11, 1020–1031 (2011).Article 

Google Scholar 
39.Jin, T. et al. Diversity and quantity of ammonia-oxidizing archaea and bacteria in sediment of the Pearl River Estuary, China. Appl. Microbiol. Biotechnol. 90, 1137–1145 (2011).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
40.Meinhardt, K. A. et al. Evaluation of revised polymerase chain reaction primers for more inclusive quantification of ammonia-oxidizing archaea and bacteria. Environ. Microbiol. Rep. 7, 354–363 (2015).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
41.Marshall, A. et al. Primer selection influences abundance estimates of ammonia oxidizing archaea in coastal marine sediments. Mar. Environ. Res. 140, 90–95 (2018).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
42.Koops, H. P. & Pommerening-Roser, A. Distribution and ecophysiology of the nitrifying bacteria emphasizing cultured species. FEMS Microbiol. Ecol. 37, 1–9 (2001).CAS 
Article 

Google Scholar 
43.Hillebrand, H. On the generality of the latitudinal diversity gradient. Am. Nat. 163, 192–211 (2004).PubMed 
Article 
PubMed Central 

Google Scholar 
44.Pommier, T. et al. Global patterns of diversity and community structure in marine bacterioplankton. Mol. Ecol. 16, 867–880 (2007).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
45.Fierer, N. & Jackson, R. B. The diversity and biogeography of soil bacterial communities. Proc. Natl Acad. Sci. 103, 626–631 (2006).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
46.Hendershot, J. N., Read, Q. D., Henning, J. A., Sanders, N. J. & Classen, A. T. Consistently inconsistent drivers of microbial diversity and abundance at macroecological scales. Ecology 98, 1757–1763 (2017).PubMed 
Article 
PubMed Central 

Google Scholar 
47.Hitchcock, J. N., Mitrovic, S. M., Kobayashi, T. & Westhorpe, D. P. Responses of estuarine bacterioplankton, phytoplankton and zooplankton to dissolved organic carbon (DOC) and inorganic nutrient additions. Estuaries Coast. 33, 78–91 (2010).CAS 
Article 

Google Scholar 
48.Guo, X. -P. et al. Bacterial community structure in response to environmental impacts in the intertidal sediments along the Yangtze Estuary, China. Mar. Pollut. Bull. 126, 141–149 (2018).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
49.Howarth, R. W. Nutrient limitation of net primary production in marine ecosystems. Annu. Rev. Ecol. 19, 89–110 (1988).Article 

Google Scholar 
50.Murrell, M. C. et al. Evidence that phosphorus limits phytoplankton growth in a Gulf of Mexico estuary: Pensacola Bay, Florida, USA. Bull. Mar. Sci. 70, 155–167 (2002).
Google Scholar 
51.Johnson, M. W., Heck, K. L. Jr & Fourqurean, J. W. Nutrient content of seagrasses and epiphytes in the northern Gulf of Mexico: evidence of phosphorus and nitrogen limitation. Aquat. Bot. 85, 103–111 (2006).CAS 
Article 

Google Scholar 
52.Rysgaard, S., Thastum, P., Dalsgaard, T., Christensen, P. B. & Sloth, N. P. Effects of salinity on NH4+ adsorption capacity, nitrification, and denitrification in Danish estuarine sediments. Estuaries 22, 21–30 (1999).CAS 
Article 

Google Scholar 
53.Peng, X. et al. Long-term fertilization alters the relative importance of nitrate reduction pathways in salt marsh sediments. J. Geophys. Res. Biogeosci. 121, 2082–2095 (2016).CAS 
Article 

Google Scholar 
54.Brown, J. H., Gillooly, J. F., Allen, A. P., Savage, V. M. & West, G. B. Toward a metabolic theory of ecology. Ecology 85, 1771–1789 (2004).Article 

Google Scholar 
55.Taylor, A. E., Giguere, A. T., Zoebelein, C. M., Myrold, D. D. & Bottomley, P. J. Modeling of soil nitrification responses to temperature reveals thermodynamic differences between ammonia-oxidizing activity of archaea and bacteria. ISME J. 11, 896–908 (2017).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
56.Ouyang, Y., Norton, J. M. & Stark, J. M. Ammonium availability and temperature control contributions of ammonia oxidizing bacteria and archaea to nitrification in an agricultural soil. Soil Biol. Biochem. 113, 161–172 (2017).CAS 
Article 

Google Scholar 
57.Mukhtar, H., Lin, Y. -P., Lin, C. -M. & Lin, Y. -R. Relative abundance of ammonia oxidizing archaea and bacteria influences soil nitrification responses to temperature. Microorganisms 7, 526 (2019).
Google Scholar 
58.Fierer, N., Carney, K. M., Horner-Devine, M. C. & Megonigal, J. P. The biogeography of ammonia-oxidizing bacterial communities in soil. Microb. Ecol. 58, 435–445 (2009).PubMed 
Article 
PubMed Central 

Google Scholar 
59.Park, H.-D., Lee, S.-Y. & Hwang, S. Redundancy analysis demonstration of the relevance of temperature to ammonia-oxidizing bacterial community compositions in a full-scale nitrifying bioreactor treating saline wastewater. J. Microbiol. Biotechnol. 19, 346–350 (2009).PubMed 
Article 
CAS 
PubMed Central 

Google Scholar 
60.Avrahami, S., Liesack, W. & Conrad, R. Effects of temperature and fertilizer on activity and community structure of soil ammonia oxidizers. Environ. Microbiol. 5, 691–705 (2003).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
61.Avrahami, S. & Conrad, R. Patterns of community change among ammonia oxidizers in meadow soils upon long-term incubation at different temperatures. Appl. Environ. Microbiol. 69, 6152–6164 (2003).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
62.Seitzinger, S. P., Gardner, W. S. & Spratt, A. K. The effect of salinity on ammonium sorption in aquatic sediments—implications for benthic nutrient recycling. Estuaries 14, 167–174 (1991).CAS 
Article 

Google Scholar 
63.Dollhopf, S. L. et al. Quantification of ammonia-oxidizing bacteria and factors controlling nitrification in salt marsh sediments. Appl. Environ. Microbiol. 71, 240–246 (2005).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
64.Beman, J. M., Bertics, V. J., Braunschweiler, T. & Wilson, J. M. Quantification of ammonia oxidation rates and the distribution of ammonia-oxidizing archaea and bacteria in marine sediment depth profiles from Catalina Island, California. Front. Microbiol. 3, 263 (2012).CAS 
PubMed 
PubMed Central 

Google Scholar 
65.Nicol, G. W., Leininger, S., Schleper, C. & Prosser, J. I. The influence of soil pH on the diversity, abundance and transcriptional activity of ammonia oxidizing archaea and bacteria. Environ. Microbiol. 10, 2966–2978 (2008).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
66.Lehtovirta, L. E., Prosser, J. I. & Nicol, G. W. Soil pH regulates the abundance and diversity of group 1.1c Crenarchaeota. FEMS Microbiol. Ecol. 70, 367–376 (2009).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
67.Bello, M. O., Thion, C., Gubry-Rangin, C. & Prosser, J. I. Differential sensitivity of ammonia oxidising archaea and bacteria to matric and osmotic potential. Soil Biol. Biochem. 129, 184–190 (2019).CAS 
Article 

Google Scholar 
68.Fuchslueger, L. et al. Effects of drought on nitrogen turnover and abundances of ammonia-oxidizers in mountain grassland. Biogeosciences. 11, 6003–6015 (2014).Article 

Google Scholar 
69.Thion, C. & Prosser, J. I. Differential response of nonadapted ammonia-oxidising archaea and bacteria to drying-rewetting stress. FEMS Microbiol. Ecol. 90, 380–389 (2014).CAS 
PubMed 
PubMed Central 

Google Scholar 
70.Fowler, S. J., Palomo, A., Dechesne, A., Mines, P. D. & Smets, B. F. Comammox Nitrospira are abundant ammonia oxidizers in diverse groundwater-fed rapid sand filter communities. Environ. Microbiol. 20, 1002–1015 (2018).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
71.How, S. W., Chua, A. S. M., Ngoh, G. C., Nittami, T. & Curtis, T. P. Enhanced nitrogen removal in an anoxic-oxic-anoxic process treating low COD/N tropical wastewater: low-dissolved oxygen nitrification and utilization of slowly-biodegradable COD for denitrification. Sci. Total Environ. 693, 133526 (2019).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
72.Gonzalez-Martinez, A., Rodriguez-Sanchez, A., van Loosdrecht, M. C. M., Gonzalez-Lopez, J. & Vahala, R. Detection of comammox bacteria in full-scale wastewater treatment bioreactors using tag-454-pyrosequencing. Environ. Sci. Pollut. Res. 23, 25501–25511 (2016).CAS 
Article 

Google Scholar 
73.Wang, D. -Q., Zhou, C. -H., Nie, M., Gu, J. -D. & Quan, Z. -X. Abundance and niche specificity of different types of complete ammonia oxidizers (comammox) in salt marshes covered by different plants. Sci. Total Environ. 768, 144933 (2021).
Google Scholar 
74.Xia, F. et al. Ubiquity and diversity of complete ammonia oxidizers (comammox). Appl. Environ. Microbiol. 84, e01390 (2018).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
75.Yu, C. et al. Evidence for complete nitrification in enrichment culture of tidal sediments and diversity analysis of clade a comammox Nitrospira in natural environments. Appl. Microbiol. Biotechnol. 102, 9363–9377 (2018).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
76.Zhao, Z. et al. Abundance and community composition of comammox bacteria in different ecosystems by a universal primer set. Sci. Total Environ. 691, 146–155 (2019).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar  More

1.Pomella, A. W. V., Barreto, R. W. & Charudattan, R. Nimbya alternantherae a potential biocontrol agent for alligatorweed, Alternanthera philoxeroides. Biocontrol 52, 271–288 (2007).CAS 
Article 

Google Scholar 
2.Barreto, R. W. & Torres, A. N. L. Nimbya alternantherae and Cercospora alternantherae: two new records of fungal pathogens on Alternanthera philoxeroides (alligatorweed) in Brazil, Australas. Plant Pathol 28, 103–107 (1999).
Google Scholar 
3.Ridenour, W. M. & Callaway, R. M. The relative importance of allelopathy in interference: the effects of an invasive weed on a native bunchgrass. Oecologia 126, 444–450 (2001).ADS 
Article 

Google Scholar 
4.Tanveer, A., Ali, H. H. & Manalil, S. Eco-biology and management of alligator weed [Alternanthera philoxeroides (Mart.) Griseb.] a review. Wetlands 38, 1067–1079 (2018).Article 

Google Scholar 
5.Garbari, F. & Pedullà, M. L. Alternanthera philoxeroides (Mart) Griseb (Amaranthaceae), a new species for the exotic flora of Italy. J. Plant Taxon Geogr 56, 139–143 (2001).
Google Scholar 
6.Chen, X., Wang, R. & Cao, Q. The relationship between the distribution of invasive plant Alternanthera philoxeroides and soil properties is scale-dependent. Pol. J. Environ. Stud. 24, 1931–1938 (2015).Article 

Google Scholar 
7.Wang, T., Hu, J. & Miao, L. The invasive stoloniferous clonal plant Alternanthera philoxeroides outperforms its co-occurring non-invasive functional counterparts in heterogeneous soil environments-invasion implications. Sci. Rep. 6, 38036 (2016).ADS 
CAS 
PubMed Central 
Article 
PubMed 

Google Scholar 
8.Wang, B., Li, W. & Wang, J. Genetic diversity of Alternanthera philoxeroides in China. Aquat. Bot. 81, 277–283 (2005).Article 

Google Scholar 
9.Shen, J., Shen, M. & Wang, X. Effect of environmental factors on shoot emergence and vegetative growth of alligatorweed (Alternanthera philoxcroides). Weed Sci. 53, 471–478 (2005).CAS 
Article 

Google Scholar 
10.Bassett, I., Paynter, Q. & Hankin, R. Characterising alligator weed (Alternanthera philoxeroides; Amaranthaceae) invasion at a northern New Zealand lake. N. Z. J. Ecol. 36, 216–222 (2012).
Google Scholar 
11.Pan, X. Y. Invasive Alternanthera philoxeroides: biology, ecology and management. Acta Phytotaxonomica Sinica 45, 884–900 (2007).Article 

Google Scholar 
12.Phung, T., Xuan, T. & Tu, A. T. Weed suppressing potential and isolation of potent plant growth inhibitors from Castanea crenata Sieb. et Zucc. Molecules 23, 345 (2018).Article 
CAS 

Google Scholar 
13.Yu, Z. & Bi, H. Status Quo of research on ecosystem services value in China and suggestions to future research. Energy Procedia 5, 1044–1048 (2011).Article 

Google Scholar 
14.Yang, S., Wang, Q. & Hu, T. Physiological responses to allelopathy of decomposing Cinnamomum septentrionale leaf litter of three crops (corn, cucumber, and cowpea). Chin. J. App. Environ. Biol. 29, 292–298 (2018).
Google Scholar 
15.Dong, B. C., Fu, T. & Luo, F. L. Herbivory-induced maternal effects on growth and defense traits in the clonal species Alternanthera philoxeroides. Sci. Total Environ. 605–606, 114–123 (2017).ADS 
Article 
CAS 

Google Scholar 
16.Dugdale, T. M., Clements, D. & Hunt, T. D. Alligatorweed produces viable stem fragments in response to herbicide treatment. J. Aquat. Plant Manag. 48, 84–91 (2010).
Google Scholar 
17.Clements, D., Dugdale, T. M. & Butler, K. L. Management of aquatic alligator weed in an early stage of invasion. Manag. Biol. Invas. 5, 327–339 (2014).Article 

Google Scholar 
18.Clements, D., Dugdale, T. M. & Butler, K. L. Herbicide efficacy for aquatic Alternanthera philoxeroides management in an early stage of invasion: integrating above-ground biomass, below-ground biomass and viable stem fragmentation. Weed Res. 57, 257–266 (2017).CAS 
Article 

Google Scholar 
19.Bond, W. & Grundy, A. Non-chemical weed management in organic farming systems. Weed Res. 41, 383–405 (2001).Article 

Google Scholar 
20.Schooler, S., Cook, T. & Bourne, A. Selective herbicides reduce alligator weed (Alternanthera Philoxeroides) biomass by enhancing competition. Weed Sci. 56, 259–264 (2008).CAS 
Article 

Google Scholar 
21.Sainty, G., Mccorkelle, G., & Julien, M. Control and spread of alligator weed Alternanthera philoxeroides (Mart.) Griseb., in Australia: Lessons for other regions. Wetlands Ecol. Manage. 5, 195–201 (1997).Article 

Google Scholar 
22.Annett, R., Habibi, H. R. & Hontela, A. Impact of glyphosate and glyphosate-based herbicides on the freshwater environment. J. Appl. Toxicol. 34, 458–479 (2014).CAS 
Article 

Google Scholar 
23.Bais, H. P., Vepachedu, R. & Gilroy, S. Allelopathy and exotic plant invasion: from molecules and genes to species interactions. Science 301, 1377–1380 (2003).ADS 
CAS 
Article 

Google Scholar 
24.Zhang, Z., Deng, L. L. & Wang, L. C. Allelopathic potential of Phragmites australis extracts on the growth of invasive plant Alternanthera philoxeroides. Allelopath. J. 45, 54–63 (2018).Article 

Google Scholar 
25.Jabran, K., Mahajan, G. & Sardana, V. Allelopathy for weed control in agricultural systems. Crop Prot. 72, 57–65 (2015).Article 

Google Scholar 
26.Kumbhar, B. A. & Patel, D. D. Allelopathic effects of different weed species on crop. J. Pharm. Sci. Biosci. Res. 6, 801–805 (2016).
Google Scholar 
27.Weston, L. A. & Duke, S. O. Weed and crop allelopathy. Crit. Rev. Plant Sci. 22, 367–389 (2003).CAS 
Article 

Google Scholar 
28.Rice, E.L., Allelopathy, 2nd Ed, A. Press, Editor. 1-50 (1984).29.Da Silva, I. F. & Vieira, E. A. Phytotoxic potential of Senna occidentalis (L.) Link extracts on seed germination and oxidative stress of Ipe seedlings. Plant Biol. 21, 770–779 (2019).Article 
CAS 

Google Scholar 
30.Zeng, R. S. Allelopathy—the solution is indirect. J. Chem. Ecol. 40, 515–516 (2014).CAS 
Article 

Google Scholar 
31.Otusanya, O. O., Ilori, O. J. & Adelusi, A. A. Allelopathic effects of tithonia diversifolia (Hemsl) A. gray on germination and growth of Amaranthus cruentus. Res. J. Environ. Sci. 1, 285–293 (2007).CAS 
Article 

Google Scholar 
32.Chen, Z. & Meng, S. Research progress of Humulus scandens. Chin. Pharm. Affairs 24, 73–77 (2011).CAS 

Google Scholar 
33.Cao, Y., Wang, T. & Xiao, Y. A. The interspecific competition between Humulus scandens and Alternanthera philoxeroides. J. Plant Interactions 9, 194–199 (2013).Article 
CAS 

Google Scholar 
34.Huang, Y. M., Zhang, Y. & Liu, Q. Research on allelopathy of aqueous extract from tagetes patula to four garden plants. Acta pratacultural sinica (Chinese) 24, 150–158 (2015).
Google Scholar 
35.Li, W., Luo, J. & Tian, X. A new strategy for controlling invasive weeds: selecting valuable native plants to defeat them. Sci. Rep. 5, 11004 (2015).ADS 
CAS 
PubMed Central 
Article 
PubMed 

Google Scholar 
36.Cheng, T. S. The toxic effects of diethyl phthalate on the activity of glutamine synthetase in greater duckweed (Spirodela polyrhiza L.). Aquat. Toxicol. 124, 171–178 (2012).Article 
CAS 

Google Scholar 
37.Makoi, J. H. & Ndakidemi, P. A. Biological, ecological and agronomic significance of plant phenolic compounds in rhizosphere of the symbiotic legumes. Afr. J. Biotech. 6, 739–748 (2007).
Google Scholar 
38.Zhang, K. M., Shen, Y. & Yang, J. The defense system for Bidens pilosa root exudate treatments in Pteris multifida gametophyte. Ecotoxicol. Environ. Saf. 173, 203–213 (2019).CAS 
Article 

Google Scholar 
39.Noctor, G., Reichheld, J.-P. & Foyer, C. H. ROS-related redox regulation and signaling in plants. Semin. Cell Dev. Biol. 80, 3–12 (2017).Article 
CAS 

Google Scholar 
40.Kaur, N., Chugh, V. & Gupta, A. K. Essential fatty acids as functional components of foods-a review. J. Food Sci. Technol. 51, 2289–2303 (2014).CAS 
Article 

Google Scholar 
41.Sun, C. H., Li, Y. & He, H. Y. Physiological and biochemical responses of Chenopodium album to drought stresses. Atca Ecologica Sinica 25, 2556–2561 (2005).CAS 

Google Scholar 
42.Li, P., Wang, X. & Li, Y. The contents of phenolic acids in continuous cropping peanut and their allelopathy. Acta Ecol. Sin. 30, 2128–2134 (2010).MathSciNet 
Article 

Google Scholar 
43.Wang, Y. X., Sun, G. R. & Wang, J. B. Relationships among MDA content, plasma membrane permeability and the chlorophyll fluorescence parameters of Puccinellia tenuiflora seedlings under NaCl stress. Acta Ecol. Sin. 26, 122–129 (2006).CAS 

Google Scholar 
44.Li, Z. Q., Li, J. T. & Bing, J. The role analysis of APX gene family in the growth and developmental processes and in response to abiotic stresses in Arabidopsis thaliana. Hereditas (Beijing) 41, 534–547 (2019).
Google Scholar 
45.Tang, K., Ming, L. & Shan, D. Allelopathy autotoxicity effects of aquatic extracts from rhizospheric soil on rooting and growth of stem cuttings in Pogostemon cablin. J. Chin. Med. Mater. 37, 935–939 (2014).CAS 

Google Scholar 
46.Coelho, E. M. P., Barbosa, M. C. & Mito, M. S. The activity of the antioxidant defense system of the weed species Senna obtusifolia L and its resistance to allelochemical stress. J. Chem. Ecol. 43, 725–738 (2017).CAS 
Article 

Google Scholar 
47.Li, J. & Wang, X. Advance of research on Humulus scandens. Qilu Pharm. Affairs 26, 353–355 (2007).
Google Scholar 
48.Xu, B., Jin, Y. & Yihan, W. Chemical constituents from stems and leaves of Humulus scandens. Chin. Tradit. Herb. Drugs 45, 1228–1231 (2014).CAS 

Google Scholar 
49.Zhang, J., Liu, J. & Dai, L.-F. Unlocking the potential antioxidant and anti-inflammatory activities of Rhododendron molle G. Don. Pak. J. Pharm. Sci. 32, 2375–2383 (2019).CAS 

Google Scholar 
50.Chen, Z., Guo, Q. & Huang, K. Analysis of volatile components of solidago canadensis by SPME/GC-MS. Acta Agriculturae Jiangxi (Chinese) 26, 1 (2008).
Google Scholar 
51.Wang, Y., Yu, J. & Zhang, Y. Effects of two allelochemicals on growth and physiological characteristics of eggplant seedlings. J. Gansu Agric. Univ. (Chinese) 3, 47–50 (2007).
Google Scholar 
52.Yu, J., Zhang, Y. & Niu, C. Effects of two kinds of allelochemicals on photosynthesis and chlorophyll fluorescence parameters of Solanum melongena L. seedlings. J. Appl. Ecol. 17, 1629–1632 (2006).ADS 
CAS 

Google Scholar 
53.Lande, M. L., Kanedi, M. & Zulkifli, Z. Supperssive effexts of lantana camara leaf extracts on the growth of red chillli (Carsicum annuum). World J. Pharm. Life Sci. 3, 543–551 (2017).
Google Scholar 
54.Erida, G. & Saidi, N. Allelopathic screening of several weed species as potential bioherbicides. IOP Conf. Ser. Earth Environ. Sci. 334, 12–34 (2019).Article 

Google Scholar 
55.Cimmino, A., Masi, M. & Rubiales, D. Allelopathy for parasitic plant management. Nat. Prod. Commun. 13, 289–294 (2018).
Google Scholar 
56.Deng, J., Zhang, Y. & Hu, J. Autotoxicity of phthalate esters in tobacco root exudates: Effects on seed germination and seedling growth. Pedosphere 27, 1073–1082 (2017).CAS 
Article 

Google Scholar 
57.Gu, S., Zheng, H. & Xu, Q. Comparative toxicity of the plasticizer dibutyl phthalate to two freshwater algae. Aquat. Toxicol. 191, 122–130 (2017).CAS 
Article 

Google Scholar 
58.Perveen, S., Yousaf, M. & Zahoor, A. F. Extraction, isolation, and identification of various environment friendly components from cock’s comb (Celosia argentea) leaves for allelopathic potential. Toxicol. Environ. Chem. Rev. 96, 1523–1534 (2014).CAS 
Article 

Google Scholar 
59.Alara, O. R., Abdurahman, N. H. & Ukaegbu, C. I. Extraction and characterization of bioactive compounds in Vernonia amygdalina leaf ethanolic extract comparing soxhlet and microwave-assisted extraction techniques. J. Taibah Univ. Sci. 13, 414–422 (2019).Article 

Google Scholar 
60.Wei, W., Hou, Y. & Peng, S. Effects of light intensity on growth and biomass allocation of invasive plants Mikania micrantha and Chromolaena odorata. Acta Ecol. Sin. 37, 6021–6028 (2017).Article 

Google Scholar 
61.Williamson, G. B. & Richardson, D. Bioassays for allelopathy: measuring treatment responses with independent controls. J. Chem. Ecol. 14, 181–187 (1988).Article 

Google Scholar 
62.Gao, Y. B., Li, G. P. & Shi, H. Allelopathic effect of endophyte-infected achnatherum sibiricum on stipa grandis. Acta Ecol. Sin. 37, 1063–1073 (2017).Article 

Google Scholar 
63.Zhang, L., Wang, X. & Guo, J. Metabolic profiling of Chinese tobacco leaf of different geographical origins by GC-MS. Agric. Food Chem. 61, 2597–2605 (2013).CAS 
Article 

Google Scholar  More

Two species unidirectional coupling ecosystemThis bio-inspired ecosystem S((beta _{xy}=0), (beta _{yx})) describing the unidirectional coupling is run for 1000 time steps for reaching stationarity, generating a set of 1000 points long time-series dependent on (beta _{yx}). This means that species X has an increasing effect on Y with the increase of (beta _{yx}), but Y has no effect on X. Both CCM and the proposed OIF model are separately used to quantify the potential causality between species X and Y. The inferred causality dependent on (beta _{yx}) only (as a physical interaction) is shown in Fig. 2A. (beta), (rho) and TE have different units; specifically (beta) and (rho) are dimensionless while TE is measured in bits or nats (a logarithmic unit of information or entropy). Therefore, any comparison is done considering gradients of change when these variables vary together rather than making comparisons between absolute values which are meaningless. Figure 2A shows that under the condition of (beta _{xy})=0, results of “Y to X” (i.e. the estimated effect on Y on X) is close to 0 for the OIF model ((TE_{Y rightarrow X}(beta _{yx}))) that precisely describe the no-effect of Y on X. “X to Y” ((TE_{X rightarrow Y}(beta _{yx}))) well tracks the increasing strength of the effect of X on Y for increasing values of the physical interaction (beta _{yx}) embedded into the mathematical model. However, considering results of the CCM model, “Y to X” ((rho _{Y rightarrow X}(beta _{yx}))) presents non obvious (and likely wrong) non-zero values with higher fluctuations compared to (TE_{Y rightarrow X}(beta _{yx})) especially for lower values of (beta _{yx}). This erroneous estimates of CCM is likely related to the need of CCM for convergence. For CCM, “X to Y” (((rho _{X rightarrow Y}(beta _{yx})))) shows an increasing trend for increasing values of (beta _{yx}) and decreasing when (beta _{yx}) is greater than (sim)0.5 non-trivially. In consideration of these results for the unidirectional coupling ecosystem, the OIF model performs better over CCM in terms of unidirectional causality inference.Two species bidirectional coupling ecosystemIn this case, the effect between two species is bidirectional. Species X has an effect on species Y and vice versa. The univariate dynamical systems S(0.2/0.5/0.8, (beta _{yx})) are run for 1000 time steps under the same conditions determined by (beta _{xy}). Certainly this situation is fictional since in real ecosystems the interaction strength is changing when other interacting species change their interactions.Thus, keeping one interaction fixed around one value is a strong unrealistic simplification (analogous of one-factor at-a-time sensitivity analyses) but it is a toy model that allows to verify the power of network inference models. These models generate three sets of 1000 points long time-series dependent of (beta _{yx}) for each fixed (beta _{xy}). OIF and CCM are used to infer “causality” between X and Y—in the form of (rho) and TE—and compare that against the real embedded interaction (beta _{yx}) and (beta _{xy}) shown in Fig. 2B,C,D. Considering all results of Fig. 2 corresponding to fixed (beta _{xy})s, the correlation coefficient (rho) yielded from CCM and TE from OIF are both able to track the strength of causal trajectories. However, TE seems to perform better in term of ability to infer fine-scale changes in interactions. In particular, considering Fig. 2D (right plot), higher (TE_{yx}) higher for low (beta _{yx}) makes sense because (beta _{xy} > beta _{yx}) that means Y has a larger influence on X than vice versa and then Y is able to predict X. Additionally, TE does not suffer of convergence problems; specifically, considering Fig. 2A (left plot), higher (rho) for small (beta _{yx}) is not sensical and that is likely related to convergence problems of CCM.Considering all results of Fig. 2 corresponding to fixed (beta _{xy})s, the correlation coefficient (rho) yielded from CCM and TE from OIF are both able to track the strength of causal trajectories. Ideally, the causality from Y to X is a constant since (beta _{xy}) is a fixed value for each case. In this figure, the red curve in the right panel representing the OIF-inferred (TE-based) causality from Y to X is higher for greater (beta _{xy})s, while red curves representing CCM-inferred ((rho)) causality in the left panel present higher fluctuations especially for lower (beta _{yx}). For the causality from X to Y determined by (beta _{yx}) in the mathematical model, theoretically speaking, the causality from X to Y should monotonously grow when (beta _{yx}) increases from 0 to 1. In Fig. 2, blue curves in the right panel representing the OIF-inferred (TE-based) causality from X to Y present monotonously increasing features as a whole with the increasing (beta _{yx}), while those from CCM model ((rho)) do not and show considerable fluctuations.Therefore, OIF outperforms CCM in terms of the ability to infer the fine-scale changes in causality. In particular, considering Fig. 2D (right plot), higher (TE_{yx}) higher for low (beta _{yx}) makes sense because (beta _{xy} > beta _{yx}) that means Y has a larger influence on X than vice versa and then Y is able to predict X. Additionally, TE does not suffer of convergence problems; specifically, considering Fig. 2A (left plot), higher (rho) for small (beta _{yx}) is not sensical and that is likely related to convergence problems of CCM.Additionally, (rho _{Y rightarrow X}(beta _{yx})) shows higher fluctuations on average especially for the condition of lower (beta _{yx})s compared to (TE_{Y rightarrow X}(beta _{yx})). When considering the effect of X on Y that is a function of (beta _{yx}) for CCM, (rho _{X rightarrow Y}) reaches an extreme value at around (beta _{yx}) = 0.5 and then declines for larger values of (beta _{yx}). This is not consistent with the expected effect of X on Y that should be proportional to (beta _{yx}) embedded into the mathematical model. The ability of (rho) to reflect the proportional relationship between the effect of X on Y (manifested by (beta _{yx})) vanishes for high (beta _{xy})s due to unexpected and somewhat inconspicuous changes in (rho _{X rightarrow Y}) for larger (beta _{yx}). In simple words, the expected increasing trend of (rho) is lost for larger (beta _{xy}) that is counterintuitive. On the other side, (TE_{X rightarrow Y}(beta _{yx})) invariably maintains an increasing trend for increasing values of (beta _{xy}). OIF is also performing better than CCM when predicting higher average values of (TE_{Y rightarrow X}) for increasing values of (beta _{xy}) (red curves in Fig. 2A–D, right plots) as expected by the fixed effect in the mathematical model of Y on X. These results suggest that when compared to (rho) of CCM, TE can track well the causal interactions over (beta _{yx}) with higher performance and without considering the convergence requirement of CCM. CCM needs to consider the length of time series that makes (rho _{X rightarrow Y}(beta _{yx})) convergent to a stable value, but uncertain for large differences in time-series length of (X,Y) and sensitive to short time series.In more realistic settings for real ecosystems (and in analogy to global sensitivity analyses) when (beta _{xy}) and (beta _{yx}) are both considered as arguments of the two-variable (X,Y) bio-inspired model, the simulated ecosystem becomes a truly bivariate system, yet yielding complexity but more interest into the causality inference (Fig. 3). The dynamical system S((beta _{xy}), (beta _{yx})) was generated for 800 time steps under the same conditions mentioned above. We generated the datasets that allowed us to study linear and non-linear predictability indicators for inferring the embedded physical interactions. Specifically, we measure undirected linear correlation coefficient (corr_{X;Y}(beta _{xy},beta _{yx})), non-linear undirected mutual information (MI_{X;Y}(beta _{xy},beta _{yx})), directed non-linear correlation coefficient (rho _{X rightarrow Y}(beta _{xy},beta _{yx})) and (rho _{Y rightarrow X}(beta _{xy},beta _{yx})), and non-linear directed transfer entropy (TE_{X rightarrow Y}(beta _{xy},beta _{yx})) and (TE_{Y rightarrow X}(beta _{xy},beta _{yx})) as shown in Fig. 3. These 2D phase-space maps in Fig. 3 show strikingly similar patterns for classical linear correlation coefficients, MI, (rho) of CCM and TE of OIF which underline the fact that all methods are able to infer the interdependence patterns of interacting variables explicitly defined by (beta _{xy}) and (beta _{yx}). The color of phase-space maps is proportional to the inferred interaction between X and Y when the mutual physical interactions are varying according to the mathematical model in Eq. (1). In Fig. 3, even though phase-space maps of undirected (corr_{X;Y}(beta _{xy},beta _{yx})) and (MI_{X;Y}(beta _{xy},beta _{yx})) present similar patterns (in value organization and not value range) to those of directed (rho) and TE, neither (corr_{X;Y}((beta _{xy},beta _{yx})) and (MI_{X;Y}((beta _{xy},beta _{yx})) provide information about the direction of causality. As expected MI shows the opposite pattern of the average TE due to the fact that MI is the amount of shared information (or similarity) versus the amount of divergent information (divergence and asynchronicity) between X and Y.Figure 3Phase-space maps of normalized coupling predictive causation via correlation, mutual information, CCM and OIF for varying true causal interactions. Both true causal interactions (beta _{xy}) and (beta _{yx}) are free varying within the range [0, 1], indicating a bivariate model S((beta _{xy}),(beta _{yx})) where both species (or variables more generally) are interacting with each other with different strength. (A) normalized correlation coefficient, (B) normalized mutual information, (C) and (E) normalized CCM correlation coefficient ((rho)) for interaction directions of (X rightarrow Y) and (Y rightarrow X), (D) and (F) normalized transfer entropy (TE) from OIF model for interaction directions of (X rightarrow Y) and (Y rightarrow X).Full size imageFigure 4Dynamics of abundance and predictability for the bidirectional two species ecosystem model. (A) plots refer to the species abundance in time for the mathematical model in Eq. 1 for different predictability regimes associated to different interaction dynamics from low to high complexity ecosystem associated to low and high predictability. Blue, green and red refer to a range of predictable interactions as in Fig. 3: specifically, Blue is for ((beta _yx), (beta _xy))=(0.18, 0.39) (small mutual interaction, and predominant effect of Y on X), Green is for (0.64, 0.57) (high mutual interactions, and slightly predominant effect of X on Y), and Red for (0.94, 0.34) (high mutual interactions, and predominant effect of X on Y). (B) phase-space plots showing the non-time delayed associations between X and Y corresponding to synchronous and homogeneous, mildly asynchronous and divergent, and asynchronous and divergent dynamics. The transition from synchronous/small interactions to asynchronous/high interaction leads to a transition from modular to nested ecosystem interactions when more than one species exist (Fig. 6).Full size imageIn a biological sense TE should be interpreted as the probability of likely uncooperative dynamics (leading to or driven by environmental or biological heterogeneity) while MI as the probability of cooperative dynamics (leading to or driven by homogeneity). Here we refer to cooperative and uncooperative interactions based on the similarity or dissimilarity in pair dynamics manifested by species abundance fluctuations. For instance divergence and asynchronicity (that define TE) in pair species dynamics manifest uncooperative interactions. The balance of cooperative and uncooperative interactions can result into net interactions at the ecosystem scale manifesting neutral patterns, or net interactions may lead to niche patterns biased toward strong environmental or biological factors52. Certainly, cooperation in a biological sense should be interpreted on a case by case basis. In a broader uncertainty propagation perspective49, “cooperation” between variables means that variables contribute similarly to the uncertainty propagation, while “competition” means that one variable is predominant over the other in terms of magnitude of effects since TE is proportional to the magnitude rather than the frequency of effects. For the former case the total entropy of the system is higher than the latter case. Interestingly, correlation (text{ corr }(beta )), (rho) and TE show similar patterns in both organization and value range (but not in singular values of course), which sheds some important conclusions about the similarity and divergence of these methods as well as their capacity and limitations in characterizing non-linear systems.When comparing the phase-space patterns from CCM and OIF (displaying (rho) and TE) a more colorful and informative pattern is revealed by OIF. This means that TE gives a better gradient when tracking the increasing strength of causality for increasing values of (beta _{xy}) and (beta _{yx}). When comparing the phase-space patterns for the two causal directions of (“X rightarrow Y”) and (“Y rightarrow X”), phase-space maps from CCM are very similar, while those from TE present apparent differences in the strength of effects for the two opposite direction of interaction. Therefore, OIF is more sensitive to the direction of interaction compared to CCM when detecting directional causality.These results imply that TE performs better to distinguish directional embedded physical interactions (that are dependent on direct interactions (beta)-s, species growth rate (r_x) and (r_y), and contingent values X(t) and Y(t) determining the total interaction as seen in the model of Eq. (1)) in the species causal relationships. It should be emphasized how all linear and non-linear interaction indicators are inferring the total interaction and not only those exerted by (beta)-s. In a broad uncertainty purview49 the importance of these three factors ((beta)-s, r-s and X(t)/Y(t)) depends on their values and probability distributions that define the dynamics of the system; dynamics such as defined by the regions identified by patterns in Fig. 3 for the predator-prey system in Eq. (1). In principle, the higher the difference between these three interaction factors in the species considered, the higher the predictability and sensitivity of OIF. Figure 4 highlights three different dynamics corresponding to the TE blue, green and red regions in Fig. 3.In all dynamical states represented by Fig. 3, species are interacting with different magnitudes and this defines distinct network topologies. Three prototypical dynamics are shown in Fig. 4 with colors representative of (rho) and TE in Fig. 3. The “blue” deterministic dynamics has very high synchronicity and no divergence considering variable fluctuation range (the gap is deterministic and related to the numerically imposed (u=1)), as well as no linear correlation between non-lagged variables. In perfect synchrony one would have one point in the phase-space. Thus, absence of correlation does not imply complete decoupling of species but it can be a sign of small interactions. The “green” dynamics shows a relatively high synchronicity and medium divergence. In the phase-space of synchronous values of X and Y a correlation is observed with relatively small fluctuations because the divergence is small. Lastly, the “red” dynamics shows a relatively high asynchronicity and divergence. The stochasticity is higher than previous dynamics and the “mirage correlation” in the phase space has higher variance. Time-dependent mirage correlations in sign and magnitude mean that correlation (that may suggest common dynamics in a linear framework) does not imply similarity in dynamics for the two species. Non-linearity is higher from blue to red dynamics as well as predictability but lower absolute information entropy. Then, it is safe to say that linear dynamics (or small stochasticity) does not imply higher predictability.Real-world sardine–anchovy-temperature ecosystemCCM and proposed OIF model are also used for a real-world fishery ecosystem to infer potential causal interactions between Pacific sardines (Sardinops sagax) landings, Northern anchovies (Engraulis mordax) and sea-surface temperature (SST) recorded at Scripps Pier and Newport Pier, California. Sardines and anchovies do not interact physically (or the interaction is low in number), while both of them are influenced by the external environmental SST that is the external forcing. To quantify the likely causal interactions between species and SST based on real data, we use CCM considering the length of time series for convergence of (rho), as well as OIF considering a set of time delays for acquiring stable values of inferred interactions TEs.Figure 5Inferred predictive causality for the sardine-anchovy-Sea Surface Temperature ecosystem. CCM correlation coefficient ((rho)) and OIF predictor (TE) are shown in the left and middle plots for different pairs considered (sardine–anchovy, sardine and SST, anchovy and SST from top to bottom).Full size imageResults from CCM in Fig. 5A (plots from top to bottom) show that no significant interaction can be claimed between sardines and anchovies, as well as from sardines or anchovies in the SST manifold which expectedly indicates that neither sardines nor anchovies affect SST. This latter results, considering its biological plausibility should be taken as one validation criteria of predictive models, or complimentary as a test for anomaly detection of spurious interactions. The reverse effect of SST on sardines and anchovies can be quantitatively detected with the correlation coefficient (rho) as well as TE. Although the calculated causations between SST and sardines or anchovies are moderate, CCM is able to provide a good performance in causality inference when the length of time series used is long enough due to convergence requirement.Figure 5B shows OIF’s results of inferred causal interactions between sardines, anchovies and SST dependent on the time delay u. For sardines and anchovies, OIF exposes bidirectional interactions that are actually biologically plausible, especially when both populations coexist in the same habitat, versus the results of CCM that infer (rho =0). Ecologically speaking, even though fish populations do not directly influence sea temperature, we can find some clues about SST in fish populations influenced by SST. These clues can be interpreted as information of SST encoded in fish populations over abundance time records. So, observations of fish populations can be used to inversely predict the change of SST; this can be interpreted as “reverse predictability” (or “biological hindcasting”) in a similar way of when predicting historical climate change from ice cores. This information is captured by OIF, leading to nonzero values of TE from fish populations to SST. In this regard, we emphasize the distinction between direct and indirect (reverse) information flow, where direct information flow is most of the time larger and signifies causality (e.g. of SST for sardine and anchovies), and indirect (reverse) information flow that is typically smaller and signifies predictability (e.g. sardine and anchovies for ocean fluctuations). It is possible—especially for linear systems where an effect is observed immediately after a change—that information of SST encoded in fish populations is high if the interdependence, represented by the functional time delay u, of the environment-biota is small. However, for highly non-linear systems such as fishes and the ocean, changes in temperature may take a while before being encoded into fish population abundance53. Thus, it is correct that the highest values of TE are for high u. Values of TE for small u-s are numerical artifacts related to systematic errors leading to overestimation of interactions that are time-delayed eventually. One way to circumvent this problem, largely present for short time series, would be to extend time series by conserving their dynamics (see Li and Convertino39) or to bound the calculation of TE only for the u that maximizes the Mutual Information; this would provide an average u within a range where TE is approximately invariant. Thus, for the effect of external SST on sardines and anchovies, OIF model gives unstable causal interactions with bias for lower time delays due to known dependencies of TE on u (such as cross-correlation for instance) that establishes the temporal lag on which the dependency between X and Y is evaluated. In a sense, plots in Fig. 5B are like cross-variograms for the pairs of variables considered. TE becomes stable when the time delay is located in an appropriate range. It means that OIF requires an optimal time delay that makes results of the causality inference robust and that is related to optimal TEs (as highlighted in Li and Convertino39 and Servadio and Convertino49) that defines the most likely interdependency between variables for the u with the highest predictability. The fact that TE of sardine and anchovies to SST is high for same small ranges of u may be also a byproduct of data sampling, i.e., fish and SST sampling locations are different (fish abundance is actually about fish landings) and that can introduce spurious correlations/causation. Overall, these findings suggest that the OIF model provides more plausible results, but it requires careful selection of optimal time delays.Figure S1 shows the relationships between normalized (rho) and TE estimated for all selected values of L and u of pairs in Fig. 5 (sardine-anchovy, sardine and SST, anchovy and SST). These plots show opposite results than the proportionality between (rho) and TE in Fig. 3 because non-optimal values are used, that is non-convergent (rho)-s and suboptimal TE during the interaction inference procedure (Fig. S1). TE for too small u-s determines overestimation of interactions due to the implicit assumptions that variables have an immediate effect on each other and that is not always the case as highlighted by the vast time-lagged determined non-linear regions in Fig. 3. If “transitory” values of (rho) for small L are disregarded, as well as TEs for small u-s, the relationship between (rho) and TE shows a correct linear proportionality.Real-world multispecies ecosystemInteractions between fish species living in the Maizuru bay are intimately related to external environmental factors of the ecosystem where they live, the number of species living in this region considering also the unreported ones and biological species interactions, which lead to a complex dynamical nonlinear system. In Fig. 6 the network of observed fish species (Table S1) is reported where only the interactions considered in Ushio et al.37 for the CCM are reported. This is because the goal is to compare the CCM inferred network to the TE-based one based on abundance. Figure 7 shows the temporal fluctuations of abundance and the functional interaction matrices of (rho) and TE. In this paper we study and compare average ecosystem networks for the whole time period considered but dynamical networks can also be extracted via time-fluctuating (rho) and TE as shown in Fig. S3. These dynamical networks can be useful for studying how diversity is changing over time and ecosystem stability (Figs. S4, S6–S7) as well as understanding the relationship between (rho) and TE (Fig. S5). In the network of Fig. 6 the color and width of links are proportional to the magnitude of TE (Table S2); for the former a red-blue scale is adopted where the red/blue is for the highest/lowest TEs. The diameter is proportional to the Shannon entropy of the species abundance pdf (Table S3). The color of nodes is proportional to the structural node degree, i.e. how many species are interconnected to others. Therefore, the network in Fig. 6 is focusing on uncooperative species whose divergence and/or asynchronicity (that is a predominant factor in determining TE over divergence) is large. Yet, the connected species are rarely but strongly interacting in magnitude rather than frequently and weakly (i.e., cooperative or similar dynamics). Additionally, the species with the smallest variance in abundance are characterized by the smallest Shannon entropy (smallest nodes) and more power-law distribution although the latter is not a stringent requirement since both pdf shape and abundance range (in particular maximum abundance) play a role in the magnitude of entropy. Average entropy such as average abundance are quantities with limited utility in understanding the dynamics of an ecosystem as well as ecological function. Nonetheless, species with high average abundance (e.g. species 5) have a very regular seasonal oscillations and the largest number of interactions with divergent species. This dynamics is expected considering the population size of these dominant species and their synchrony with regular environmental fluctuations.Figure 6Part of the estimated species interaction network for the Maizuru Bay ecosystem. Species properties are reported in Table S1. The color and width of links are proportional to the magnitude of TE (Table S2); for the former a red-blue scale is adopted where the red/blue is for the highest/lowest TEs. The diameter is proportional to the Shannon entropy of the species abundance (Table S3) that is directly proportional to the degree of uniformity of the abundance pdf and the diversity of abundance values (e.g., the higher the zero abundance instances the lower the entropy). The color of nodes is proportional to the structural node degree, i.e. how many species are interconnected to others after considering only the CCM derived largest interactions (see Ushio et al.37 and Fig. 7). Other interactions exist between species as reported in Fig. 7. TE is on average proportional to (rho) (Figs. S4 and S5). Freely available fish images are from FishBase https://www.fishbase.in/search.php; the network was created in Matlab and the composition of network and images was made in Adobe Illustrator version 21 (2017) https://www.adobe.com/products/illustrator.html.Full size imageFigure 7Normalized species interactions matrices inferred by CCM and OIF models for Maizuru Bay ecosystem. In the census of the aquatic community, 15 fish species were counted in total. Interaction inferential models use time lagged abundance magnitude (CCM) or pdfs of abundance (OIF) shown in (A). (B) normalized CCM correlation coefficients ((rho)) between all possible pairs of species. (C) normalized transfer entropies (TEs) between all pairs of species from the OIF model. Both CCM and OIF predict that the most interacting species (in terms of magnitude rather than frequency) are 7, 8 and 9 on average. Thus, interaction matrices are more proportional to the asynchronicity than the divergence of species in terms of abundance pdf, although abundance value range defines the uncertainty (and diversity) for each species that ultimately affects entropy and interactions (e.g., if one species have many zero abundance instances or many equivalent values, such as species 2, TEs of that species are expected to be low due to lower uncertainty despite the asynchrony and divergence).Full size imageFigure S2 shows that the strongest linear correlation is for the most divergent and asynchronous species (from species 4–9) for which both (rho) and TE are the highest (Fig. 7B,C). This confirms the results of Fig. 3 and the fact that competition (or dynamical diversity more generally) increases predictability. This also highlights the fact that linear correlation among state variables does not imply synchronicity or dynamic similarity as commonly assumed. The interaction matrices in Fig. 7B, C confirm that TE has the ability to infer a larger gradient of interactions than (rho) and the total entropy of the TE matrix is lower than (rho). Pairwise the inferred interaction values by CCM and OIF are different but (rho) and TE patterns appear clearly similar and yet proportional to each other.CCM and OIF models are applied to calculate the potential interactions between all pairs of species. Figure 7B,C show interaction matrices describing the normalized (rho) from CCM and TE from OIF model of all pairwise species, respectively. The greater the strength of likely interaction, the warmer the color. These results demonstrate that CCM and OIF model present similar patterns for the interaction matrices in terms of interaction distribution, gradient and magnitude in order of similarity. This indicates that both CCM and OIF are able to infer the potentially causal relationships between species. Compared to the CCM interaction heatmap the OIF heatmap presents larger gradients of inferred interactions that highlight the divergence and asynchrony in fish populations of species 4–9 from other species. This difference can be observed in Figure S2 that shows the strongest linear correlations for the most divergent and asynchronous species (4–9). It is worth noting that species 4-8 are all native species (See Table SS11). Therefore, despite the patterns of interactions of CCM and TE are similar, CCM allows one a better identification of clusters of species with similar or distinct interaction ranges. Additionally TE estimates some weak observed interactions such as of species 2 (E. japonicus) with others, while CCM essentially considers null interactions for these species.Precisely, the most interacting species (4–9) are the most divergent and asynchronous species (with respect to the whole community) as well as diverse in terms of values of abundance; these species form the ”collective core” that is likely determining the stability of the ecosystem. Interestingly the number of these species is relatively small and it confirms results of other studies (see39, 54,55,56,57) showing that the number of species with weak interactions is much larger. Theory suggests that this pattern promotes stability as weak interactors dampen the destabilizing potential of strong interactors54. Mediated cooperation (e.g. by many “weakly” interacting competitors) as shown by Tu et al.52 promotes biodiversity and diversity increases stability. When considering abundance values (at same time steps) of collective core species (Fig. S2) these species are linearly related and this increases their mutual predictability by either using linear or non-linear models based on correlation coefficient and TE. This proves that non-linearity increases predictability.The choice of the optimal u that maximizes MI leads to the optimal TE model and resultant interaction network. The observed u over time is really small (Fig. S9) and this signifies how likely the ecosystem has small memory and responds quickly to rapid changes, or the information of change is carried over time by ecosystem’s interactions which lead to accurate short-term forecasting. In other words, temperature-induced changes may take long time but the information of change is replicated at short time periods. The chosen time delay (u =1) corresponds to the species sampling of two weeks. Note that values of u are also dependent on the data resolution and they are strongly related to fluctuations rather than absolute (alpha)-diversity value. Thus, while biodiversity may fluctuate rapidly in time, value of (alpha)-diversity for seasons or longer time periods can be more stable and manifesting higher memory (representative of u for the whole ecosystem) than the one between species pairs (related to pair’s u). Short-term catastrophic dynamics (for instance related to dramatic habitat change, sudden invasions, extinctions or rapid adaptations) may lead to irreversible shifts in interactions (strength and sign); this, in turn can affect biodiversity patterns that are completely uninformed by past dynamics. Thus, there is certainly a limit to predictability and to the validity of time delays which can change very rapidly. However, we insist in emphasizing that models are predictive tools and predictions are not necessarily causality reflecting the many and highly complex underlying processes. Yet, interpretation of results must be done with care.We also study temporally dynamical networks for the fish ecosystem community (see “Real-world sardine–anchovy-temperature ecosystem” section). CCM and OIF model are applied to quantify the causality between all possible pairs of species at each time period by calculating (rho) and TE, respectively. Estimated effective (alpha)-diversity (Eq. S1.2) from CCM- and TE-based inferred networks at each time point can be obtained and then compared to the taxonomic (or ”real”) (alpha)-diversity. Results are shown in Fig. 8 and Fig. S6. In the whole time period, the estimated (alpha)-diversity from CCM is constant, whereas the global trend of the estimated (alpha)-diversity from OIF model slightly decreases over time that is consistent with the global trend of real (alpha)-diversity. CCM always predicts a non-zero interaction for all species (including negative values) whereas OIF predicts zero interactions for some species that are then not making part of the estimated effective (alpha)-diversity.Figure 8Predicted (alpha) -diversity via optimal interaction threshold for CCM’s (rho) and OIF’s TE versus taxonomic diversity. Effective (alpha)-diversity from CCM and OIF are shown (blue and red) for an optimal threshold of (rho) and TE (i.e., 0.2 and 0.3) that maximizes the correlation coefficient and mutual information (MI) between (alpha _{CCM}) or (alpha _{TE}) and the taxonomic (alpha), respectively. The maximization of the correlation coefficient and MI guarantees that the estimated effective (alpha) are the closest to the taxonomic (alpha). g is the resolution of the network inference determined by the minimum number of points required to construct pdfs and infer TE robustly (see Supplementary Information section S1.3).Full size imageFigure 8 shows the effective (alpha) diversity from CCM and OIF for an optimal threshold of (rho) and TE (i.e., 0.2 and 0.3) that maximizes the correlation coefficient and Mutual Information (MI) between (alpha _{CCM}) or (alpha _{TE}) and the taxonomic (alpha), respectively. The maximization of the correlation coefficient and MI guarantees that the estimated effective (alpha) are the closest to the taxonomic (alpha). Figure S6 shows effective (alpha) for unthresholded interactions and other thresholds. Note that the threshold on TE does not coincide with the value of TE that maximizes the total network entropy (Fig. 8) and then some of the reported species may not be part of the ecosystem strongly. Thus, this threshold method is also useful to identify species that are truly forming local diversity versus transient species. Considering the pattern of fluctuations of effective (alpha)-diversity from CCM, they are poorly unrelated to the real (alpha)-diversity, while those from OIF are much more synchronous with seasonal fluctuations of real (alpha)-diversity. However, (alpha _{TE}) is a bit higher than the average taxonomic (alpha). Both CCM and TE predict a decrease in (alpha) in time that corresponds to an increase in SST. As shown in Fig S6, OIF is attributing higher sensitivity to SST for small interaction species because (alpha) fluctuations show seasonality that happens when species follow environmental dynamics closely. Vice versa, CCM predicts a broader sensitivity for all positively interacting species. These results reveal that OIF gives an effective tool to measure meaningful interdependence relationships between species for constructing temporally dynamical networks where the number of nodes over time [estimated (alpha (t))] can reflect closely the taxonomic (alpha)-diversity. This allows us to find more reliably how changes of environmental factors (e.g. SST) affect biodiversity in ecosystems. The establishment of thresholds on interactions is also useful for exploring ranges of interdependencies and associated effective (alpha)-diversity with respect to the average taxonomic diversity. (beta) effective diversity is another very important macroecological indicator informing about ecosystem changes; for instance in Li and Convertino39 (beta)-diversity identified distinct ecosystem health states. However, (alpha) and (beta) (effective diversity) variability are highly linked to each other and yet looking into one or another would provide equivalent results. The difference between taxonomic and effective (beta)-diversity may provide some information about invasive or rare species have weak influence on the ecosystem since they are characterized by low TEs. Supplementary Information contains further elaborations on results. More

1.Bradford MA, Wieder WR, Bonan GB, Fierer N, Raymond PA, Crowther TW. Managing uncertainty in soil carbon feedbacks to climate change. Nat Clim Chang. 2016;6:751–8.Article 
CAS 

Google Scholar 
2.Cavicchioli R, Ripple WJ, Timmis KN, Azam F, Bakken LR, Baylis M, et al. Scientists’ warning to humanity: microorganisms and climate change. Nat Rev Microbiol. 2019;17:569–86.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
3.Crowther TW, Hoogen JVD, Wan J, Mayes MA, Keiser AD, Mo L, et al. The global soil community and its influence on biogeochemistry. Science. 2019;365:eaav0550.CAS 
PubMed 
Article 

Google Scholar 
4.Heimann M, Reichstein M. Terrestrial ecosystem carbon dynamics and climate feedbacks. Nature. 2008;451:289–92.CAS 
PubMed 
Article 

Google Scholar 
5.Wieder WR, Bonan GB, Allison SD. Global soil carbon projections are improved by modelling microbial processes. Nat Clim Chang. 2013;3:909–12.CAS 
Article 

Google Scholar 
6.Walker TWN, Kaiser C, Strasser F, Herbold CW, Leblans NIW, Woebken D, et al. Microbial temperature sensitivity and biomass change explain soil carbon loss with warming. Nat Clim Chang. 2018;8:885–9.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
7.Crowther TW, Todd-Brown KEO, Rowe CW, Wieder WR, Carey JC, Machmuller MB, et al. Quantifying global soil carbon losses in response to warming. Nature. 2016;540:104–8.CAS 
PubMed 
Article 

Google Scholar 
8.Li JQ, Pei JM, Pendall E, Fang CM, Nie M. Spatial heterogeneity of temperature sensitivity of soil respiration: a global analysis of field observations. Soil Biol Biochem. 2020;141:107675.CAS 
Article 

Google Scholar 
9.Wang QK, Zhao XC, Chen LC, Yang QP, Chen S, Zhang WD, et al. Global synthesis of temperature sensitivity of soil organic carbon decomposition: latitudinal patterns and mechanisms. Funct Ecol. 2019;33:514–23.Article 

Google Scholar 
10.Nottingham AT, Baath E, Reischke S, Salinas N, Meir P. Adaptation of soil microbial growth to temperature: using a tropical elevation gradient to predict future changes. Glob Chang Biol. 2019;25:827–38.PubMed 
PubMed Central 
Article 

Google Scholar 
11.Ye JS, Bradford MA, Dacal M, Maestre FT, García-Palacios P. Increasing microbial carbon use efficiency with warming predicts soil heterotrophic respiration globally. Glob Chang Biol. 2019;25:3354–64.PubMed 
Article 

Google Scholar 
12.Smith TP, Thomas TJH, Garcia-Carreras B, Sal S, Yvon-Durocher G, Bell T, et al. Community-level respiration of prokaryotic microbes may rise with global warming. Nat Commun. 2019;10:5124.PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
13.Schipper LA, Hobbs JK, Rutledge S, Arcus VL. Thermodynamic theory explains the temperature optima of soil microbial processes and high Q10 values at low temperatures. Glob Chang Biol 2014;20:3578–86.PubMed 
Article 

Google Scholar 
14.Pietikainen J, Pettersson M, Baath E. Comparison of temperature effects on soil respiration and bacterial and fungal growth rates. FEMS Microbiol Ecol. 2005;52:49–58.PubMed 
Article 
CAS 

Google Scholar 
15.Bárcenas-Moreno G, Gómez-Brandón M, Rousk J, Bååth E. Adaptation of soil microbial communities to temperature: comparison of fungi and bacteria in a laboratory experiment. Glob Chang Biol. 2009;15:2950–7.Article 

Google Scholar 
16.Engqvist MKM. Correlating enzyme annotations with a large set of microbial growth temperatures reveals metabolic adaptations to growth at diverse temperatures. BMC Microbiol. 2018;18:177.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
17.Oliverio AM, Bradford MA, Fierer N. Identifying the microbial taxa that consistently respond to soil warming across time and space. Glob Chang Biol. 2017;23:2117–29.PubMed 
Article 

Google Scholar 
18.Bier RL, Bernhardt ES, Boot CM, Graham EB, Hall EK, Lennon JT, et al. Linking microbial community structure and microbial processes: an empirical and conceptual overview. FEMS Microbiol Ecol. 2015;91:fiv113.PubMed 
Article 
CAS 

Google Scholar 
19.Dubey A, Malla MA, Khan F, Chowdhary K, Yadav S, Kumar A, et al. Soil microbiome: a key player for conservation of soil health under changing climate. Biodivers Conserv. 2019;28:2405–29.Article 

Google Scholar 
20.Hungate BA, Mau RL, Schwartz E, Caporaso JG, Dijkstra P, van Gestel N, et al. Quantitative microbial ecology through stable isotope probing. Appl Environ Microbiol. 2015;81:7570–81.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
21.Koch BJ, McHugh TA, Hayer M, Schwartz E, Blazewicz SJ, Dijkstra P, et al. Estimating taxon-specific population dynamics in diverse microbial communities. Ecosphere. 2018;9:e02090.Article 

Google Scholar 
22.Hamdi S, Moyano F, Sall S, Bernoux M, Chevallier T. Synthesis analysis of the temperature sensitivity of soil respiration from laboratory studies in relation to incubation methods and soil conditions. Soil Biol Biochem. 2013;58:115–26.CAS 
Article 

Google Scholar 
23.Martiny AC, Treseder K, Pusch G. Phylogenetic conservatism of functional traits in microorganisms. ISME J. 2013;7:830–8.CAS 
PubMed 
Article 

Google Scholar 
24.DeAngelis KM, Pold G, Topçuoğlu BD, van Diepen LTA, Varney RM, Blanchard JL, et al. Long-term forest soil warming alters microbial communities in temperate forest soils. Front Microbiol. 2015;6:104.PubMed 
PubMed Central 
Article 

Google Scholar 
25.Euskirchen ES, Bret-Harte MS, Shaver GR, Edgar CW, Romanovsky VE. Long-term release of carbon dioxide from Arctic Tundra ecosystems in Alaska. Ecosystems. 2017;20:960–74.CAS 
Article 

Google Scholar 
26.Reed SC, Reibold R, Cavaleri MA, Alonso-Rodríguez AM, Berberich ME, Wood TE. Chapter six—soil biogeochemical responses of a tropical forest to warming and hurricane disturbance. In: Dumbrell AJ, Turner EC, Fayle TM, editors. Advances in ecological research. (Academic Press, Cambridge MA, 2020) pp 225–52.27.Witt C, Gaunt JL, Galicia CC, Ottow JCG, Neue HU. A rapid chloroform-fumigation extraction method for measuring soil microbial biomass carbon and nitrogen in flooded rice soils. Biol Fertil Soils. 2000;30:510–9.CAS 
Article 

Google Scholar 
28.Berry D, Ben Mahfoudh K, Wagner M, Loy A. Barcoded primers used in multiplex amplicon pyrosequencing bias amplification. Appl Environ Microbiol. 2012;78:612.CAS 
PubMed Central 
Article 
PubMed 

Google Scholar 
29.Apprill A, McNally S, Parsons R, Weber L. Minor revision to V4 region SSU rRNA 806R gene primer greatly increases detection of SAR11 bacterioplankton. Aquat Micro Ecol. 2015;75:129–37.Article 

Google Scholar 
30.Parada AE, Needham DM, Fuhrman JA. Every base matters: assessing small subunit rRNA primers for marine microbiomes with mock communities, time series and global field samples. Environ Microbiol. 2016;18:1403–14.CAS 
PubMed 
Article 

Google Scholar 
31.Rohland N, Reich D. Cost-effective, high-throughput DNA sequencing libraries for multiplexed target capture. Genome Res. 2012;22:939–46.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
32.Aronesty E. ea-utils: “Command-line tools for processing biological sequencing data”. 2011. https://github.com/ExpressionAnalysis/ea-utils.33.Caporaso JG, Kuczynski J, Stombaugh J, Bittinger K, Bushman FD, Costello EK, et al. QIIME allows analysis of high-throughput community sequencing data. Nat Methods. 2010;7:335–6.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
34.Caporaso JG, Bittinger K, Bushman FD, Desantis TZ, Andersen GL, Knight R. PyNAST: a flexible tool for aligning sequences to a template alignment. Bioinformatics. 2010;26:266–7.CAS 
PubMed 
Article 

Google Scholar 
35.Edgar RC. Search and clustering orders of magnitude faster than BLAST. Bioinformatics. 2010;26:2460–1.CAS 
PubMed 
Article 

Google Scholar 
36.Quast C, Pruesse E, Yilmaz P, Gerken J, Schweer T, Yarza P, et al. The SILVA ribosomal RNA gene database project: improved data processing and web-based tools. Nucleic Acids Res. 2013;41:D590–6.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
37.Bokulich NA, Subramanian S, Faith JJ, Gevers D, Gordon JI, Knight R, et al. Quality-filtering vastly improves diversity estimates from Illumina amplicon sequencing. Nat Methods. 2013;10:57–9.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
38.Morrissey EM, Mau RL, Schwartz E, McHugh TA, Dijkstra P, Koch BJ, et al. Bacterial carbon use plasticity, phylogenetic diversity and the priming of soil organic matter. ISME J. 2017;11:1890–9.PubMed 
PubMed Central 
Article 

Google Scholar 
39.Li J, Mau RL, Dijkstra P, Koch BJ, Schwartz E, Liu X-JA, et al. Predictive genomic traits for bacterial growth in culture versus actual growth in soil. ISME J. 2019;13:2162–72.PubMed 
PubMed Central 
Article 

Google Scholar 
40.Gross N, Bagousse-Pinguet YL, Liancourt P, Berdugo M, Gotelli NJ, Maestre FT. Functional trait diversity maximizes ecosystem multifunctionality. Nat Ecol Evol. 2017;1:0132.PubMed 
PubMed Central 
Article 

Google Scholar 
41.Laliberté E, Norton DA, Scott D. Contrasting effects of productivity and disturbance on plant functional diversity at local and metacommunity scales. J Veg Sci. 2013;24:834–42.Article 

Google Scholar 
42.Plass-Johnson JG, Taylor MH, Husain AAA, Teichberg MC, Ferse SCA. Non-random variability in functional composition of coral reef fish communities along an environmental gradient. PLOS ONE. 2016;11:e0154014.PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
43.Götzenberger L, Botta-Dukát Z, Lepš J, Pärtel M, Zobel M, de Bello F. Which randomizations detect convergence and divergence in trait-based community assembly? A test of commonly used null models. J Veg Sci. 2016;27:1275–87.Article 

Google Scholar 
44.Delgado-Baquerizo M, Trivedi P, Trivedi C, Eldridge DJ, Reich PB, Jeffries TC, et al. Microbial richness and composition independently drive soil multifunctionality. Funct Ecol. 2017;31:2330–43.Article 

Google Scholar 
45.R Core Team. R: a language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing; 2020. https://www.R-project.org/.46.Bruelheide H, Dengler J, Purschke O, Lenoir J, Jiménez-Alfaro B, Hennekens SM, et al. Global trait–environment relationships of plant communities. Nat Ecol Evol. 2018;2:1906–17.PubMed 
Article 

Google Scholar 
47.Piton G, Legay N, Arnoldi C, Lavorel S, Clément J-C, Foulquier A. Using proxies of microbial community-weighted means traits to explain the cascading effect of management intensity, soil and plant traits on ecosystem resilience in mountain grasslands. J Ecol. 2020;108:876–93.CAS 
Article 

Google Scholar 
48.Alster CJ, von Fischer JC, Allison SD, Treseder KK. Embracing a new paradigm for temperature sensitivity of soil microbes. Glob Chang Biol. 2020;26:3221–9.PubMed 
Article 

Google Scholar 
49.Letunic I, Bork P. Interactive Tree Of Life (iTOL) v4: recent updates and new developments. Nucleic Acids Res. 2019;47:W256–9.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
50.Li J, Nie M, Pendall E, Reich PB, Pei J, Noh NJ, et al. Biogeographic variation in temperature sensitivity of decomposition in forest soils. Glob Chang Biol. 2020;26:1873–85.PubMed 
Article 

Google Scholar 
51.Lipson DA. The complex relationship between microbial growth rate and yield and its implications for ecosystem processes. Front Microbiol. 2015;6:615.PubMed 
PubMed Central 

Google Scholar 
52.Buckeridge KM, Mason KE, McNamara NP, Ostle N, Puissant J, Goodall T, et al. Environmental and microbial controls on microbial necromass recycling, an important precursor for soil carbon stabilization. Commun Earth Environ. 2020;1:36.Article 

Google Scholar 
53.Ali A, Yan E-R, Chang SX, Cheng J-Y, Liu X-Y. Community-weighted mean of leaf traits and divergence of wood traits predict aboveground biomass in secondary subtropical forests. Sci Total Environ. 2017;574:654–62.CAS 
PubMed 
Article 

Google Scholar 
54.Buzzard V, Michaletz ST, Deng Y, He Z, Ning D, Shen L, et al. Continental scale structuring of forest and soil diversity via functional traits. Nat Ecol Evol. 2019;3:1298–308.PubMed 
Article 

Google Scholar 
55.Luo Y-H, Cadotte MW, Burgess KS, Liu J, Tan S-L, Zou J-Y, et al. Greater than the sum of the parts: how the species composition in different forest strata influence ecosystem function. Ecol Lett. 2019;22:1449–61.PubMed 
Article 

Google Scholar 
56.Díaz S, Lavorel S, de Bello F, Quétier F, Grigulis K, Robson TM. Incorporating plant functional diversity effects in ecosystem service assessments. Proc Natl Acad Sci USA. 2007;104:20684–9.PubMed 
Article 

Google Scholar 
57.Bradford MA. Thermal adaptation of decomposer communities in warming soils. Front Microbiol. 2013;4:333.PubMed 
PubMed Central 
Article 

Google Scholar 
58.Morrissey EM, Mau RL, Schwartz E, Koch BJ, Hayer M, Hungate BA. Taxonomic patterns in the nitrogen assimilation of soil prokaryotes. Environ Microbiol. 2018;20:1112–9.CAS 
PubMed 
Article 

Google Scholar 
59.Coskun OK, Ozen V, Wankel SD, Orsi WD. Quantifying population-specific growth in benthic bacterial communities under low oxygen using H218O. ISME J. 2019;13:1546–59.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
60.Zhou G, Zhou X, Liu R, Du Z, Zhou L, Li S, et al. Soil fungi and fine root biomass mediate drought-induced reductions in soil respiration. Funct Ecol. 2020;34:2634–43.Article 

Google Scholar 
61.Melillo JM, Frey SD, Deangelis KM, Werner WJ, Bernard MJ, Bowles FP, et al. Long-term pattern and magnitude of soil carbon feedback to the climate system in a warming world. Science. 2017;358:101–5.CAS 
PubMed 
Article 

Google Scholar 
62.Johnston ASA, Sibly RM. The influence of soil communities on the temperature sensitivity of soil respiration. Nat Ecol Evol. 2018;2:1597–602.PubMed 
Article 

Google Scholar  More

We demonstrate our approach using an agricultural field under active commercial cultivation of soybeans located in the Arkansas Delta. For this study, we acquired time-series datasets of UAV and geophysical data during the 2018 growing season. This provided the opportunity to investigate the relationship between plant growth and spatially heterogeneous soil physical properties, as well as to analyse the temporal dynamics of such a relationship. All the methodology and data collection comply with relevant institutional and national legislation.Site description and crop managementThe investigated study area (34° 24.458′ N, 91° 40.462′ W), located in Humphrey, Arkansas, is a crop field of about 20 hectars that is under active commercial cultivation of soybeans (Fig. 1).Figure 1Study area located in Humphrey, Arkansas: (a) map produced by cropScape (USDA) showing part of the Delta region with the main agricultural land covers; (b) map of the field under study reporting the locations of the ground data collection; (c) map showing the topographic elevation derived by an UAV-based DEM, acquired acquisition prior to planting; (d) gSSurgo map provided by USDA, showing the soil types. Black vertical lines separate the east, center, and west sections used in the statistical analysis. Maps made with QGIS (v 3.6, https://www.qgis.org)36.Full size imageThe study area is part of the Lower Mississippi River Basin’s (LMRB) alluvial plain37. This region is characterized by a wide range of cropland, such as soybean, corn, rice, grass, and pasture, supported by a humid subtropical climate with hot, humid summers and mild, slightly drier winters.The crop field is adjacent to the oxbow Glenwood Lake (shown in Fig. 1) and is characterized by an alternating of Hebert silt-loam and Rilla silt-loam, according to the description provided by the gSSURGO map38. The Hebert profile consists of a very deep, poorly drained, moderately slowly permeable soil that formed in silty alluvium, while the Rilla profile consists of a very deep, well drained, moderately permeable soils that formed in reddish silty and loamy alluvium.High-quality soybean seed requires three appropriate conditions for germination—soil moisture, temperature and oxygen. Provided that the soil is not saturated, oxygen concentration is not a limitation, but germination will not occur in flooded soil dueto lack of oxygen. Research shows that soybeans are in general very susceptible to soil saturation and anoxic conditions, which would provoke damage to the root system and in some case the termination of the plant growth39.Planting started on April 21st in rows (raised beds) with a spacing of 0.9 m. Plants reached their first vegetative stage of emergence (VE) on the first days of May and developed their first node with fully developed unifoliate leaves around mid-May (V1). The stage of initial bloom, which represents the starting of the reproductive stage (R1–R8) after the last vegetative stage, started in mid-June, reaching the full maturity (R8) about the end of September. Subsequent to R8, the leaves dry then fall from the plant at the time of harvest, which was performed on October 24th.A common channeled (furrowed) surface irrigation approach was used as an irrigation system, which consists of a plastic pipe (polytube) located along one side of the field to deliver water into the channels between the rows40. At the site, the plastic pipe was located on the east side of the field with the rows and furrows running east to west following the main slope gradient (Fig. 1) from east to west, in order to facilitate the water flow. Smoothing of the field’s ground surface was performed prior to the planting to ensure a good surface drainage. Irrigation started on June 12th, after the rain period, and performed every 14 days.Temperature and precipitation were measured by a weather station that we installed in the southest corner of the field (Fig. 1). Based on the data collected, the average daily temperature increased from 14.6 to 24.7 °C between planting and the first vegetative stage (April–May), reached 27 °C during the vegetative period through the initial reproductive period (June–July), and decreased from 27 to 24 °C during the final period of maturity (August–September). The 2018 growing season was characterized by sporadic heavy-rain events, especially during the vegetative period (from April to mid-June), with a total cumulative precipitation (period April–September) of about 650 mm (Supplementary Figure S1). As a result of rain events, the soil moisture in April ranged from 0.17 (m3/m3) to 0.3 (m3/m3) (Supplementary Figure S2).Instrumentation, data acquisition, and preprocessingUAV-data acquisitionMultitemporal image acquisitions were performed over the field using UAV-mounted sensors. The images were acquired on: May 8th, a few days after planting, May 28th, during the early vegetative development, and June 25th, July 23rd, September 19th, the beginning, middle, and end of the reproductive period.We used the DJI Matrice 600 flying platform with the DJI FC350 on-board RGB (red, green, blue) camera as the primary image acquisition sensor. A flight plan to acquire high-resolution images at centimeter resolution was designed using this acquisition system. Flights were performed during the maximum Sun’s elevation and with clear sky to avoid possible changes in light conditions, at an altitude of about 60 m, with a 75% image overlap, reaching a ground resolution of 2 cm.In order to minimize the impact of possible variations in illumination conditions during each acquisition, two normalized mosaics were computed by using the pictures of a reference panel taken before and after each flight. The final mosaic was obtained by averaging the two normalized mosaics.A total of 13 Ground Control Points (GCP) were positioned along the field perimeter and within the field for the mosaic reconstruction (Fig. 1). These consisted of flat targets and were surveyed at each acquisition with a Topcon real‐time kinematic differential global positioning system (Hyper-V RTK-GPS). The system is composed of a rover and a fixed base station, ensuring high-resolution geolocation information at centimeter resolution. All the collected geolocation data were post-processed obtaining an average planar error of about 2 cm and a vertical error of about 3 cm.We used the commercial PhotoScan (Agisoft) software for photogrammetry to construct the final mosaics. The software uses photogrammetry techniques based on structure from motion from multi-view stereo (SfM-MVS) images42. The processing produced mosaics with spatial resolution less than 10 cm, but the mosaics used in this work have a have a pixel size of 10 cm. Additionally, we used the software to compute digital elevation models (DEMs) from the RGB images with a final vertical resolution less than 5 cm.Electromagnetic Induction (EMI) systemSoil ECa was measured with an electromagnetic induction (EMI) system (CMD Mini-Explorer, GF Instruments), which consists of a transmitter and three receiver coils, allowing an effective acquisition of soil apparent electrical conductivity (ECa) averaged across a depth range of 0.5 m, 1.0 m, and 1.8 m from the surface. The instrument provides soil ECa measurements in milliSiemes per meter (mS/m) and the in-phase in parts per thousand (not considered in this study). The system was mounted in a sled pulled by an all-terrain vehicle and supported by a Differential Global Positioning System (DGPS) in order to tag each measurement with high-resolution positioning information. Time-lapse soil ECa data collection was performed four times, including before and during the growing season: April 21st, May 12th , June 11th, and July 9th. EMI acquisitons were planned according to the irrigation plan and UAV acquisitions. Since the EMI instrument is sensitive to changes in soil moisture, we avoided the collection during or right after irrigation events. The data were acquired consistently every 10 rows with a distance between traverses of about 9 m.In situ soil monitoring, ground imagery, and soil samplingSoil sensors (5TE, METER) were used for continuous acquisition of volumetric water content (VWC) and soil temperature. The sensors were deployed in five areas (A, B …, E, see Fig. 1) based on two factors: a) spatial varibaility of geophysical properties, b) differences in yield during the 2017 growing season (historical data available for this study). Each area has one logger that accomodates 4 sensors within a maximum radius of 4 m. To capture both shallower and deeper soil dynamics, we positioned 2 sensors vertically aligned at the depths of 12 cm and 25 cm on both the south (side1) and north (side2) of the logger. The areas A and D were positioned in high and low soil EC zones, respectively, while B was in an area of transition. Sensors A, B, and D were positioned along the same rows to forms two transects, one on the south and one on the north. Sensor E was positioned within a a very low 2017 yield in the west side of the field but at a mid-range of soil EC values. The logger C was positioned across a soil EC boundarie, but it stopped recording data at early season and was discarded from this analysis. The collected temprature data are used to correct the EMI acquisitions, whereas soil moisture is used to track the irrigation events and possible saturations. As a support system for the sensors’ data transmission, we used EM60G (METER) data loggers and ZENTRA cloud (METER) for cloud data storage.Ground RGB images were taken at several plots (Fig. 1) for ground-based plant spatial abundance assessment. We used a wooden-frame of 1 m by 1 m to delimit the plot’s area and took a picture from above (1.4 m from the ground). This assessment will be used to validate UAV-based plant spatial abundance estimates.Soil samples were collected at twelve locations (Fig. 1). Physical laboratory testing was performed at the Fayetteville Agriculture Diagnostic Laboratory of the University fo Arkanas to characterize soil properties, such as clay density, porosity, dry bulk density, volumetric water content, and organic matter. Such information was used for statistical analysis and to provide an interpretation of the soil ECa measurements. No plant tissues or seeds were collected in this study.Yield data from a combine harvesterCrop yield was harvested on October 24th. The combine harvester recorded yield information as cloud-point data in bushels per acre (bu/ac) and reported here in kilogram per hectare (kg/ha), and the relativre geolocation information. The combine harvester is able to collect 11 rows at once and records a point measurement every 2 s. Areas with artefacts were determined based on the swath parameter recorded during the harvesting. The wrong choice of the swath parameter on the combine harvester would produce unreliable yield values. Such areas were not considered in further analysis.MethodologyWe designed a data analysis framework to combine the UAV optical images and soil geophysical measurements for characterizing plant phenological and physiological properties and soil spatial heterogeneity, respectively. We then performed statistical analysis of their co-variability to quantify the influence of soil properties on plant development. This framework is constructed by the following three steps:

1.

UAV-based monitoring of plant dynamics: This step presents the data processing pipeline used for computing UAV products to characterize plant phenology and structure, such as plant spatial abundance, plant-specific vigor, and plant height.

2.

EMI-based monitoring of soil characterization: This step focuses on quantifying the soil spatial variability. We perform statistical analyses to identify the relationship between soil ECa signal and textural information derived by soil samples collected during the ground data acquisition.

3.

Quantifying soil–plant spatiotemporal co-variability: This step focuses on the statistical analysis of plant physiological and structural properties and near-surface properties, as well as identifying key environmental factors.

UAV-based plant monitoringFigure 2 depicts the data processing pipeline developed to assess the spatial variability of plant characteristics during the growing season. The steps to obtain the plant spatial abundance and plant vigor maps are presented in detail in the following sections.Figure 2Pipeline developed for the UAV data processing.Full size imagePlant spatial abundancePlant spatial abundance is estimated by the RGB mosaics. In the first step of the procedure, a vegetation map was obtained by performing a binary classification to discriminate the two classes of vegetation and soil. We used a supervised spectral-spatial image classification approach43 that combines both the spectral and the contextual (i.e., geometrical) information. The contextual information represents the pixel spatial arrangement and the spatial reletionship between adjacent pixels. By modelling such infromation, we can extract structures within the scene. For this step, we use mathematical operators, such as attibute profiles44,45,46 that belong to the image processing sub-field of mathematical morpohlogy. Such operators are 2D filters that act on regions of connected pixels based on the evaluation of an attribute (e.g., scale, defined as number of pixel composing the region). The advantage of using this operator is the ability to preserve the geometrical detail of the structures in the scene. In a filtered image (or feature), regions with similar properties (i.e., scale) are preserved or otherwise merged to their surroundings (i.e., filtered). A multiscale contextual characterization is obtained by applying a filter recursively and relaxing the scale parameter of the filter at each iteration, resulting in a stack of filtered images. An automatic procedure developed in45 is used to perform such a characterization. Such operation serves to extract homogeneous structures present in the scene at different spatial scales, allowing us to better delineate the soybean rows. The original RGB image, together with the filtered images, compose the feature space used as input to a support vector machine (SVM) algorithm with a radial basis function (RBF) kernel. The algorithm is based on the LIBSVM library47 developed for the MATLAB environment, using a one-against-one multiclass strategy. In a second step, we compute the spatial abundance of the detected plants as the percentage of pixels that belong to the vegetation class within a grid unit. In our study we chose a grid unit of 2 by 2 m, which would allow to obtain a good approximation of the local change, while capturing large scale patterns. The quality of the plant spatial abundance map is mainly dependent on the performance of the classification algorithm used to separate the plants from the soil in the background. To validate this product, we computed the classification performance based on ground-truth data and compared the UAV-based estimates of plant abundance to ground-based assessments. The ground assessment of the plant spatial abundance was performed by selecting 31 plots in conjunction of the UAV acquisition occurring in May 28th, 2018. The plots were located mainly along the east and west sides of the field within a distance that ranges between 12 and 25 m from the edge of the field towards the field center; and some plots were collected along the south and north sides within a distance of 50 m from the field edge. These plots were selected to captured representative spatial abundances, ranging from 0% (almost bare) to 60%. Ground images were systematically taken at a height of 1.4 m from the ground. From these images, we computed the GCC and identified a threshold to separate the plant covered area from the soil backgroun. The plant spatial abundance was then computed as the percentage of pixels identified as soybean over the entire area of the plot. All the ground-data were geolocated using the RTK-GPS system previously described.Plant-specific vigor estimationRemote-sensing-derived VIs have been extensively used as a proxy to estimate plant vigor and investigate plant physiology and response to possible ecosystem changes. Widely used VIs for plant characterization are in general derived by the near-infrared (NIR) channel (e.g., NDVI, SAVI, etc.), which is sensitive to changes in chlorophyll and leaf area. However, RGB-based indices have been developed and used in several phenology studies, showing their effectiveness in estimating plant characteristics compared to NIR-based VIs48. The GCC was chosen based on authors’ recent studies, showing that this RGB-based VI is less affected by saturation compared to NDVI31. GCC is formally defined as follows:$${text{GCC}} = frac{{{text{G}}_{{{text{DN}}}} }}{{{text{B}}_{{{text{DN}}}} + {text{G}}_{{{text{DN}}}} + {text{R}}_{{{text{DN}}}} }},$$
(1)

where R, G, and B are the red, green, and blue channels, and DN refers to the pixel intensity values as digital numbers41. The plant-specific vigor was computed considering only those pixels identified as vegetation, the information for which is provided by the binary classification. In order to have products spatially comparable, the plant-specific vigor was transformed into a 2 by 2 m grid by computing the average within each grid unit. This procedure has the advantage to minimize the soil component in estimating plant vigor, avoiding the nonuniqueness of lower productivity versus increased soil coverage. Specifically, it allows a more accurate way of capturing the spatial variability of the plant’s vigor, in particular during the early vegetative stage, when the soil component is overabundant compared to the plant spatial abundance. This is a clear advantage over a satellite-derived VI with a coarser resolution that would provide an underestimation of plant vigor due to soil coverage.Plant height estimationPlant height was computed by subtracting the digital elevation model (DEM) of bare soil (i.e., reference DEM) from the digital surface model computed at each acquisition. The reference DEM was computed from acquisition occurred on May 8th, which was planned a week after planting so that it is not affected by surface disturbances by the planter. DEMs were computed using the PhotoScan software (Agisoft) during the mosaic reconstruction phase. We assume that the ground surface elevation stays the same during the growing season.EMI-based soil characterizationWe used the ECa data to evaluate the spatial heterogeneity of soil properties and the temporal variability in soil moisture. We primarily used the EMI measurements associated with the 1 m effective depth, which is the depth range expected to encompass the soybean roots. Temperature corrections were applied to the point-cloud data for correcting the measured soil ECa to the reference temperature of 25 °C. We employed the exponential model presented in Corwin and Lesch49 and developed by Sheets and Hendrickx50, and formally described as follows:$${text{E}}C_{25} = {text{E}}C_{a} cdot left[ {0.4470 + 1.4034e^{{left( { – T/26.815} right)}} } right],$$
(2)

where ({text{E}}C_left{ {25} right}) represents the corrected soil ECa at the reference temperature of 25 °C, ({text{E}}C_a) represents the measured soil ECa, and ({text{T}}) represents the soil temperature at the time of the acquisition. The model was chosen based on the on the methodological comparison51, which shows that the exponential model provides very similar estimations compared to the widely used ratio model31,52 in the temperature range of 3–43 °C, but with a smaller total residual.A spatial correction was applied to the data to resolve a spatial shift caused by the distance from the EMI instrument and the actual position of the GPS antenna, which were 4.5 m apart due to the system setup. We performed statistical analysis between EMI data and soil samples collected during the April data acquisition to identify the main control on the spatial variability of soil ECa.Co-variability among plant and soil signaturesWe performed a suite of statistical analyses, based on the high-resolution data layers, on soil and plant characteristics (i.e., plant spatial abundance, plant vigor, plant height, soil ECa, and crop yield) to investigate the impact of soil properties on plant development during the growing season. For the analysis, we considered the yield point-cloud data as geo-reference.Since the EMI and crop-yield data had a lower vertical resolution (along the y axis) of about 10 m (distance between rows selected for the data collection), we computed the average of each metric within a window of 20 m by 20 m. The exploratory data analysis included run charts and scatter plots to investigate the temporal evolution of plant development during the growing season, and the co-variability between plant characteristics and soil ECa. We evaluated the correlation between variables using Pearson’s correlation, as well as their spatial association. To control for the spatial non-independence, we derived a corrected Perason’s correlation by using the hypothesis testing procedure propsed by Clifford, Richardson, and Hémon (1989)53. The method uses Moran’s I54 to compute the spatial autocorrelation in the spatial data sets and estimate the correct degrees of freedom, which are then used to assess the significance of the correlation. The method uses the Sturges’ formula to identify the number of distance classes55. More