Ecology

1.
Faust K, Raes J, Faust K, Raes J. Microbial interactions: from networks to models. Nat Rev Microbiol. 2012;10:538–50.
CAS  PubMed  Article  Google Scholar 
2.
Phelan VV, Liu WT, Pogliano K, Dorrestein PC. Microbial metabolic exchange—the chemotype-to-phenotype link. Nat Chem Biol. 2011;8:26–35.
PubMed  Article  CAS  Google Scholar 

3.
Natale P, Brüser T, Driessen AJM. Sec- and Tat-mediated protein secretion across the bacterial cytoplasmic membrane: Distinct translocases and mechanisms. Biochim Biophys Acta. 2007;1778:1735–56.
PubMed  Article  CAS  Google Scholar 

4.
Holland IB. The extraordinary diversity of bacterial protein secretion mechanisms. Meth Mol Biol. 2010;619:1–20.
CAS  Article  Google Scholar 

5.
Guerrero-Mandujano A, Hernández-Cortez C, Ibarra JA, Castro-Escarpulli G. The outer membrane vesicles: Secretion system type zero. Traffic. 2017;18:425–32.
CAS  PubMed  Article  Google Scholar 

6.
Orench‐Rivera N, Kuehn MJ. Environmentally controlled bacterial vesicle‐mediated export. Cell Microbiol. 2016;18:1525–36.
PubMed  PubMed Central  Article  CAS  Google Scholar 

7.
Kim JH, Lee J, Park J, Gho YS, editors. Gram-negative and Gram-positive bacterial extracellular vesicles. Semin Cell Dev Biol. 2015;40:97–104.

8.
Schwechheimer C, Kuehn MJ. Outer-membrane vesicles from Gram-negative bacteria: biogenesis and functions. Nat Rev Microbiol. 2015;13:605–19.
CAS  PubMed  PubMed Central  Article  Google Scholar 

9.
McBroom AJ, Kuehn MJ. Release of outer membrane vesicles by Gram‐negative bacteria is a novel envelope stress response. Mol Microbiol. 2007;63:545–58.
CAS  PubMed  PubMed Central  Article  Google Scholar 

10.
Arntzen MO, Varnai A, Mackie RI, Eijsink VGH, Pope PB. Outer membrane vesicles from Fibrobacter succinogenes S85 contain an array of carbohydrate-active enzymes with versatile polysaccharide-degrading capacity. Environ Microbiol. 2017;19:2701–14.
CAS  PubMed  Article  Google Scholar 

11.
Nordstrom T, Blom AM, Tan TT, Forsgren A, Riesbeck K. Ionic binding of C3 to the human pathogen Moraxella catarrhalis is a unique mechanism for combating innate immunity. J Immunol. 2005;175:3628–36.
PubMed  Article  Google Scholar 

12.
Fulsundar S, Harms K, Flaten GE, Johnsen PJ, Chopade B, Nielsen KM. Gene transfer potential of outer membrane vesicles of Acinetobacter baylyi and effects of stress on vesiculation. Appl Environ Microb. 2014;80:3469–83.
Article  CAS  Google Scholar 

13.
Mashburn LM, Whiteley M. Membrane vesicles traffic signals and facilitate group activities in a prokaryote. Nature. 2005;437:422–5.
CAS  PubMed  Article  Google Scholar 

14.
Toyofuku M, Morinaga K, Hashimoto Y, Uhl J, Shimamura H, Inaba H, et al. Membrane vesicle-mediated bacterial communication. ISME J. 2017;11:1504–9.
PubMed  PubMed Central  Article  Google Scholar 

15.
Lee EY, Choi DY, Kim DK, Kim JW, Park JO, Kim S, et al. Gram‐positive bacteria produce membrane vesicles: proteomics‐based characterization of Staphylococcus aureus‐derived membrane vesicles. Proteomics. 2009;9:5425–36.
CAS  PubMed  Article  Google Scholar 

16.
Prados-Rosales R, Baena A, Martinez LR, Luque-Garcia J, Kalscheuer R, Veeraraghavan U, et al. Mycobacteria release active membrane vesicles that modulate immune responses in a TLR2-dependent manner in mice. J Clin Investig. 2011;121:1471–83.
PubMed  Article  CAS  Google Scholar 

17.
Prados-Rosales R, Brown L, Casadevall A, Montalvo-Quiros S, Luque-Garcia JL. Isolation and identification of membrane vesicle-associated proteins in Gram-positive bacteria and mycobacteria. MethodsX. 2014;1:124–9.
PubMed  PubMed Central  Article  Google Scholar 

18.
White DW, Elliott SR, Odean E, Bemis LT, Tischler AD. Mycobacterium tuberculosis Pst/SenX3-RegX3 regulates membrane vesicle production independently of ESX-5 activity. mBio. 2018;9:e00778–18.
CAS  PubMed  PubMed Central  Google Scholar 

19.
Hoffmann C, Leis A, Niederweis M, Plitzko JM, Engelhardt H. Disclosure of the mycobacterial outer membrane: cryo-electron tomography and vitreous sections reveal the lipid bilayer structure. PNAS. 2008;105:3963–7.
CAS  PubMed  Article  Google Scholar 

20.
Ganz T, Nemeth E. Iron homeostasis in host defence and inflammation. Nat Rev Immunol. 2015;15:500–10.
CAS  PubMed  PubMed Central  Article  Google Scholar 

21.
Huber DL. Synthesis, properties, and applications of iron nanoparticles. Small. 2005;1:482–501.
CAS  PubMed  Article  Google Scholar 

22.
Wandersman C, Delepelaire P. Bacterial iron sources: from siderophores to hemophores. Annu Rev Microbiol. 2004;58:611–47.
CAS  PubMed  Article  Google Scholar 

23.
Morel FM, Price N. The biogeochemical cycles of trace metals in the oceans. Science. 2003;300:944–7.
CAS  PubMed  Article  PubMed Central  Google Scholar 

24.
Ram RJ, VerBerkmoes NC, Thelen MP, Tyson GW, Baker BJ, Blake RC, et al. Community proteomics of a natural microbial biofilm. Science. 2005;308:1915–20.
CAS  PubMed  Article  Google Scholar 

25.
Cao B, Shi L, Brown RN, Xiong Y, Fredrickson JK, Romine MF, et al. Extracellular polymeric substances from Shewanella sp. HRCR‐1 biofilms: characterization by infrared spectroscopy and proteomics. Environ Microbiol. 2011;13:1018–31.
CAS  PubMed  Article  Google Scholar 

26.
Vong L, Laës A, Blain S. Determination of iron–porphyrin-like complexes at nanomolar levels in seawater. Anal Chim Acta. 2007;588:237–44.
CAS  PubMed  Article  Google Scholar 

27.
Létoffé S, Nato F, Goldberg ME, Wandersman C. Interactions of HasA, a bacterial haemophore, with haemoglobin and with its outer membrane receptor HasR. Mol Microbiol. 1999;33:546–55.
PubMed  Article  Google Scholar 

28.
Tong Y, Guo M. Bacterial heme-transport proteins and their heme-coordination modes. Arch Biochem Biophys. 2009;481:1–15.
CAS  PubMed  Article  Google Scholar 

29.
Pilpa RM, Robson SA, Villareal VA, Wong ML, Phillips M, Clubb RT. Functionally distinct NEAT (NEAr Transporter) domains within the Staphylococcus aureus IsdH/HarA protein extract heme from methemoglobin. J Biol Chem. 2009;284:1166–76.
CAS  PubMed  PubMed Central  Article  Google Scholar 

30.
Gat O, Zaide G, Inbar I, Grosfeld H, Chitlaru T, Levy H, et al. Characterization of Bacillus anthracis iron‐regulated surface determinant (Isd) proteins containing NEAT domains. Mol Microbiol. 2008;70:983–99.
CAS  PubMed  PubMed Central  Google Scholar 

31.
Choby JE, Skaar EP. Heme synthesis and acquisition in bacterial pathogens. J Mol Biol. 2016;428:3408–28.
CAS  PubMed  PubMed Central  Article  Google Scholar 

32.
Allen CE, Schmitt MP. HtaA is an iron-regulated hemin binding protein involved in the utilization of heme iron in Corynebacterium diphtheriae. J Bacteriol. 2009;191:2638–48.
CAS  PubMed  PubMed Central  Article  Google Scholar 

33.
Allen CE, Schmitt MP. Novel hemin binding domains in the Corynebacterium diphtheriae HtaA protein interact with hemoglobin and are critical for heme iron utilization by HtaA. J Bacteriol. 2011;193:5374–85.
CAS  PubMed  PubMed Central  Article  Google Scholar 

34.
Duckworth AW, Grant S, Grant WD, Jones BE, Meijer D. Dietzia natronolimnaios sp. nov., a new member of the genus Dietzia isolated from an East African soda lake. Extremophiles. 1998;2:359–66.
CAS  PubMed  Article  PubMed Central  Google Scholar 

35.
Mayilraj S, Suresh K, Kroppenstedt R, Saini H. Dietzia kunjamensis sp. nov., isolated from the Indian Himalayas. Int J Syst Evol Microbiol. 2006;56:1667–71.
CAS  PubMed  Article  PubMed Central  Google Scholar 

36.
Li J, Chen C, Zhao G-Z, Klenk H-P, Pukall R, Zhang Y-Q, et al. Description of Dietzia lutea sp. nov., isolated from a desert soil in Egypt. Syst Appl Microbiol. 2009;32:118–23.
CAS  PubMed  Article  PubMed Central  Google Scholar 

37.
Fang H, Qin X-Y, Zhang K-D, Nie Y, Wu X-L. Role of the Group 2 Mrp sodium/proton antiporter in rapid response to high alkaline shock in the alkaline-and salt-tolerant Dietzia sp. DQ12-45-1b. Appl Microbiol Biotechnol. 2018;102:3765–77.
CAS  PubMed  Article  PubMed Central  Google Scholar 

38.
Wang X-B, Chi C-Q, Nie Y, Tang Y-Q, Tan Y, Wu G, et al. Degradation of petroleum hydrocarbons (C6–C40) and crude oil by a novel Dietzia strain. Bioresour Technol. 2011;102:7755–61.
CAS  PubMed  Article  PubMed Central  Google Scholar 

39.
Rédei GP M9 Bacterial Minimal Medium. In: Rédei GP, editors. Encyclopedia of genetics, genomics, proteomics and informatics, 3rd edn. Dordrecht: Springer Group; 2008. pp. 484–6.

40.
Van Kessel JC, Hatfull GF. Recombineering in Mycobacterium tuberculosis. Nat Methods. 2007;4:147–52.
PubMed  Article  CAS  Google Scholar 

41.
Liang J, Jiangyang J, Nie Y, Wu X. Regulation of the alkane hydroxylase CYP153 gene in a Gram-positive alkane-degrading bacterium, Dietzia sp. strain DQ12-45-1b. Appl Environ Microbiol. 2016;82:608–19.
CAS  PubMed  PubMed Central  Article  Google Scholar 

42.
Lu S, Nie Y, Tang Y-Q, Xiong G, Wu X-L. A critical combination of operating parameters can significantly increase the electrotransformation efficiency of a Gram-positive Dietzia strain. J Microbiol Methods. 2014;103:144–51.
CAS  PubMed  Article  Google Scholar 

43.
Szvetnik A, Bihari Z, Szabo Z, Kelemen O, Kiss I. Genetic manipulation tools for Dietzia spp. J Appl Microbiol. 2010;109:1845–52.
CAS  PubMed  Google Scholar 

44.
Deininger PL. Molecular cloning: a laboratory manual. Anal Biochem. 1990;186:182–3.
Article  Google Scholar 

45.
McBroom AJ, Johnson AP, Vemulapalli S, Kuehn MJ. Outer membrane vesicle production by Escherichia coli is independent of membrane instability. J Bacteriol. 2006;188:5385–92.
CAS  PubMed  PubMed Central  Article  Google Scholar 

46.
Prados-Rosales R, Weinrick BC, Pique DG, Jacobs WR Jr, Casadevall A, Rodriguez GM. Role for Mycobacterium tuberculosis membrane vesicles in iron acquisition. J Bacteriol. 2014;196:1250–6.
PubMed  PubMed Central  Article  CAS  Google Scholar 

47.
Bligh EG, Dyer WJ. A rapid method of total lipid extraction and purification. Biochem Cell Biol. 1959;37:911–7.
CAS  Google Scholar 

48.
Keddie RM, Cure GL. The cell wall composition and distribution of free mycolic acids in named strains of coryneform bacteria and in isolates from various natural sources. J Appl Microbiol. 1977;42:229–52.
CAS  Google Scholar 

49.
Liu Y, Zhang Q, Hu M, Yu K, Fu J, Zhou F, et al. Proteomic analyses of intracellular Salmonella enterica serovar Typhimurium reveal extensive bacterial adaptations to infected host epithelial cells. Infect Immun. 2015;83:2897–906.
CAS  PubMed  PubMed Central  Article  Google Scholar 

50.
Calderoncelis F, Encinar JR, Sanzmedel A. Standardization approaches in absolute quantitative proteomics with mass spectrometry. Mass Spectrom Rev. 2018;37:715–37.
CAS  Article  Google Scholar 

51.
Liang J-L, Gao Y, He Z, Nie Y, Wang M, JiangYang J-H, et al. Crystal structure of TetR family repressor AlkX from Dietzia sp. strain DQ12-45-1b implicated in biodegradation of n-alkanes. Appl Environ Microbiol. 2017;83:e01447–17.
CAS  PubMed  PubMed Central  Article  Google Scholar 

52.
Tashiro Y, Hasegawa Y, Shintani M, Takaki K, Ohkuma M, Kimbara K, et al. Interaction of bacterial membrane vesicles with specific species and their potential for delivery to target cells. Front Microbiol. 2017;8:571.
PubMed  PubMed Central  Article  Google Scholar 

53.
Marchler-Bauer A, Derbyshire MK, Gonzales NR, Lu S, Chitsaz F, Geer LY, et al. CDD: NCBI’s conserved domain database. Nucleic Acids Res. 2014;43:D222–D6.
PubMed  PubMed Central  Article  CAS  Google Scholar 

54.
Tamura K, Stecher G, Peterson D, Filipski A, Kumar S. MEGA6: molecular evolutionary genetics analysis version 6.0. Mol Biol Evol. 2013;30:2725–9.
CAS  PubMed  PubMed Central  Article  Google Scholar 

55.
Yu NY, Wagner JR, Laird MR, Melli G, Rey S, Lo R, et al. PSORTb 3.0: improved protein subcellular localization prediction with refined localization subcategories and predictive capabilities for all prokaryotes. Bioinformatics. 2010;26:1608–15.
CAS  PubMed  PubMed Central  Article  Google Scholar 

56.
Song H, Sandie R, Wang Y, Andrade-Navarro MA, Niederweis M. Identification of outer membrane proteins of Mycobacterium tuberculosis. Tuberculosis. 2008;88:526–44.
CAS  PubMed  Article  Google Scholar 

57.
Daffé M, Quémard A, Marrakchi H. Mycolic acids: from chemistry to biology. In: Geiger O, editors. Biogenesis of fatty acids, lipids and membranes. Cham: Springer; 2017. p. 1–36.

58.
Choi D, Kim D, Choi SJ, Lee J, Choi J, Rho S, et al. Proteomic analysis of outer membrane vesicles derived from Pseudomonas aeruginosa. Proteomics. 2011;11:3424–9.
CAS  PubMed  Article  Google Scholar 

59.
Marchand CH, Salmeron C, Bou Raad R, Meniche X, Chami M, Masi M, et al. Biochemical disclosure of the mycolate outer membrane of Corynebacterium glutamicum. J Bacteriol. 2012;194:587–97.
CAS  PubMed  PubMed Central  Article  Google Scholar 

60.
Daffe M, Marrakchi H. Unraveling the structure of the mycobacterial envelope. Microbiol Spectr. 2019;7:1087–95.
Article  Google Scholar 

61.
Nishiuchi Y, Baba T, Yano I. Mycolic acids from Rhodococcus, Gordonia, and Dietzia. J Microbiol Methods. 2000;40:1–9.
CAS  PubMed  Article  Google Scholar 

62.
Collins M, Goodfellow M, Minnikin D. A survey of the structures of mycolic acids in Corynebacterium and related taxa. Microbiology. 1982;128:129–49.
CAS  Article  Google Scholar 

63.
Rath P, Saurel O, Czaplicki G, Tropis M, Daffé M, Ghazi A, et al. Cord factor (trehalose 6, 6′-dimycolate) forms fully stable and non-permeable lipid bilayers required for a functional outer membrane. Biochim Biophys Acta-Biomemb. 2013;1828:2173–81.
CAS  Article  Google Scholar 

64.
Caruana JC, Walper SA. Bacterial membrane vesicles as mediators of microbe – microbe and microbe – host community interactions. Front Microbiol. 2020;11:432.
PubMed  PubMed Central  Article  Google Scholar 

65.
Rich PR, Maréchal A 8.5 electron transfer chains: structures, mechanisms and energy coupling. In: Egelman EH, editor. Comprehensive biophysics. Amsterdam: Elsevier; 2012. p. 72–93.

66.
Butaitė E, Baumgartner M, Wyder S, Kümmerli R. Siderophore cheating and cheating resistance shape competition for iron in soil and freshwater Pseudomonas communities. Nat Commun. 2017;8:1–12.
Article  CAS  Google Scholar 

67.
Zuber B, Chami M, Houssin C, Dubochet J, Griffiths G, Daffé M. Direct visualization of the outer membrane of mycobacteria and corynebacteria in their native state. J Bacteriol. 2008;190:5672–80.
CAS  PubMed  PubMed Central  Article  Google Scholar 

68.
Sani M, Houben ENG, Geurtsen J, Pierson J, De Punder K, Van Zon M, et al. Direct visualization by cryo-EM of the mycobacterial capsular layer: a labile structure containing ESX-1-secreted proteins. PLoS Pathog. 2010;6:e1000794.
PubMed  PubMed Central  Article  CAS  Google Scholar 

69.
Kramer J, Özkaya Ö, Kümmerli R. Bacterial siderophores in community and host interactions. Nat Rev Microbiol. 2020;18:152–63.
CAS  PubMed  Article  Google Scholar 

70.
Rakoff-Nahoum S, Coyne MJ, Comstock LE. An ecological network of polysaccharide utilization among human intestinal symbionts. Curr Biol. 2014;24:40–9.
CAS  PubMed  Article  PubMed Central  Google Scholar  More

1.
Food and Agriculture Authority. Microbiome the Missing Link. http://www.fao.org/3/ca6767en/CA6767EN.pdf (2020).
2.
Singh, B. K. & Trivedi, P. Microbiome and the future for food and nutrient security. Microbial. Biotechnol. 10, 50–53 (2017).
Article  Google Scholar 

3.
National Academies of Science, Engineering and Medicine. New Report Identifies Five Breakthroughs to Address Urgent Challenges and Advance Food and Agricultural Sciences by 2030. http://www8.nationalacademies.org/onpinews/newsitem.aspx?RecordID=25059 (2018).

4.
Stratistics Market Research Consulting. Agricultural Microbials – Global Market Outlook (2017–2026). https://www.researchandmarkets.com/research/vdth4t/global?w=4 (2020).

5.
European Commission. Farm to fork fact sheet https://ec.europa.eu/commission/presscorner/api/files/attachment/865559/factsheet-farm-fork_en.pdf.pdf (2020).

6.
Trivedi, P., Leach, J. E. & Tringe, S. G. et al. Plant–microbiome interactions: from community assembly to plant health. Nat. Rev. Microbiol. https://doi.org/10.1038/s41579-020-0412-1 (2020).
Article  PubMed  Google Scholar 

7.
Qiu, Z. et al. New frontiers in agriculture productivity: optimised microbial inoculants and in situ microbiome engineering. Biotechnol. Adv. 37, 107371 (2019).
CAS  Article  Google Scholar 

8.
Kaminsky, L. M. et al. The inherent conflicts in developing soil microbial noculants. Trends Biotechnol. 37, 140–151 (2019).
CAS  Article  Google Scholar 

9.
Busby, P. E. et al. Research priorities for harnessing plant microbiomes in sustainable agriculture. PLOS Biol. 15, e2001793 (2017).
Article  Google Scholar  More

The study area
The study was carried out on Egbokodo River (longitude 5° 38′ and 5° 41′ and latitude 5° 36′ and 5° 33′) in Warri South Local Government Area of Delta State, Southern Nigeria (Fig. 1), between the periods of September 2008–May 2009. The river is a brackish and tidal River that serves as a source of water for drinking, washing, and fishing to the communities in the catchment area. Three (3) Stations (tagged Stations A, B, and C) were selected about 150 m apart, based on distinct anthropogenic activities. Station A was located at a vandalized oil pipeline, while Stations B and C were located downstream at points of dredging and municipal waste disposal respectively. Station A was 6.3–9.3 m in depth, Station B was 11.4–16.1, and Station C was 7.5–19.5 m during the study duration.
Figure 1

Map of the study area showing sampled stations. Map designed using QGIS software version 3.10.1 ‘A Coruña’ (QGIS Development Team29). https://qgis.org/en/site/forusers/download.html#.

Full size image

The study area comprises of coarse and interspersed soil with lignite and patches of laterite and sandy clay soil. The climate of the study area is typically tropical. It is characterized by the humid tropical wet and dry climate which is primarily regulated by rainfall. The wet season lasts a period of 7 months (April to October). Rainfall ranged from 15 to 91 mm during this period. The driest months are December to January; with a mean monthly rainfall of 15 mm. The bank of the river was densely shaded by a thick canopy of vegetation, dominated by mangrove plants, Nypa palm, and Rhizophora sp.
Collection of samples (water and zooplankton)
Water samples were collected from the 3 stations using a 1 L sampling bottle which was pre-cleaned with the deionized water at each station. This sampling procedure was repeated monthly from September 2008 to May 2009. The samples were preserved in a cooler and transported to the laboratory where they were refrigerated at − 10 °C before the physiochemical analysis. Preservation and analysis of water samples were according to standard methods of the American Public Health Association (APHA).
Samples of zooplankton were collected at the 3 stations between 0800 and 1100 h by towing a hydrobios plankton net (mesh size 25 µm) with a speed boat at 2 knots, just below the water surface for 5 min at every station. At each station, the filtered zooplankton samples were condensed in a 25 mL plankton bottle and preserved using buffered 4% formalin. Each plankton bottle was properly labeled indicating the stations and dates of collection. This procedure was repeated for 9 months (September 2008–May 2009).
Analysis of water
Determination of pH
The pH was estimated using a PH meter—Orion Model 290A (ASTM D 1293B) and recorded accordingly every month.
Measurement of temperature (°C)
A mercury-in-glass thermometer was used to measure surface water temperature. A stable initial reading was ensured by shaking it the thermometer carefully. Afterward, the thermometer was left inside the water for about 3 min till a stable reading was observed and recorded.
Determination of phosphate
Five (5) mL antimony molybdate was added to 40 mL of water sample was in a 50 mL measuring cylinder. Afterward, 2 mL of Ascorbic acid was added to the mixture. It was left to stand for 30 min for full colour formation2. The absorbance was measured with a UV–visible spectrophotometer at 680 nm.
Phosphate was then calculated thus;

$${text{Phosphate}},({text{mg}}/{text{l}}) = frac{{{text{Y}} – {text{C}}}}{{text{M}}}$$
(1)

In Eq. (1) above, Y = absorbance of the sample.
C = absorbance of blank

$${text{M}} = {text{Gradient}}frac{{({text{B}} – {text{A}})}}{{text{X}}}$$
(2)

B = absorbance of standard (Eq. 2)
A = absorbance of blank
X = concentration of the standard.
Determination of nitrate
Nitrate was tested using the diazotization method—Alpha 419 C/ASTM D3867. 0.5 mL of (0.1% W/V) NaN3 was added to the water sample to remove any NO2 present. 3.0 mL of (2.6% W/V) NH4Cl solution was added. One (1) mL of (2.1% W/V) Borax solution was added. 0.5–0.6 g of spongy cadmium was added. It was then covered and shaken for some 15 to 20 min. Afterward, 7 mL of the solution was transferred to a 25 mL measuring cylinder. 1 mL of (1.0% W/V in 10% HCl) sulphanilamide reagent and was mixed by swirling. After about 3 min, 1.0 mL N-1—naphthalene diamine dihydrochloride (0.1% W/V) was added and mixed thoroughly2. The mark was made-up with distilled water. The blank solution was also subjected to the same treatment as the sample. After about 10–20 min, the absorbance of both the water sample and the blank solutions were measured with a UV–visible spectrophotometer at a wavelength of 543 nm.
Analysis of total petroleum hydrocarbons (TPH)
HP-5 capillary column coated with 5% phenyl methyl siloxane (30 m length × 0.32 mm diameter × 0.25 µm film thickness) (Agilent Technologies) was used as a stationary phase of separation of hydrocarbons from water samples. 1µL of the samples was injected in splitless mode at an injection temperature of 300 °C, and pressure of 13.74psi and a total flow of 21.364 mL/min. Purge flow to split vent was set at 15 mL/min at 0.75 min. The oven was initially programmed at 40 °C (1 min) then ramped at 12 °C/min to 300 °C for 10 min. The temperature of the flame ionization detector was regulated to 300 °C using hydrogen gas. Airflow was at 30 mL/min while nitrogen was used as makeup gas at a flow of 22 mL/min. Agilent 7890B gas chromatography coupled to flame ionization detector (GC-FID) was used for the determination of TPH at 254 nm. After calibration, water samples were analyzed and corresponding TPH concentrations were obtained3,10.
Analysis of oil and grease (OG)
One (1) L separating funnels with retort stand, 100 mL volumetric flask, glass jar, xylene, and anhydrous sodium sulfate were used in determining the concentrations of oil and grease (OG) in the water.
Extraction
Twenty (20) mL xylene was put in a glass jar containing a water sample. The content of the jar was shaken, poured into the separating funnel and shaken again. It was allowed for phase separation and the bottom layer xylene was drained into a 100 mL volumetric flask through a funnel with a plug of glass wool and about 2/3 full with anhydrous Na2SO4.
Another 20 mL xylene was added to the content in the separating funnel, agitated thoroughly and xylene layer was again drained into the same flask. Water was drained into a measuring cylinder and the volume was noted. Separating funnel was rinsed with 20 mL xylene into the same flask as done earlier. It was made up to mark of the extract in the 100 mL volumetric flask with pure xylene.
The oil and grease (OG) was calculated thus:
The concentration of oil reported as OG (mg/L)

$$= frac{{{text{Conc}}. , left( {{text{mg}}/{text{L}},{text{extract}}} right) times {text{DF}} times {text{EV}},{text{(mL)}}}}{{{text{The}},{text{volume}},{text{of}},{text{water}},,{text{(mL)}}}}$$
(3)

In Eq. (3), DF = Dilution factor
CF = Conversion factor from absorbance to mg/L extract
EV = Extraction volume of solvent in (mL).
Analysis of trace metals
Ten (10) mL of water sample was put in a beaker and 2 mL concentrated nitric acid was added to the sample. The mixture was then heated to evaporation and allowed to cool afterward and then transferred into a volumetric flask. It was then allowed to stand for 24 h, after when it was centrifuged at 3000 rpm until clear. The sample was screened for suspended solids which were filtered off before further analysis. The trace metals in the mixture were then read using an atomic absorption spectrophotometer (AAS, Philips model PU 9100) at a wavelength range of 250–350 V using the ABS knob10.
The experimental procedures were conducted as described by Estefan et al.30 and modified by Jones Jr.31.
Quality control and quality assurance
Validation of trace metals
The precision of the AAS was validated by repeating every experimental procedure 3 times. Certified reference materials (CRM) and standard reference materials (SRM) published by the Federal Environmental Protection Agency32 were employed as a guide. The recovery rates ranged from 87 to 95%. The calculated relative standard deviation (SD) was  More

1.
Zhang, Y. et al. Variation of soil microbial community along elevation in the Shennongjia Mountain. For. Sci. 50, 161–166 (2014).
ADS  CAS  Google Scholar 
2.
Davidson, E. A. & Janssens, I. A. Temperature sensitivity of soil carbon decomposition and feedbacks to climate change. Nature 440, 165–173 (2006).
ADS  CAS  PubMed  Article  Google Scholar 

3.
Mariotte, P. et al. Plant-soil feedback: bridging natural and agricultural sciences. Trends Ecol. Evol. 33, 129–136 (2018).
PubMed  Article  Google Scholar 

4.
Mommer, L. et al. Lost in diversity: The interactions between soil-borne fungi, biodiversity and plant productivity. New Phytol. 218, 542–553 (2018).
PubMed  PubMed Central  Article  Google Scholar 

5.
Seneviratne, S. I. et al. Land radiative management as contributor to regional-scale climate adaptation and mitigation. Nat. Geosci. 11, 88–96 (2018).
ADS  CAS  Article  Google Scholar 

6.
Song, X. P. et al. Global land change from 1982 to 2016. Nature 560, 639 (2018).
ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

7.
Tripathi, B. M. et al. Soil pH mediates the balance between stochastic and deterministic assembly of bacteria. ISME J. 12, 1072–1083 (2018).
CAS  PubMed  PubMed Central  Article  Google Scholar 

8.
Zhang, Y. G. et al. Soil bacterial endemism and potential functional redundancy in natural broadleaf forest along a latitudinal gradient. Sci. Rep. 6, 28–39 (2016).
ADS  Article  CAS  Google Scholar 

9.
Yang, Y. F. et al. The microbial gene diversity along an elevation gradient of the Tibetan grassland. ISME J. 8, 430–440 (2014).
CAS  PubMed  Article  Google Scholar 

10.
Zhang, Y. G. et al. Soil bacterial diversity patterns and drivers along an elevational gradient on Shennongjia Mountain, China. Microb. Biotechnol. 8, 739–746 (2015).
PubMed  PubMed Central  Article  Google Scholar 

11.
Shen, C. C. et al. Soil pH drives the spatial distribution of bacterial communities along elevation on Changbai Mountain. Soil Biol. Biochem. 57, 204–211 (2013).
CAS  Article  Google Scholar 

12.
Zhang, Y. G. et al. The microbially mediated soil organic carbon loss under degenerative succession in an alpine meadow. Mol. Ecol. 26, 3676–3686 (2017).
CAS  PubMed  Article  Google Scholar 

13.
Davidson, E. A., Janssens, I. A. & Luo, Y. Q. On the variability of respiration in terrestrial ecosystems: moving beyond Q10. Glob. Change Biol. 12, 154–164 (2006).
ADS  Article  Google Scholar 

14.
Walther, G. R. et al. Ecological responses to recent climate change. Nature 416, 389–395 (2002).
ADS  CAS  PubMed  Article  Google Scholar 

15.
Zhang, Q. et al. Effects of clipping frequency on the relationships between species diversity and productivity in temperate steppe. Int. J. Agric. Biol. 20, 2325–2328 (2018).
Google Scholar 

16.
Guo, Z. G., Cheng, G. D. & Wang, G. X. Plant Diversity of Alpine Kobresia Meadow in the Northern Region of the Tibetan Plateau. J. Glaciol. Geocryol. 26, 95–100 (2004).
Google Scholar 

17.
Fu, B. et al. Hydrogeomorphic ecosystem responses to natural and anthropogenic changes in the Loess Plateau of China. Annu. Rev. Earth Planet. Sci. 45, 223–243 (2017).
ADS  CAS  Article  Google Scholar 

18.
Yang, Y., Dou, Y. & An, S. Testing association between soil bacterial diversity and soil carbon storage on the loess plateau. Sci. Total Environ. 626, 48–58 (2018).
ADS  CAS  PubMed  Article  Google Scholar 

19.
Yang, Y., Cheng, H., Liu, L. X., Dou, Y. X. & An, S. S. Comparison of soil microbial community between planted woodland and natural grass vegetation on the Loess Plateau. For. Ecol. Manage. 460, 117–128 (2020).
Article  Google Scholar 

20.
Tong, X. W. et al. Rasmus Fensholt. Increased vegetation growth and carbon stock in China karst via ecological engineering. Nat. Sustain. 1, 44–50 (2018).
Article  Google Scholar 

21.
Liu, Y. et al. Temporal and spatial succession and dynamics of soil fungal communities in restored grassland on the Loess Plateau in China. Land Degrad. Dev. 30, 1273–1287 (2019).
Article  Google Scholar 

22.
Bardgett, R. D. et al. Below-ground microbial community development in a high temperature world. Oikos 85, 193–203 (1999).
Article  Google Scholar 

23.
Nave, L. E. et al. Reforestation can sequester two petagrams of carbon in us top soils in a century. Proc. Natl. Acad. Sci. USA https://doi.org/10.1073/pnas.1719685115 (2018).
Article  PubMed  Google Scholar 

24.
Chen, J. S., Zhu, R. F. & Zhang, Y. X. The effect of nitrogen addition on seed yield and yield components of Leymus chinensis. J. Soil Sci. Plant Nutr. 13, 329–339 (2013).
Google Scholar 

25.
Lange, M. et al. Plant diversity increases soil microbial activity and soil carbon storage. Nat. Commun. 6, 6707 (2015).
ADS  CAS  PubMed  Article  Google Scholar 

26.
Lal, R. Digging deeper: a holistic perspective of factors affecting soil organic carbon sequestration in agroecosystems. Glob. Change Biol. 24, 3285–3301 (2018).
ADS  Article  Google Scholar 

27.
Gao, Q. Z. et al. Alpine grassland degradation index and its response to recent climate variability in Northern Tibet, China. Quatern. Int. 226, 143–150 (2010).
Article  Google Scholar 

28.
Li, N. et al. Short-term effects of temperature enhancement on community structure and biomass of alpine meadow in the Qinghai-Tibet Plateau. Acta Ecol. Sin. 31, 0895–0905 (2011).
CAS  Google Scholar 

29.
Gao, Y. H. et al. Vegetation net primary productivity and its response to climate change during 2001–2008 in the Tibetan Plateau. Sci. Total Environ. 444, 356–362 (2013).
ADS  CAS  PubMed  Article  Google Scholar 

30.
Chen, J. S. et al. Effects of clipping and fertilizing n on the relationship between diversity and productivity of Leymus Chinensis Meadow. Acta Agrestia Sin. 24, 910–914 (2016).
Google Scholar 

31.
Peng, F. et al. Intensified plant N and C pool with more available nitrogen under experimental warming in an alpine meadow ecosystem. Ecol. Evol. 6, 8546–8555 (2016).
PubMed  PubMed Central  Article  Google Scholar 

32.
Li, L. et al. Precipitation overrides warming in mediating soil nitrogen pools in an alpine grassland ecosystem on the Tibetan Plateau. Sci. Rep. 6, 31–38 (2016).
Google Scholar 

33.
Niu, S. Q. et al. Microbial diversity in saline alkali soil from Hexi Corridor analyzed by Illumina MiSeq high-throughput sequencing system. Microbiology China 9, 66–72 (2017).
Google Scholar 

34.
Chaffron, S. et al. A global network of coexisting microbes from environmental and whole-genome sequence data. Genome Res. 20, 947–959 (2010).
CAS  PubMed  PubMed Central  Article  Google Scholar 

35.
Steele, J. A. et al. Marine bacterial, archaeal and protest an association networks reveal ecological linkages. ISME J. 5, 1414–1425 (2011).
PubMed  PubMed Central  Article  Google Scholar 

36.
Kreimer, A. et al. NetCmpt: a network-based tool for calculating the metabolic competition between bacterial species. Bioinformatics 28, 2195–2197 (2012).
CAS  PubMed  Article  Google Scholar 

37.
O’Brien, J. D. et al. A Bayesian approach to inferring the phylogenetic structure of communities from metagenomic data. Genetics 197, 925–937 (2014).
PubMed  PubMed Central  Article  Google Scholar 

38.
Deng, Y. et al. Elevated carbon dioxide accelerates the spatial turnover of soil microbial communities. Glob. Change Biol. 22, 957–964 (2016).
ADS  Article  Google Scholar 

39.
Zhou, J. Z. et al. Functional molecular ecological networks. mBio 1, e00169-10 (2010).
PubMed  PubMed Central  Article  CAS  Google Scholar 

40.
Zhou, J. Z. et al. Phylogenetic molecular ecological network of soil microbial communities in response to elevated CO2. mBio 2, e00122-e211 (2011).
PubMed  PubMed Central  Article  Google Scholar 

41.
Deng, Y. et al. Molecular ecological network analyses. BMC Bioinformatics 13, 113 (2012).
PubMed  PubMed Central  Article  Google Scholar 

42.
Wang, J. F. & Wu, Q. B. Influences of the vegetaion degradation on the shallow cryic soil environment in the wet meadow areas on the Qinghai-Tibetan plateau. J. Lanzhou Univ. 47, 39–45 (2011).
Google Scholar 

43.
Wang, Y. et al. Artificial reforestation produces less diverse soil nitrogen-cycling genes than natural restoration. Ecosphere 10, e02562 (2019).
Google Scholar 

44.
Bao, S. D. Soil and Agricultural Chemical Analysis (China Agriculture Press, Beijing, 2000).
Google Scholar 

45.
Zhou, J. Z., Bruns, M. A. & Tiedje, J. M. DNA recovery from soils of diverse composition. Appl. Environ. Microbiol. 62, 316–322 (1996).
CAS  PubMed  PubMed Central  Article  Google Scholar 

46.
Shen, C. C. et al. Dramatic increases of soil microbial functional gene diversity at the treeline ecotone of Changbai Mountain. Front. Microbiol. 7, 1184 (2016).
PubMed  PubMed Central  Google Scholar 

47.
Chen, W., Koide, R. T. & Eissenstat, D. M. Nutrient foraging by mycorrhizas: From species functional traits to ecosystem processes. Funct. Ecol. 32, 858–869 (2018).
Article  Google Scholar 

48.
Guo, X. et al. Climate warming leads to divergent succession of grassland microbial communities. Nat. Glob. Change 8, 813–818 (2018).
ADS  Article  Google Scholar 

49.
Jiao, S., Xu, Y., Zhang, J. & Lu, Y. Soil bacterial community dynamics reflect changes in plant community and soil properties during the secondary succession of abandoned farmland in the Loess Plateau. Microbiome 6, 146–156 (2018).
PubMed  PubMed Central  Article  Google Scholar 

50.
Brabcová, V., Štursová, M. & Baldrian, P. Nutrient content affects the turnover of fungal biomass in forest topsoil and the composition of associated microbial communities. Soil Biol. Biochem. 118, 187–198 (2018).
Article  CAS  Google Scholar 

51.
He, J. Z. & Ge, Y. Recent advances in soil microbial biogeography. Acta Ecol. Sin. 28, 5571–5582 (2008).
CAS  Google Scholar 

52.
Horner-Devine, M. C., Lage, M., Hughes, J. B. & Bohannan, B. J. M. A taxa-area relationship for bacteria. Nature 432, 750–753 (2004).
ADS  CAS  PubMed  Article  Google Scholar 

53.
McCann, K., Hastings, A. & Huxel, G. R. Weak trophic interactions and the balance of nature. Nature 395, 794–798 (1998).
ADS  CAS  Article  Google Scholar 

54.
Pimm, S. L., Lawton, J. H. & Cohen, J. E. Food web patterns and their consequences. Nature 350, 669–674 (1991).
ADS  Article  Google Scholar 

55.
Fuhrman, J. A. Microbial community structure and its functional implications. Nature 459, 193–199 (2009).
ADS  CAS  PubMed  Article  Google Scholar 

56.
Sun, X., Gao, Y. & Yang, Y. F. Recent advancement in environmental research with metagenomics tools. Biodivers. Sci. 21, 393–400 (2013).
CAS  Google Scholar 

57.
Barabási, A. L. & Oltvai, Z. N. Network biology: understan Mommer the cell’s functional organization. Nat. Rev. Genet. 5, 101–113 (2004).
PubMed  Article  CAS  Google Scholar 

58.
Ding, J. J. et al. Integrated metagenomics and network analysis of soil microbial community of the forest timberline. Sci. Rep. 5, 7994 (2015).
CAS  PubMed  PubMed Central  Article  Google Scholar 

59.
Brundrett, M. C. & Tedersoo, L. Evolutionary history of mycorrhizal symbioses and global host plant diversity. New Phytol. 220, 55–61 (2018).
Article  Google Scholar 

60.
Xiong, J. B. et al. Geographic distance and pH drive bacterial distribution in alkaline lake sediments across Tibetan Plateau. Environ. Microbiol. 14, 2457–2466 (2012).
CAS  PubMed  PubMed Central  Article  Google Scholar  More

1.
Ege, V. A revision of the Genus Anguilla Shaw. Dana Rep. 16, 8–256 (1939).
Google Scholar 
2.
Arai, T. Taxonomy and distribution. In Biology and Ecology of Anguillid Eels (ed. Arai, T.) 1–20 (CRC Press, Boca Raton, 2016).
Google Scholar 

3.
Arai, T. & Chino, N. Diverse migration strategy between freshwater and seawater habitats in the freshwater eels genus Anguilla. J. Fish Biol. 81, 442–455 (2012).
CAS  PubMed  Article  PubMed Central  Google Scholar 

4.
Arai, T., Limbong, D., Otake, T. & Tsukamoto, K. Recruitment mechanisms of tropical eels, Anguilla spp., and implications for the evolution of oceanic migration in the genus Anguilla. Mar. Ecol. Prog. Ser. 216, 253–264 (2001).
ADS  CAS  Article  Google Scholar 

5.
Shiao, J. C., Tzeng, W. N., Collins, A. & Iizuka, Y. Role of marine larval duration and growth rate of glass eels in determining the distribution of Anguilla reinhardtii and A. australis on Australian eastern coasts. Mar. Freshw. Res. 53, 1–10 (2002).
Article  Google Scholar 

6.
Arai, T., Abdul Kadir, S. R. & Chino, N. Year-round spawning by a tropical catadromous eel Anguilla bicolor bicolor. BMar. Biol. 163, 37 (2016).
Google Scholar 

7.
Arai, T. & Abdul Kadir, S. R. Opportunistic spawning of tropical anguillid eels Anguilla bicolor bicolor and A. bengalensis bengalensis. Sci. Rep. 7, 41649 (2017).
ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

8.
Beumer, J. & Sloane, R. Distribution and abundance of glass-eels Anguilla spp. in east Australian waters. Int. Rev. Hydrobiol. 75, 721–736 (1990).
Article  Google Scholar 

9.
Arai, T., Aoyama, J., Limbong, D. & Tsukamoto, K. Species composition and inshore migration of the tropical eels Anguilla spp. recruiting to the estuary of the Poigar River, Sulawesi Island. Mar. Ecol. Prog. Ser. 188, 299–303 (1999).
ADS  Article  Google Scholar 

10.
Sugeha, H., Arai, T., Miller, M. J., Limbong, D. & Tsukamoto, K. Inshore migration of the tropical eels Anguilla spp. recruiting to the Poigar River estuary on North Sulawesi Island. Mar. Ecol. Prog. Ser. 221, 233–243 (2001).
ADS  Article  Google Scholar 

11.
Hewavitharane, C.A., Pickering, T.D., Ciro, R., Mochioka. N. Species composition, abundance and seasonal recruitment patterns of freshwater eels (Anguilla spp.) to Viti Levu, Fiji Islands, in the western South Pacific. Mar. Freshw. Res. 69,1704–1711 (2018).
Article  Google Scholar 

12.
Leander, N. J., Shen, K. N., Chen, R. T. & Tzeng, W. N. Species composition and seasonal occurrence of recruiting glass eels (Anguilla spp.) in the Hsiukuluan Rive, River Eastern Taiwan. Zool. Stud. 51, 59–71 (2012).
CAS  Google Scholar 

13.
Arai, T., Chino, N., Zulkifli, S. Z. & Ismail, A. Notes on the occurrence of the tropical eel Anguilla bicolor bicolor in Peninsular Malaysia, Malaysia. J. Fish Biol. 80, 692–697 (2012).
CAS  PubMed  Article  Google Scholar 

14.
Arai, T. First record of a tropical mottled eel, Anguilla bengalensis bengalensis (Actinopterygii: Anguillidae) from the Langkawi Islands, Peninsular Malaysia, Malaysia. Mar. Biodivers. Rec. 7, e38 (2014).
Article  Google Scholar 

15.
Arai, T., Chin, T. C., Kwong, K. O. & Siti Azizah, M. N. Occurrence of the tropical eels, Anguilla bengalensis bengalensis and A. bicolor bicolor in Peninsular Malaysia, Malaysia and implications for the eel taxonomy. Mar. Biodivers. Rec. 8, e28 (2015).
Article  Google Scholar 

16.
Arai, T. & Wong, L. L. Validation of the occurrence of the tropical eels, Anguilla bengalensis bengalensis and A. bicolor bicolor at Langkawi Island in Peninsular Malaysia, Malaysia. Tropic. Ecol. 57, 23–31 (2016).
Google Scholar 

17.
Abdul Kadir, S. R., Abdul Rasid, M. H. F., Kwong, K. O., Wong, L. & Arai, T. Occurrence and the ecological implication of a tropical anguillid eel Anguilla marmorata from peninsular Malaysia. Zookeys 695, 103–110 (2017).
Article  Google Scholar 

18.
Arai, T. & Abdul Kadir, S. R. Diversity, distribution and different habitat use among the tropical freshwater eels of genus Anguilla. Sci. Rep. 7, 7593 (2017).
ADS  PubMed  PubMed Central  Article  CAS  Google Scholar 

19.
Arai, T. Evidence of local short-distance spawning migration of tropical freshwater eels, and implications for the evolution of freshwater eel migration. Ecol. Evol. 4, 3812–3819 (2014).
PubMed  PubMed Central  Article  Google Scholar 

20.
Arai, T., Kotake, A. & McCarthy, T. K. Habitat use by the European eel Anguilla auguilla in Irish waters. Estuar. Coast. Shelf Sci. 67, 569–578 (2006).
ADS  Article  Google Scholar 

21.
Shiao, J. C., Ložys, L., Iizuka, Y. & Tzeng, W. N. Migratory patterns and contribution of stocking to the population of European eel in Lithuanian waters as indicated by otolith Sr:Ca ratios. J. Fish Biol. 69, 749–769 (2006).
CAS  Article  Google Scholar 

22.
Lin, Y. J., Yalçin-Özdilek, S., Iizuka, Y., Gümüş, A. & Tzeng, W. N. Migratory life history of European eel Anguilla anguilla from freshwater regions of the River Asi, southern Turkey and their high otolith Sr:Ca ratios. J. Fish Biol. 78, 860–868 (2011).
PubMed  Article  Google Scholar 

23.
Arai, T., Kotake, A., Harrod, C., Morrissey, M. & Mccarthy, T. K. Ecological plasticity of the European eel Anguilla anguilla in a tidal Atlantic lake system in Ireland. J. Mar. Biol. Assoc. UK 99, 1189–1195 (2019).
CAS  Article  Google Scholar 

24.
Thibault, I. et al. Facultative catadromy in American eels: testing the conditional strategy hypothesis. Mar. Ecol. Prog. Ser. 344, 219–229 (2007).
ADS  Article  Google Scholar 

25.
Jessop, B. M. et al. Otolith Sr:Ca and Ba:Ca may give inconsistent indications of estuarine habitat use for American eels (Anguilla rostrata). Environ. Biol. Fish. 93, 193–207 (2012).
Article  Google Scholar 

26.
Tsukamoto, K. & Arai, T. Facultative catadromy of the eel, Anguilla japonica, between freshwater and seawater habitats. Mar. Ecol. Prog. Ser. 220, 365–376 (2001).
Article  Google Scholar 

27.
Tzeng, W. N., Shiao, J. C. & Iizuka, Y. Use of otolith Sr:Ca ratios to study the riverine migratory behaviors of the Japanese eel Anguilla japonica. Mar. Ecol. Prog. Ser. 245, 213–221 (2002).
ADS  Article  Google Scholar 

28.
Shiao, J. C., Iizuka, Y., Chang, C. W. & Tzeng, W. N. Disparities in habitat use and migratory behaviour between tropical eel Anguilla marmorata and temperate eel A. japonica in four Taiwanese rivers. Mar. Ecol. Prog. Ser. 261, 233–242 (2003).
ADS  Article  Google Scholar 

29.
Chino, N. & Arai, T. Relative contribution of migratory type on the reproduction of migrating silver eels, Anguilla japonica, collected off Shikoku Island, Japan. Mar. Biol. 156, 661–668 (2009).
CAS  Article  Google Scholar 

30.
Arai, T., Kotake, A., Lokman, P. M., Miller, M. J. & Tsukamoto, K. Evidence of different habitat use by New Zealand freshwater eels, Anguilla australis and A. dieffenbachii, as revealed by otolith microchemistry. Mar. Ecol. Prog. Ser. 266, 213–225 (2004).
ADS  Article  Google Scholar 

31.
Briones, A. A., Yambot, A. V., Shiao, J. C., Iizuka, Y. & Tzeng, W. N. Migratory pattern and habitat use of tropical eels Anguilla spp. (teleostei: Anguilliformes: Anguillidae) in the Philippines, as revealed by otolith microchemistry. Raffl. Bull. Zool. 14, 141–149 (2007).
Google Scholar 

32.
Chino, N. & Arai, T. Migratory history of the giant mottled eel (Anguilla marmorata) in the Bonin Islands of Japan. Ecol. Freshw. Fish. 19, 19–25 (2010).
Article  Google Scholar 

33.
Lin, Y. J. et al. Regional variation in otolith Sr:Ca ratios of African longfinned eel Anguilla mossambica and mottled eel Anguilla marmorata: a challenge to the classic tool for reconstructing migratory histories of fishes. J. Fish Biol. 81, 427–441 (2012).
PubMed  Article  PubMed Central  Google Scholar 

34.
Arai, T., Chino, N. & Le, D. Q. Migration and habitat use of the tropical eels Anguilla marmorata and A. bicolor pacifica in Vietnam. Aquat. Ecol. 47, 57–65 (2013).
Article  Google Scholar 

35.
Arai, T. & Chino, N. Opportunistic migration and habitat use of the giant mottled eel Anguilla marmorata (Teleostei: Elopomorpha). Sci. Rep. 8, 5666 (2018).
ADS  PubMed  PubMed Central  Article  CAS  Google Scholar 

36.
Arai, T. & Chino, N. Variations in the migratory history of the tropical catadromous eels Anguilla bicolor bicolor and A. bicolor pacifica in south-east Asian waters. J. Fish Biol. 94, 752–758 (2019).
CAS  PubMed  Article  PubMed Central  Google Scholar 

37.
Chino, N. & Arai, T. Occurrence of marine resident tropical eel Anguilla bicolor bicolor in Indonesia. Mar. Biol. 157, 1075–1081 (2010).
Article  Google Scholar 

38.
Chino, N. & Arai, T. Habitat use and habitat transitions in the tropical eel Anguilla bicolor bicolor. Environ. Biol. Fish. 89, 571–578 (2010).
Article  Google Scholar 

39.
Tzeng, W. N. Effects of salinity and ontogenetic movements on strontium:calcium ratios in the otoliths of the Japanese Eel, Anguilla japonica Temminck and Schlegel. J. Exp. Mar. Biol. Ecol. 199, 111–122 (1996).
CAS  Article  Google Scholar 

40.
Lin, S. H., Chang, C. W., Iizuka, Y. & Tzeng, W. N. Salinities, not diets, affect strontium/calcium ratios in otoliths of Anguilla japonica. J. Exp. Mar. Biol. Ecol. 34, 254–263 (2007).
Article  CAS  Google Scholar 

41.
Arai, T. & Chino, N. Influence of water salinity on the strontium:calcium ratios in otoliths of the giant mottled eel Anguilla marmorata. Environ. Biol. Fish. 100, 281–286 (2017).
Article  Google Scholar 

42.
Jellyman, D. J. & Todd, P. R. Why are migrating male shortfinned eels (Anguilla australis) in Lake Ellesmere, New Zealand, getting smaller but not younger?. Bull. Fr. Peche. Piscic. 349, 141–152 (1998).
Article  Google Scholar 

43.
Jellyman, D. J., Glova, G. J., Sagar, P. M. & Sykes, J. R. E. Spatiotemporal distribution of fish in the Kakanui River estuary, South Island, New Zealand. NZ J. Mar. Freshw. Res. 30, 103–118 (2001).
Google Scholar 

44.
McDowall, R. M. New Zealand Freshwater Fishes: A Natural History and Guide (Heinemann Reed, Auckland, 1990).
Google Scholar 

45.
Jowett, I. G. & Richardson, J. Habitat preferences of common, riverine New Zealand native fishes and implication for flow management. NZ J. Mar. Freshw. Res. 29, 13–23 (1995).
Article  Google Scholar 

46.
Glova, G. J., Jellyman, D. J. & Bonnett, M. L. Factors associated with the distribution and habitat of eels (Anguilla spp.) in three New Zealand lowland streams. NZ J. Mar. Freshw. Res. 32, 255–269 (1998).
Article  Google Scholar 

47.
Bozeman, E. K., Helfman, G. S. & Richardson, T. Population size and home range of American eels in a Georgia tidal creek. Trans. Am. Fish. Soc. 114, 821–825 (1985).
Article  Google Scholar 

48.
Ford, T. E. & Mercer, E. Density, size distribution and home range of American eels, Anguilla rostrata, in a Massachusetts salt marsh. Environ. Biol. Fish. 17, 309–314 (1986).
Article  Google Scholar 

49.
Chisnall, B. L. & Kalish, J. M. Age validation and movement of freshwater eels (Anguilla dieffenbachii and A. australis) in a New Zealand pastoral stream. NZ J. Mar. Freshw. Res. 27, 333–338 (1993).
Article  Google Scholar 

50.
Oliveira, K. Movement and growth rates of yellow phase American eels in the Annaquatucket River, Rhode Island. Trans. Am. Fish. Soc. 126, 638–646 (1997).
Article  Google Scholar 

51.
Jellyman, D. J. & Sykes, J. R. E. Diel and seasonal movements of radio-tagged freshwater eels, Anguilla spp., in two New Zealand streams. Environ. Biol. Fish. 66, 143–154 (2003).
Article  Google Scholar 

52.
Poole, W. R. & Reynolds, J. D. Variability in growth rate in European eel Anguilla anguilla (L.) in a Western Irish catchment. Proc. R. Irish Acad. 98B, 141–145 (1998).
Google Scholar 

53.
Jessop, B. M. Geographic effects on American eel (Anguilla rostrata) life history characteristics and strategies. Can. J. Fish. Aquat. Sci. 2010(67), 326–346 (2010).
Article  Google Scholar 

54.
Kotake, A. et al. Ecological aspects of Japanese eels, Anguilla japonica, collected from coastal areas of Japan. Zool. Sci. 24, 1213–1221 (2007).
PubMed  Article  Google Scholar 

55.
Jellyman, D. J. Status of New Zealand fresh-water eel stocks and management initiatives. ICES J. Mar. Sci. 64, 1379–1386 (2007).
Article  Google Scholar 

56.
Chisnall, B. L. & Hicks, B. J. Age and growth of longfinned eels (Anguilla dieffenbachii) in pastoral and forested streams in the Waikato River basin, and in two hydro-electric lakes in the North Island, New Zealand. NZ J. Mar. Freshw. Res. 27, 317–332 (1993).
Article  Google Scholar 

57.
Iizuka, Y. Electron microprobe study of otolith: Migratory behavior and habitat of three major temperate species of eels. JEOL News 47, 33–50 (2012).
Google Scholar 

58.
Arai, T., Otake, T. & Tsukamoto, K. Drastic changes in otolith microstructure and microchemistry accompanying the onset of metamorphosis in the Japanese eel Anguilla japonica. Mar. Ecol. Prog. Ser. 161, 17–22 (1997).
ADS  CAS  Article  Google Scholar 

59.
Oliveira, K. Field validation of annular growth rings in the American eel Anguilla rostrata, using tetracycline-marked otolith. Fish. Bull. 94, 186–189 (1996).
Google Scholar 

60.
Graynoth, E. Improved otolith preparation, ageing and back-calculation techniques for New Zealand freshwater eels. Fish. Res. 42, 137–146 (1999).
Article  Google Scholar 

61.
Arkhipkin, A. I., Schuchert, P. C. & Danyushevsky, L. Otolith chemistry reveals fine population structure and close affinity to the Pacific and Atlantic oceanic spawning grounds in the migratory southern blue whiting (Micromesistius australis australis). Fish. Res. 96, 188–194 (2009).
Article  Google Scholar 

62.
Schuchert, P. C., Arkhipkin, A. I. & Koenig, A. E. Traveling around Cape Horn: Otolith chemistry reveals a mixed stock of Patagonian hoki with separate Atlantic and Pacific spawning grounds. Fish. Res. 102, 80–86 (2010).
Article  Google Scholar  More

SARS-CoV-2 genome data sets and associated travel history
To focus on the early stage of COVID-19 spread, we analyzed SARS-CoV-2 genome sequences and metadata available in GISAID on March 10th8. We curated a data set of 305 genomes by removing error-prone sequences, keeping only genomes with appropriate metadata, and a single genome from patients with multiple genomes available. We assigned each genome a global lineage designation based on the nomenclature scheme outlined in Rambaut et al.28 using pangolin v1.1.14 (https://github.com/hCoV-2019/pangolin), lineages data release 2020-05-19 (https://github.com/hCoV-2019/lineages). We aligned the remaining genomes using MAFFT v.7.45329 and partially trimmed the 5′ and 3′ ends. All sequences were associated with exact sampling dates in their meta-information, except for one genome from Anhui with known month of sampling. Upon visualizing root-to-tip divergence as a function of sampling time, using TempEst v.1.5.330 based on an ML tree inferred with IQ-TREE v.2.0-rc131, we removed one potential outlier. The root-to-tip plots without the outlier are shown in Supplementary Fig. S3. We formally tested for temporal signal using BETS32. The final 282 genomes were sampled from 28 different countries, with Chinese samples originating from 13 provinces, one municipality (Beijing), and one special administrative area (Hong Kong), which we considered as separate locations in our (discrete) phylogeographic analyses. Phylogenetic signal in the data set was explored through likelihood mapping analysis33 (Supplementary Fig. S4).
We searched for travel history data associated with the genomes in the GISAID records, media reports, and publications and retrieved recent travel locations for 64 genomes (22.5%, Supplementary Table 2): 43 traveled/returned from Hubei (Wuhan), 1 from Beijing, 3 from China without further detail (which we associated with an appropriate ambiguity code in our phylogeographic analysis that represents all sampled Chinese locations), 2 from Singapore, 1 from Southeast Asia (which we also associated with an ambiguity code that represents all sampled Southeast Asian locations), 7 from Italy, and 7 from Iran. In this data set, Italy is better represented by recent travel locations than actual samples (n = 4) and Iran is exclusively represented by travelers returning from this country. For 46 out of the 64 genomes, we retrieved the date of travel, which represents the most recent time point at which the ancestral lineage circulated in the travel location.
In order to examine (i) to what extent our reconstructions could be updated by the genome data that has become available retrospectively for the same locations and the same time period before March 10 ~4 months after this date, and (ii) how sampling bias can be mitigated by downsampling from the larger collection of available genomes, we assembled an additional data set of 500 genomes. For this purpose, SARS-CoV-2 genomes were downloaded from GISAID on June 23, 2020 and processed according to the COG-UK pre-analysis pipeline (https://github.com/COG-UK/grapevine). Briefly, sequences were aligned to the reference sequence Wuhan-Hu-1 (Genbank accession number NC_045512) using Minimap2 v.2.1734. Problematic sites were masked (https://virological.org/t/issues-with-sars-cov-2-sequencing-data/473), and sequences with More

SCALE model
Combining the processes of biogeochemical transformation, soil erosion/deposition and resultant landscape evolution, and bioturbation by soil fauna, the governing equations, which conserves SOC mass, can be expressed as24:

$$begin{aligned} frac{ partial }{ partial t } int _{0}^{Z} mathbf{C } dz =int _{0}^{Z} mathbf{g } dz – nabla cdot mathbf{q }_C + int _{0}^{Z} nabla cdot big [D(z) nabla mathbf{C } big ] dz end{aligned}$$
(1)

where (mathbf{C }) is the SOC concentration ([M L^{-3}]), (mathbf{C } = [ C_l,C_h,C_b]^T) represents the fast (or litter), slow (or humus), and microbial biomass pool, respectively; (mathbf{g }) is the rate of the biogeochemical transformation process, which is a function of the fast (or litter), slow (or humus), and microbial biomass pool; (nabla cdot mathbf{q }_C) is the surface SOC flux associated with soil transport; and (D(z)) is the bioturbation diffusivity. When the vertical column is discretized into several layers, the equations can be written as:

$$begin{aligned} text {Top soil layer: }&frac{ partial big ( mathbf{C }_1 z_1 big ) }{ partial t } = mathbf{g }_1 z_1 – nabla cdot mathbf{q }_C end{aligned}$$
(2a)

$$begin{aligned} text {Sublayers: }&frac{ partial mathbf{C }_n }{ partial t } = mathbf{g }_n + nabla cdot big [D(z) nabla mathbf{C } big ] dz end{aligned}$$
(2b)

where the subscripts 1 and n denote the surface soil layer and the (n{mathrm{^{th}}}) layer below-surface, respectively. Other details of these equations can be found in24.
The mechanisms of soil transport and the resultant landscape evolution can be categorized into two groups: overland flow-driven transport and diffusion-driven transport from other disturbances (e.g., wind, animal, and raindrop splash). The 2-D mass conservation equation of soil transport and the resultant landscape evolution follows Exner equation:

$$begin{aligned} frac{partial eta }{partial t} = U -nabla cdot q_d – nabla cdot q_s end{aligned}$$
(3)

where (eta) is soil surface elevation [L]; U is the rate of tectonic uplift or glacial rebound ([L T^{-1}]); (q_d) is the volume flux of sediment per unit width by hillslope diffusion ([L^2 T^{-1}]); (q_s) is the volume flux of sediment per unit width by overland flow ([L^2 T^{-1}]).
The diffusion-driven transport (nabla cdot q_d) is a combination of wind erosion, animal disturbance, soil creep, raindrop splash, and biogenic transport. The 2-D equation of (q_d) is expressed as a linear relationship with slope47:

$$begin{aligned} q_d = – D_{x} frac{ partial eta }{partial x} – D_{y} frac{ partial eta }{partial y} end{aligned}$$
(4)

where (D_x) and (D_y) are the soil diffusion coefficient in x and y direction, respectively ([L^2 T^{-1}]). The overland flow-driven transport (nabla cdot q_s) is a combined form of the divergence of stream-power but limited by the detachment capacity48:

$$begin{aligned} nabla cdot q_s = min Bigg ( D_c, frac{q_{s,out} – sum q_{s, in}}{d_s} Bigg ) end{aligned}$$
(5)

where (D_c) is the detachment capacity, which is the upper limit of of local erosion rate ([L/S]); (q_{s,out}) is the sediment flux out of a cell and (sum q_{s, in}) is the total sediment flux into a cell assumed at sediment transport capacity.
The rate of change of SOC on the surface driven by soil transport is (nabla cdot q_s), which has a linear relationship with soil transport flux:

$$begin{aligned} nabla cdot mathbf{q }_{C} = nabla cdot ( k_{soc} mathbf{C }_{1} q_d) + nabla cdot ( k_{soc} mathbf{C }_{1} q_s) end{aligned}$$
(6)

where (mathbf{C } = [C_l , C_h , C_b ]^T), and the subscript 1 denotes the surface soil layer; (q_d) and (q_s) are soil transport flux of diffusion and overland flow; (k_{soc}) is an enrichment ratio, which represents a preferential transport (mobilization and deposition) of SOC. The SOC fluxes driven by diffusion and overland flow sediment transport are:

$$begin{aligned} nabla cdot ( k_{soc} mathbf{C }_{1} q_d)= & {} – frac{partial }{partial x} bigg (k_{soc} mathbf{C }_{1} D_{x} frac{ partial eta }{partial x} bigg ) – frac{partial }{partial y} bigg (k_{soc} mathbf{C }_{1} D_{y} frac{ partial eta }{partial y}bigg ) end{aligned}$$
(7)

$$begin{aligned} nabla cdot big (k_{soc} mathbf{C }_{1} q_s big )= & {} left{ begin{array}{ll} k_{soc} mathbf{C }_{1} D_c , &{}quad mathrm{if} D_c < frac{q_{s,out} - sum q_{s, in}}{d_s} \ frac{ k_{soc} mathbf{C }_{1,out} q_{s,out} - sum big ( k_{soc}mathbf{C }_{1,in} q_{s, in} big ) }{d_s} , &{}quad mathrm{otherwise} end{array} right. end{aligned}$$ (8) The initial conditions include land surface elevation, SOC profiles, soil moisture, and surface water depth for each grid box. Table S1 (in Supplementary Information) lists variables associated with initial conditions. We simulate the SOC surface transport and vertical transformation at a daily time step with 2 (times) 2 (mathrm{m}^2) spatial resolution on the surface and a range of vertical grid sizes varying from 5 to 60 (text{cm}). The initial elevation is from lidar (Light Detection and Ranging) DEM (digital elevation model) with a 2-m resolution, collected by a drone in 2015. The SOC profiles at each 2-D grid are estimated by combining the surface SOC contents from the field survey in 2016 and SOC profiles from soil cores sampled at a nearby site 2-km away (see Section below). We assume that the initial soil thickness is 1-m because the most active soil thickness49 regarding SOC and soil moisture are within the top 1 (text{m}), and the deeper SOC contents are quite uniform19,20, indicating a relatively stable condition. Also, we neglect bedrock weathering, and therefore, the soil thickness change is only from the surface soil erosion or deposition. The SOC stock and soil thickness change in the results are compared with the initial SOC stock in the top 1 m. Other initial values such as soil moisture profile and surface water depth (i.e., initial values assigned spatially uniform at each grid point) have a short-time memory in that the impacts only last for a few days to weeks and are predominantly determined by the external forcing and resulting dynamics (more information in Supplementary Information Table S1). Initial soil organic carbon profiles SOC in a natural setting exponentially decreases with soil depth. Here we assume the following relationship between SOC content and soil depth: $$begin{aligned} SOC(Z) = a e^{-bZ} + c end{aligned}$$ (9) where (text{SOC}(text{Z})) represents SOC content at depth (text{Z}); (text{Z}) is zero at the surface, and positive downward; (a), (b), and (c) are positive and constant parameters, where (a+c) represents surface SOC content, (b) represents the decay rate, and (c) represents the relative stable or immobile SOC at depth. Parameters (b) and (c) in this study are estimated from twenty deep cores—ten sites at hills of the Reference Watershed, eight sites at GLC Watershed in the hills, and two sites in the consolidated gully land (Figure S7a). The SOC contents are shown as dots in Figure S6 (in Supplementary Information). We assume the exponential decay rate ((b)) and the immobile SOC ((c)) are spatially uniform in the cornfield in the consolidated gully area and natural field, respectively. We use the least square non-linear method to fit the sample points, and the fitted curves are solid lines in Figure S6 (in Supplementary Information) with the corresponding relationships obtained as: $$begin{aligned} text {For trees/shrubs: }&SOC(Z) = 8.04 e^{-9.63Z}+2.55 end{aligned}$$ (10a) $$begin{aligned} text {For crops: }&SOC(Z) = 3.71 e^{-7.06Z}+2.27 end{aligned}$$ (10b) hence in the natural area, (b= 9.63) and (c = 2.55); and in the cornfield, (b = 7.06) and (c = 2.27). The parameter (a) varies spatially. Soil samples near the soil surface were collected across the whole area of both the GLC and Reference Watershed in 2016. To ensure that the sampling sites are uniformly distributed and represent all land cover types in each watershed, we superimposed an 80-m (times) 80-m grid on the DEM map. In each grid cell, we selected one representative site to collect soil samples; in consolidated gullies, adjacent sampling sites were spaced at an interval of  40-m because of the relatively narrow width. SOC was measured in a laboratory by using the dichromate oxidation method50. In the GLC Watershed, soil samples were collected at 89 locations with 178 samples (0–10 cm and 10–20 cm); in the Reference Watershed, soil samples were collected at 72 locations with 144 samples (Figure 2a, in Supplementary Information). We used Kriging51 to interpolate the SOC content in each 2-D grid box in the two watersheds. Then (a) is back-calculated with the given values of the two layers at 5 cm and 15 cm (the middle point of the two layers, respectively). The final surface SOC (which equals (a+c)) is shown in Figure S7a (in Supplementary Information). Forcing data To explore the co-evolution of the landscape and the vertical profiles of SOC over the decadal time scale, we target a 50-years simulation. The meteorological data is collected from the China National Field Observation Station in An’sai ((36^{circ } 51^{prime } 30^{prime prime } N), (109^{circ } 19^{prime } 23^{prime } E); data source: http://asa.cern.ac.cn), 44 km away to the northwest from the two watersheds, with a 10 years record from 2008 to 2017, which is the best available data for the simulation. The mean annual precipitation is 560 mm. These data are used to train a stochastic Weather Generator52, which is used to create an ensemble of another 40 years of data (Figure S7c in Supplementary Information). Landcover is also obtained from the field survey in 2015 (Figure S7b in Supplementary Information). It is essential for simulating surface water runoff and the input of SOC from plant residues. Different types of landcover provide different fractions that control the surface water runoff velocity, and such fraction is represented by Manning’s coefficient (Table S2 in Supplementary Information). The plant residues include dead leaves, roots, stems, and corn stover after harvest. Here, we estimate plant residues as a function of the Normalized Difference Vegetation Index (NDVI) (Figure S3a) which is obtained from Landsat satellite data (see more information in the Section below). Litter input estimation The NDVI is collected from Landsat satellite data for a full 2 years period, 2016 and 2017. It is spatially divided into three areas, the consolidated gully land, the GLC Watershed, and the Reference Watershed. The spatial distribution of NDVIs for the two watersheds are nearly identical, so we took the spatial means of NDVI for the two watersheds excluding the consolidated gully to represent the natural area (Figure S3a in Supplementary Information). During the growing season, the NDVIs in a natural area (Figure S3(a1) in Supplementary Information) are smaller than the one in the consolidated gully (Figure S3a2 in Supplementary Information). This is because the crop inside the consolidated gully has higher vegetation density. We assume the rate of litterfall has an exponentially increasing relationship with NDVI (Figure S3b in Supplementary Information). As the NDVI increases, the plant residues increase in general. When a plant’s growth slows down, the NDVI increase also slows down and is close to the maximum, but on the other hand, the litterfall rate increases much faster near the end of the growing season. This characterization allows us to fill in the gap due to the lack of litterfall data about the various land vegetation types, including the cornfield in the consolidated gully. Governing equations of SOC transformation The equation below shows the SOC transformation, which is directly affected by litter input and decomposition rate29: $$begin{aligned} frac{partial big ( C_l + C_h + C_b big ) }{partial t} = I_{litter} - big (r_r K_l C_l+r_r K_h C_h big ) end{aligned}$$ (11) where (C_l), (C_h), and (C_b) are defined in Eq. 1; (I_{litter}) is the litter input from the sum of above-ground litter fall and below-ground root-litter ([M L^{-2}T^{-1} ]); (r_r) defines the fraction of decomposed organic carbon to (hbox {CO}_2) ([-]) ((0 le r_r le 1-r_h)), which typically ranges between 0.6 and 0.8; (K_l) and (K_h) are rates of carbon decomposition in fast and slow pool, respectively ([T^{-1}]). They are regulated by soil moisture and carbon-nitrogen ((C/N)) ratio as shown below29: $$begin{aligned}&K_l = varphi (C/N) f_d(theta ) k_l C_b end{aligned}$$ (12a) $$begin{aligned}&K_h =varphi (C/N) f_d(theta ) k_h C_b end{aligned}$$ (12b) where (k_l) and (k_h) represent the rate of decomposition as a simplified term that encompasses different organic components in the litter and humus pool, respectively ([L^3 T^{-1} M^{-1}]); (varphi (C/N)) is a ratio that is from the reduction of the decomposition rate if the immobilization (controlled by nitrogen content) fails to meet the nitrogen demand by the microbes ([-]). (varphi approx 1) in agricultural fields where nitrogen supply is usually sufficient from fertilizers; (f_d(theta )) ([-]) represents the soil moisture effects on decomposition29. The optimal soil moisture condition is the field capacity which provides the highest (f_d)29, meaning that very dry or very wet conditions will result in a smaller (f_d), and hence reduce the decomposition rate. The decomposition rates for litter (or fast) pool and humus (or slow) pool from the equations are (r_r K_l) and (r_r K_h), respectively. In this study, we test the different mean residence times on the surface soil layer (5-cm) by assigning a new decay rate of the decomposition parameter (k_h) in humus (or slow) pool. A complete list of parameters can be found in Supplementary Information Table S2. More

Quantifying structural sensitivity in models with uncertain component functions
In general, we consider a system of the form:

$$begin{aligned} dot{{mathbf{x }}}={mathbf{G}} left( h_1left( {mathbf{x}} right) ,h_2left( {mathbf{x}} right) ,ldots ,h_pleft( {mathbf{x}} right) , f_1left( {mathbf{x}} right) ,f_2left( {mathbf{x}} right) ,ldots ,f_{m-p}left( {mathbf{x}} right) right) , end{aligned}$$
(1)

where ({mathbf {x}}in {mathbb {R}}) is the vector of d state variables, (h_i, f_i:{mathbb {R}}^{d_i}rightarrow {mathbb {R}}) are the m different component functions describing the inflows and outflows of biomass, energy or individuals due to certain biological processes, with ({mathbf{G}} :{mathbb {R}}^mrightarrow {mathbb {R}}^d) being a composition function describing the general topology of the system. We consider that the precise mathematical formulation of the functions (f_i) are known (or at least postulated) with the only related uncertainty being the precise choice of their parameters. The functions (h_1,ldots h_p) are considered to have unspecified functional form. Instead, they are represented by bounds on their derivatives matching the qualitative properties we would expect from such a function. For example, the per-capita reproduction rate of a population is generally decreasing, at least at large population numbers, while a feeding term described by a Holling type II functional response of a predator should be increasing and decelerating. The (h_i) may also have quantitative bounds on their values:

$$begin{aligned} h_i^{text {low}}left( {mathbf{x}} right)1)), an important question remains: which of the unknown functions contribute the most to the degree of structural sensitivity in the system? The degree of structural sensitivity does not distinguish between the various sources of uncertainty and therefore cannot quantify the relative contributions of the unknown functions to the uncertainty in the model dynamics.
To determine the contribution of each unknown function (h_i), one can allow the error terms (left( varepsilon _1,varepsilon _2,ldots ,varepsilon _pright)) to vary with the goal of investigating how the degree of sensitivity varies with them. For the purpose of this section, let us denote the initial error terms by (varepsilon _i^0). We might be tempted to use the dependence on the (varepsilon _i) to perform global optimisation under certain constraints to find the best possible reduction of (left( varepsilon _1,varepsilon _2,ldots ,varepsilon _pright)). However, one should bear in mind that this analysis would depend on the base functions (hat{h}_i) considered. While these functions are ideally fitted to experimental data, they are only accurate within the error terms (varepsilon _i^0). Excessively reducing the (varepsilon _i) will force all admissible functions to conform strongly in their shape to these base functions far beyond their demonstrated accuracy of fit.
The dependence of the degree of sensitivity on (varepsilon _i) should therefore only be evaluated locally by calculating the gradient (left( frac{partial Delta }{partial varepsilon _1},ldots ,frac{partial Delta }{partial varepsilon _p}right) |_{left( varepsilon _1^0,ldots varepsilon _p^0right) }) giving the direction for the best local reduction of the errors. To adjust for the fact that the error terms may be of different orders of magnitude, when handling the vectors of error terms we should use the norm

$$begin{aligned} left| varvec{varepsilon }right| = sqrt{sum _{i=1}^{p} left( frac{varepsilon _i }{varepsilon _i^0}right) ^2}. end{aligned}$$
(9)

Working in this norm, the gradient needs to be weighted by the initial error terms to provide the direction for the best local reduction of the error terms, this is described by the following structural sensitivity gradient.
Definition 2
The structural sensitivity gradient in a model with p unknown functions each having an error of magnitude (varepsilon ^0_i) is defined as

$$begin{aligned} left( -varepsilon _1^0cdot frac{partial Delta }{partial varepsilon _1},ldots ,-varepsilon _p^0cdot frac{partial Delta }{partial varepsilon _p}right) |_{left( varepsilon _1=varepsilon _1^0,ldots ,varepsilon _p=varepsilon _p^0right) }, end{aligned}$$
(10)

where (Delta left( varepsilon _1,ldots ,varepsilon _pright)) is the degree of structural sensitivity of the system considered as a function of the error terms (varepsilon _i).
One possible problem with the structural sensitivity gradient is that the degree of structural sensitivity in the system may not be an increasing function of the magnitude of the errors. Consider the case that the exact system is structurally unstable, e.g. at a bifurcation point. Then no matter how small the error terms are, there may still be very high levels of structural sensitivity, while larger error terms may cause the level of uncertainty to decrease5. In this case, the structural sensitivity gradient will indicate that one or more of the functions has a negative contribution to the uncertainty of the system, and cannot be taken as a basis for sensitivity analysis.
An alternative approach to quantifying the individual impact of unknown functions which avoids this issue is the computation of partial degrees of sensitivity with respect to each (h_k). To do this, we fix every unknown function except (h_k), a set which we denote ({mathbf{H}} _{sim k}), by fixing the (x_j^*), (h_ileft( {mathbf{x}} ^*right)), and (frac{partial h_i}{partial x_j} left( {mathbf{x}} ^* right)) that are consequently determined by the isocline equations. Denoting by (V_k) the cross-sections of V where only (h_k) varies, and the cross-sections for ({mathbf{H}} _{sim k}) by (V_{sim k}), the local partial degree of structural sensitivity can be defined as follows.
Definition 3
The local partial degree of structural sensitivity with respect to (h_k), is the degree of structural sensitivity in the model when (h_k) is unspecified and all other functions (h_i in {mathbf{H}} _{sim k}) are fixed:

$$begin{aligned} Delta _k({mathbf{H}} _{sim k}):= 4 cdot int _{V_{k_{text {stable}}}} rho _{mathbf{H} _kvert {mathbf{H}} _{sim k}} , d{mathbf{H}} _k cdot left( 1 – int _{V_{k_{text {stable}}}} rho _{mathbf{H} _kvert {mathbf{H}} _{sim k}} , d{mathbf{H}} _k right) , end{aligned}$$
(11)

where (rho _{mathbf{H} _kvert {mathbf{H}} _{sim k}}) is the conditional probability density function on ({mathbf{H}} _{sim k}):

$$begin{aligned} rho _{mathbf{H} _kvert {mathbf{H}} _{sim k}} = frac{rho }{int _V rho d{mathbf{H}} _{sim k }}, end{aligned}$$

with (rho) the (joint) probability distribution over V.
The local partial sensitivity (Delta _k({mathbf{H}} _{sim k})) is a function of ({mathbf{H}} _{sim k}) in that it depends upon the particular values at which the elements of (V_{sim k}) are fixed. As with the degree of structural sensitivity, it can be interpreted as either the probability that the stability of the given equilibrium will be different for two independent choices of the function (h_k) when the (h_{sim k}) are fixed at the given values, or in terms of variance as (Delta _k({mathbf{H}} _{sim k})=4cdot {text {Var}}_k(Yvert {mathbf{H}} _{sim k})) (Y is the Bernoulli variable for stability). If the joint probability distribution (rho) is uniform in V, then (Delta _k({mathbf{H}} _{sim k})) can be expressed purely in terms of the fraction of the volume of (V_k) which gives a stable equilibrium:

$$begin{aligned} Delta _k({mathbf{H}} _{sim k}) = 4 cdot frac{int _{V_{k_{text {stable}}}} d{mathbf{H}} _k}{int _{V_k} d{mathbf{H}} _k} cdot left( 1 – frac{int _{V_{k_{text {stable}}}} d{mathbf{H}} _k}{int _{V_k} d{mathbf{H}} _k} right) . end{aligned}$$
(12)

To obtain a global measure for the sensitivity of the model to (h_k), we can take the average of (Delta _k) over (V_{sim k}).
Definition 4
The partial degree of structural sensitivity with respect to (h_k) is given by

$$begin{aligned} bar{Delta }_k := int _{V_{sim k}} rho _{mathbf{H} _{sim k}} cdot Delta _k({mathbf{H}} _{sim k}) , d{mathbf{H}} _{sim k} end{aligned}$$
(13)

where (rho _{mathbf{H} _{sim k}}) is the marginal probability density function of ({mathbf{H}} _{sim k}).
Recalling the variance-based interpretation of the degree of sensitivity, we obtain (bar{Delta }_k = 4cdot E_{sim k}({text {Var}}_kleft( Yvert {mathbf{H}} _{sim k}right) )). In other words, (bar{Delta }_k) gives the scaled average variance when all functions except (h_k) are fixed. We can also relate the partial degree of structural sensitivity to indices used in conventional variance-based sensitivity analysis. Dividing (bar{Delta }_k) by the overall degree of structural sensitivity in the model gives us (frac{bar{Delta }_k}{Delta }=frac{E_{sim k}({text {Var}}_kleft( Yvert {mathbf{H}} _{sim k}right) )}{{text {Var}}(Y)}=S_{T_k}), the total effect index23 of (h_k) on the stability of the equilibrium. This is a measure of the total contribution of (h_k) to the sensitivity—both alone and in conjunction with the other functions ({mathbf{H}} _{sim k}). However, since the space of valid functions V is in general not a hypercube, the functions (h_i) are not independent factors, and a total decomposition of variance is not possible. Indeed, even if the joint probability distribution (rho) is uniform, the marginal probability distribution (rho _{mathbf{H} _{sim k}}) will generally not be: instead it will equal the volume of the corresponding cross-section (V_k) for ({mathbf{H}} _{sim k}), divided by the volume of V. An alternative to using the partial degrees of sensitivity would be to consider the first-order sensitivity indices (S_k=frac{{text {Var}}_kleft( E_{mathbf{H} _{sim k}}left( Yvert h_k right) right) }{{text {Var}}(Y)}). However, these do not take into account possible joint effects of the (h_i) on the structural sensitivity of the system, so a small (S_k) does not indicate that (h_k) is not a source of sensitivity, whereas (bar{Delta }_k=0) means that (h_k) does not contribute to the structural sensitivity in the system at all.
Similarly to the gradient of the total degree of sensitivity (Delta) as a function of the respective error tolerances, the vector (left( -bar{Delta }_{h_1},ldots , -bar{Delta }_{h_p}right)) needs to be scaled by the elements of (varepsilon ^0) to give us the optimal direction of decrease in (Delta) if the error terms (varepsilon _i) are subject to a proportional reduction. This is described by (left( -varepsilon _1^0 bar{Delta }_{h_1},ldots , -varepsilon _p^0bar{Delta }_{h_p}right)).
Outline of an iterative framework of experiments for reducing structural sensitivity
When dealing with partially specified models, an important practical task is the reduction of the overall uncertainty in the system by decreasing the uncertainty in the system processes (i.e. the unknown model functions). Here we propose an iterative process of such a reduction based on improving our empirical knowledge of the uncertain functions (h_k).
As a starting point, we assume that experiments have produced data on the unknown functions (h_1,ldots ,h_p), to which we can fit some base functions (hat{h}_1,ldots ,hat{h}_p) with initial errors (varepsilon _1^0,ldots ,varepsilon _p^0). We assume that it is possible to perform additional experiments on all uncertain processes in order to obtain more data such that the (varepsilon _i) can be decreased, but with the natural constraint that the total error can only be reduced by a magnitude of (0 More