Philippa Kaur

1.Hein, J. R., Mizell, K., Koschinsky, A. & Conrad, T. A. Deep-ocean mineral deposits as a source of critical metals for high- and green-technology applications: Comparison with land-based resources. Ore Geol. Rev. 51, 1–14 (2013).Article 

Google Scholar 
2.Petersen, S. et al. News from the seabed—Geological characteristics and resource potential of deep-sea mineral resources. Mar. Policy 70, 175–187 (2016).Article 

Google Scholar 
3.Dutkiewicz, A., Judge, A. & Müller, R. D. Environmental predictors of deep-sea polymetallic nodule occurrence in the global ocean. Geology 48, 293–297 (2020).ADS 
CAS 
Article 

Google Scholar 
4.Verlaan, P. A. & Cronan, D. S. Origin and variability of resource-grade marine ferromanganese nodules and crusts in the Pacific Ocean: A review of biogeochemical and physical controls. Geochemistry https://doi.org/10.1016/j.chemer.2021.125741 (2021).Article 

Google Scholar 
5.Radziejewska, T. & Stoyanova, V. Abyssal epibenthic megafauna of the Clarion-Clipperton area (NE Pacific): Changes in time and space versus anthropogenic environmental disturbance. Oceanol. Stud. 29, 83–101 (2000).
Google Scholar 
6.Vanreusel, A., Hilario, A., Ribeiro, P. A., Menot, L. & Arbizu, P. M. Threatened by mining, polymetallic nodules are required to preserve abyssal epifauna. Sci. Rep. 6, 1–6 (2016).Article 
CAS 

Google Scholar 
7.Simon-Lledó, E. et al. Ecology of a polymetallic nodule occurrence gradient: Implications for deep-sea mining. Limnol. Oceanogr. 64, 1883–1894 (2019).ADS 
PubMed 
PubMed Central 
Article 

Google Scholar 
8.Washburn, T. W. et al. Patterns of macrofaunal biodiversity across the Clarion-Clipperton zone: An area targeted for seabed mining. Front. Mar. Sci. 8, 626571 (2021).Article 

Google Scholar 
9.Bonifácio, P., Martinez Arbizu, P. & Menot, L. Alpha and beta diversity patterns of polychaete assemblages across the nodule province of the eastern Clarion-Clipperton Fracture Zone (equatorial Pacific). Biogeosciences 17, 865–886 (2020).ADS 
Article 

Google Scholar 
10.Ansari, Z. A. Distribution of deep-sea benthos in the proposed mining area of Central Indian Basin. Mar. Georesour. Geotechnol. 18, 201–207 (2000).Article 

Google Scholar 
11.Pasotti, F. et al. A local scale analysis of manganese nodules influence on the Clarion-Clipperton Fracture Zone macrobenthos. Deep Sea Res. Part Oceanogr. Res. Pap. 168 (2021).12.Hauquier, F. et al. Geographic distribution of free-living marine nematodes in the Clarion-Clipperton Zone: Implications for future deep-sea mining scenarios. Biogeosciences 16, 3475–3489 (2019).ADS 
CAS 
Article 

Google Scholar 
13.Kuhn, T., Uhlenkott, K., Vink, A., Rühlemann, C. & Martinez Arbizu, P. Manganese nodule fields from the Northeast Pacific as benthic habitats. In Seafloor Geomorphology as Benthic Habitat 2nd edn (eds Harris, P. T. & Baker, E.) 933–947 (Elsevier, 2020). https://doi.org/10.1016/B978-0-12-814960-7.00058-0.Chapter 

Google Scholar 
14.Miljutina, M. A., Miljutin, D. M., Mahatma, R. & Galéron, J. Deep-sea nematode assemblages of the Clarion-Clipperton Nodule Province (Tropical North-Eastern Pacific). Mar. Biodivers. 40, 1–15 (2010).Article 

Google Scholar 
15.Mahatma, R. Meiofauna Communities of the Pacific Nodule Province: Abundance, Diversity and Community Structure (University of Oldenburg, 2009).
Google Scholar 
16.Singh, R. et al. Nematode communities inhabiting the soft deep-sea sediment in polymetallic nodule fields: Do they differ from those in the nodule-free abyssal areas?. Mar. Biol. Res. 12, 1–15 (2016).Article 

Google Scholar 
17.Thiel, H., Schriever, G., Bussau, C. & Borowski, C. Manganese nodule crevice fauna. Deep Sea Res. Part Oceanogr. Res. Pap. 40, 419–423 (1993).ADS 
Article 

Google Scholar 
18.Bussau, C., Schriever, G. & Thiel, H. Evaluation of abyssal metazoan meiofauna from a manganese nodule area of the Eastern South Pacific. Vie Milieu 45, 39–48 (1995).
Google Scholar 
19.Oebius, H. U., Becker, H. J., Rolinski, S. & Jankowski, J. A. Parametrization and evaluation of marine environmental impacts produced by deep-sea manganese nodule mining. Deep Sea Res. Part II Top. Stud. Oceanogr. 48, 3453–3467 (2001).ADS 
CAS 
Article 

Google Scholar 
20.Levin, L. A. et al. Defining “serious harm” to the marine environment in the context of deep-seabed mining. Mar. Policy 74, 245–259 (2016).Article 

Google Scholar 
21.Global Sea Mineral Resources. Environmental Impact Statement—Small-scale testing of nodule collector components on the seafloor of the Clarion-Clipperton Fracture Zone and its environmental impact. 337 (2018).22.Durden, J. M. et al. A procedural framework for robust environmental management of deep-sea mining projects using a conceptual model. Mar. Policy 84, 193–201 (2017).Article 

Google Scholar 
23.Jones, D. O. B. et al. Biological responses to disturbance from simulated deep-sea polymetallic nodule mining. PLoS One 12, e0171750 (2017).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
24.Jones, D. O. B., Ardron, J. A., Colaço, A. & Durden, J. M. Environmental considerations for impact and preservation reference zones for deep-sea polymetallic nodule mining. Mar. Policy https://doi.org/10.1016/j.marpol.2018.10.025 (2018).Article 

Google Scholar 
25.Boschen, R. E. et al. A primer for use of genetic tools in selecting and testing the suitability of set-aside sites protected from deep-sea seafloor massive sulfide mining activities. Ocean Coast. Manag. 122, 37–48 (2016).Article 

Google Scholar 
26.Boucher, G. & Lambshead, P. J. D. Ecological biodiversity of marine nematodes in samples from temperate, tropical and deep-sea regions. Conserv. Biol. 9, 1594–1604 (1995).Article 

Google Scholar 
27.Ramirez-Llodra, E. et al. Deep, diverse and definitely different: Unique attributes of the world’s largest ecosystem. Biogeosciences 7, 2851–2899 (2010).ADS 
Article 

Google Scholar 
28.Rex, M. A. & Etter, R. J. Deep-Sea Biodiversity: Pattern and Scale (Harvard University Press, 2010).
Google Scholar 
29.Paterson, G. L. J. et al. Biogeography and connectivity in deep-sea habitats with mineral resource potential: A gap analysis. Deliverable 4.2. MIDAS (2014).30.Christodoulou, M., O’Hara, T. D., Hugall, A. F. & Arbizu, P. M. Dark ophiuroid biodiversity in a prospective abyssal mine field. Curr. Biol. 29, 3909–3912 (2019).PubMed 
CAS 
Article 

Google Scholar 
31.Amon, D. J. et al. Insights into the abundance and diversity of abyssal megafauna in a polymetallic-nodule region in the eastern Clarion-Clipperton Zone. Sci. Rep. 6, 30492 (2016).ADS 
PubMed 
PubMed Central 
CAS 
Article 

Google Scholar 
32.Goineau, A. & Gooday, A. J. Diversity and spatial patterns of foraminiferal assemblages in the eastern Clarion-Clipperton zone (abyssal eastern equatorial Pacific). Deep Sea Res. Part Oceanogr. Res. Pap. 149, 103036 (2019).Article 

Google Scholar 
33.Macheriotou, L., Rigaux, A., Derycke, S. & Vanreusel, A. Phylogenetic clustering and rarity imply risk of local species extinction in prospective deep-sea mining areas of the Clarion-Clipperton Fracture Zone. Proc. R. Soc. B Biol. Sci. 287, 20192666 (2020).Article 

Google Scholar 
34.Błażewicz, M., Jóźwiak, P., Menot, L. & Pabis, K. High species richness and unique composition of the tanaidacean communities associated with five areas in the Pacific polymetallic nodule fields. Prog. Oceanogr. 176, 102141 (2019).Article 

Google Scholar 
35.Janssen, A. et al. A reverse taxonomic approach to assess macrofaunal distribution patterns in abyssal pacific polymetallic nodule fields. PLoS One 10, e0117790 (2015).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
36.Soetaert, K. & Heip, C. Sample-size dependence of diversity indexes and the determination of sufficient sample size in a high-diversity deep-sea environment. Mar. Ecol. Prog. Ser. 59, 305–307 (1990).ADS 
Article 

Google Scholar 
37.Rose, A. et al. A method for comparing within-core alpha diversity values from repeated multicorer samplings, shown for abyssal Harpacticoida (Crustacea: Copepoda) from the Angola Basin. Org. Divers. Evol. 5, 3–17 (2005).Article 

Google Scholar 
38.George, K. H. et al. Community structure and species diversity of Harpacticoida (Crustacea: Copepoda) at two sites in the deep sea of the Angola Basin (Southeast Atlantic). Org. Divers. Evol. 14, 57–73 (2014).Article 

Google Scholar 
39.Mouillot, D. et al. Rare species support vulnerable functions in high-diversity ecosystems. PLoS Biol. 11, e1001569 (2013).PubMed 
PubMed Central 
CAS 
Article 

Google Scholar 
40.Naeem, S. Species redundancy and ecosystem reliability. Conserv. Biol. 12, 39–45 (1998).Article 

Google Scholar 
41.Turner, P. J., Campbell, L. M. & Van Dover, C. L. Stakeholder perspectives on the importance of rare-species research for deep-sea environmental management. Deep Sea Res. Part Oceanogr. Res. Pap. 125, 129–134 (2017).ADS 
Article 

Google Scholar 
42.Drury, W. H. Rare species. Biol. Conserv. 6, 162–169 (1974).Article 

Google Scholar 
43.Anderson, M. J., Gorley, R. N. & Clarke, K. R. PERMANOVA+ for PRIMER: Guide for Software and Statistical Methods (Primer-E Ltd, 2008).
Google Scholar 
44.Gollner, S. et al. Resilience of benthic deep-sea fauna to mining activities. Mar. Environ. Res. https://doi.org/10.1016/j.marenvres.2017.04.010 (2017).Article 
PubMed 

Google Scholar 
45.Glover, A. G. et al. Polychaete species diversity in the central Pacific abyss: Local and regional patterns, and relationships with productivity. Mar. Ecol. Prog. Ser. 240, 157–170 (2002).ADS 
Article 

Google Scholar 
46.Rosli, N., Leduc, D., Rowden, A. & Robert, K. Review of recent trends in ecological studies of deep-sea meiofauna, with focus on patterns and processes at small to regional spatial scales. Mar. Biodivers. 18, 13–34 (2018).Article 

Google Scholar 
47.Gallucci, F., Moens, T. & Fonseca, G. Small-scale spatial patterns of meiobenthos in the Arctic deep sea. Mar. Biodivers. 39, 9–25 (2009).Article 

Google Scholar 
48.Wieser, W. Die Beziehung zwischen Mundhöhlengestalt, Ernährungsweise und Vorkommen bei freilebenden marinen Nematoden Eine ökologisch-morphologische Studie. Ark. För Zool. 4, 439–483 (1953).
Google Scholar 
49.Leduc, D. Description of Oncholaimus moanae sp. nov. (Nematoda: Oncholaimidae), with notes on feeding ecology based on isotopic and fatty acid composition. J. Mar. Biol. Assoc. U. K. 89, 337–344 (2008).Article 
CAS 

Google Scholar 
50.Pape, E., van Oevelen, D., Moodley, L., Soetaert, K. & Vanreusel, A. Nematode feeding strategies and the fate of dissolved organic matter carbon in different deep-sea sedimentary environments. Deep Sea Res. Part Oceanogr. Res. Pap. 80, 94–110 (2013).ADS 
CAS 
Article 

Google Scholar 
51.Schuelke, T., Pereira, T. J., Hardy, S. M. & Bik, H. M. Nematode-associated microbial taxa do not correlate with host phylogeny, geographic region or feeding morphology in marine sediment habitats. Mol. Ecol. 27, 1930–1951 (2018).PubMed 
Article 

Google Scholar 
52.Tully, B. J. & Heidelberg, J. F. Microbial communities associated with ferromanganese nodules and the surrounding sediments. Extreme Microbiol. 4, 161 (2013).
Google Scholar 
53.Blöthe, M. et al. Manganese-cycling microbial communities inside deep-sea manganese nodules. Environ. Sci. Technol. 49, 7692–7700 (2015).ADS 
PubMed 
Article 
CAS 

Google Scholar 
54.Maybury, C. Crevice Foraminifera from abyssal South East Pacific manganese nodules. In Microfossils and Oceanic Environments (eds Moguilevsky, A. & Whatley, R.) (University of Wales, 1996).
Google Scholar 
55.Pape, E., Bezerra, T. N., Hauquier, F. & Vanreusel, A. Limited spatial and temporal variability in meiofauna and nematode communities at distant but environmentally similar sites in an area of interest for deep-sea mining. Front. Mar. Sci. 4, 205 (2017).Article 

Google Scholar 
56.Uhlenkott, K., Vink, A., Kuhn, T. & Arbizu, P. M. Meiofauna in a potential deep-sea mining area—Influence of temporal and spatial variability on small-scale abundance models. Diversity 13, 3 (2021).CAS 
Article 

Google Scholar 
57.Veillette, J., Juniper, S. K., Gooday, A. J. & Sarrazin, J. Influence of surface texture and microhabitat heterogeneity in structuring nodule faunal communities. Deep Sea Res. Part Oceanogr. Res. Pap. 54, 1936–1943 (2007).ADS 
Article 

Google Scholar 
58.Tilot, V., Ormond, R., Moreno Navas, J. & Catalá, T. S. The Benthic Megafaunal Assemblages of the CCZ (Eastern Pacific) and an approach to their management in the face of threatened anthropogenic impacts. Front. Mar. Sci. 5, 7 (2018).Article 

Google Scholar 
59.ISA. Recommendations for the guidance of contractors for the assessment of the possible environmental impacts arising from exploration for marine minerals in the Area (2020).60.ISA. Draft regulations on exploitation of mineral resources in the Area (2019).61.ISA. Environmental Management Plan for the Clarion-Clipperton Zone (2011).62.Wedding, L. M. et al. From principles to practice: A spatial approach to systematic conservation planning in the deep sea. Proc. R. Soc. B Biol. Sci. 280, 20131684 (2013).CAS 
Article 

Google Scholar 
63.ISA. Deep CCZ Biodiversity Synthesis Workshop Report. 206 (2020).64.McQuaid, K. A. et al. Using habitat classification to assess representativity of a protected area network in a large, data-poor area targeted for deep-sea mining. Front. Mar. Sci. 7, 558860 (2020).Article 

Google Scholar 
65.Mullineaux, L. S. The role of settlement in structuring a hard-substratum community in the deep sea. J. Exp. Mar. Biol. Ecol. 120, 247–261 (1988).Article 

Google Scholar 
66.Cuvelier, D. et al. Potential mitigation and restoration actions in ecosystems impacted by seabed mining. Front. Mar. Sci. 5, 467 (2018).Article 

Google Scholar 
67.De Smet, B. et al. The community structure of deep-sea macrofauna associated with polymetallic nodules in the eastern part of the Clarion-Clipperton fracture zone. Front. Mar. Sci. 4, 103 (2017).
Google Scholar 
68.Bezerra, T. N. et al. Nemys: World Database of Nematodes. http://nemys.ugent.be. https://doi.org/10.14284/366 (2021).69.George, K.-H. Gemeinschaftsanalytische Untersuchungen der Harpacticoidenfauna der Magellanregion, sowie erste similaritätsanalytische Vergleiche mit Assoziationen aus der Antarktis = Community analysis of the harpacticoid fauna of the Magellan Region, as well as first comparisons with antarctic associations, based on similarity analyses. Berichte Zur Polarforsch. Rep. Polar Res. 327, 1–187 (1999).
Google Scholar 
70.Moens, T. & Vincx, M. Observations on the feeding ecology of estuarine nematodes. J. Mar. Biol. Assoc. U. K. 77, 211–227 (1997).Article 

Google Scholar 
71.Guilini, K., Van Oevelen, D., Soetaert, K., Middelburg, J. J. & Vanreusel, A. Nutritional importance of benthic bacteria for deep-sea nematodes from the Arctic ice margin: Results of an isotope tracer experiment. Limnol. Oceanogr. 55, 1977–1989 (2010).ADS 
CAS 
Article 

Google Scholar 
72.Clarke, K. & Gorley, R. PRIMER v6: User Manual/Tutorial (Primer-E Ltd, 2006).
Google Scholar 
73.R Core Team. R: A Language and Environment for Statistical Computing (R Foundation for Statistical Computing, 2019).
Google Scholar 
74.Wickham, H. ggplot2: Elegant Graphics for Data Analysis (Springer, 2016).MATH 
Book 

Google Scholar 
75.Wilke, C. O. cowplot: Streamlined Plot Theme and Plot Annotations for ‘ggplot2’ (2019).76.Oksanen, J. et al. vegan: Community Ecology Package (2019).77.Martinez Arbizu, P. M. pairwiseAdonis: Pairwise Multilevel Comparison using Adonis (2017).78.Chao, A. et al. Rarefaction and extrapolation with Hill numbers: A framework for sampling and estimation in species diversity studies. Ecol. Monogr. 84, 45–67 (2014).Article 

Google Scholar 
79.Hsieh, T. C. & Chao, A. Package iNEXT 2.0.19: Interpolation and extrapolation of species diversity (2019).80.Schenker, N. & Gentleman, J. F. On judging the significance of differences by examining the overlap between confidence intervals. Am. Stat. 55, 182–186 (2001).MathSciNet 
Article 

Google Scholar 
81.Gehlenborg, N. UpSetR: A More Scalable Alternative to Venn and Euler Diagrams for Visualizing Intersecting Sets (2019).82.Simpson, G. L. permute: Functions for Generating Restricted Permutations of Data (2019).83.Baselga, A., Orme, D., Villeger, S., Bortoli, J. D. & Leprieur, F. betapart: Partitioning Beta Diversity into Turnover and Nestedness Components (2018). More

Sequencing of Shorea leprosula genomeSample collectionLeaf samples of S. leprosula were obtained from a reproductively mature (diameter at breast height, 50 cm) diploid tree B1_19 (DNA ID 214) grown in the Dipterocarp Arboretum, Forest Research Institute Malaysia (FRIM).DNA extractionGenomic DNA was extracted from leaf samples using the 2% cetyltrimethylammonium bromide (CTAB) method90 and purified using a High Pure PCR Template Purification kit (Roche).Library preparation and sequencingPaired-end (170, 500, and 800 bp) and mate-pair (2 kb) genomic libraries were prepared using a TruSeq DNA Library Preparation kit (Illumina) and a Mate Pair Library Preparation kit (Illumina), respectively. Mate-pair libraries with larger insert sizes were constructed using a Nextera Mate Pair Library Preparation kit (Illumina). Ten micrograms of genomic DNA were tagmented in a 400 μl reaction and fractionated using SageELF, with the recovery of 11 fractions with 3–16+ kb. Each fraction was circularized and fragmented with a Covaris S2. Biotin-containing fragments were purified using Dynabeads M-280 streptavidin beads. Sequencing adapters (KAPA TruSeq Adapter kit) were attached using a KAPA Hyper Prep kit. The libraries were amplified for 10–13 cycles and purified with 0.8× AMpure XP. DNA libraries were then sequenced (~388× coverage) using Illumina HiSeq2000 (TruSeq libraries) and HiSeq2500 (Nextera libraries) at the Functional Genomics Center Zurich (FGCZ), University of Zurich, Switzerland (Supplementary Table 1).Genome assemblyAdapters and low-quality bases for all paired-end and mate-pair reads were removed using Trimmomatic91. The filtered paired-end reads of the 170 bp library were used to identify the genome size using k-mer distribution generated by Jellyfish92 that was implemented in the scripts by Joseph Ryan42. The raw R1 reads from paired-end 170 and 800 bp libraries (clipped at 95 bp, representing about 70 genome equivalents) were used to estimate the heterozygosity using KAT43 with a k-mer size of 23 nt. De novo genome assembly of all reads was performed using ALLPATHSLG assembler v5248840.Assembly verification and assessment of the assembled genomeAssembly validationTo validate the genome assembly, we mapped (i) the short reads used for the genome assembly, (ii) scanned the assembly for the presence of single-copy orthologs, and (iii) mapped transcriptome sequences obtained from seven organs.Assembly verification by mapping of short readsFor each library used for genome assembly, all trimmed reads were aligned to the assembled S. leprosula genome using Burrows–Wheeler Aligner (BWA) v0.7.1293. Then, mapping ratio was calculated for each BAM file using Samtools94 with “flagstat” command.Identification of highly conserved single-copy orthologsBUSCO v3.1.042 was run with the Embryophyta dataset and Arabidopsis as the species for AUGUSTUS prediction (see subsection below “Protein-coding gene prediction”).Assembly verification by mapping transcriptome sequencesFor mapping transcriptome sequences, samples of seven organs (leaf bud, flower bud, flower, inner bark, small seed, large seed, and calyx) were obtained from the S. leprosula individual used for the genome sequencing (Supplementary Table 2). Total RNA was extracted from each sample using RNeasy Plant Mini Kit (Qiagen) and it was treated with Turbo DNase I (Takara). Library preparation was carried out using a TruSeq RNA Library Preparation kit (Illumina). Paired-end sequencing was conducted for all the libraries using Illumina HiSeq2000 at the FGCZ, University of Zurich, Switzerland. Adapters and low-quality bases for all paired-end reads were removed using Trimmomatic. The trimmed sequences of each library were mapped onto the assembled genome using STAR aligner v2.4.2a95, and mapping ratio was obtained from the output file of STAR.Genome annotationRepeat sequence analysisBoth homology-based and de novo prediction analyses were used to identify the repeat content in the S. leprosula assembly. For the homology-based analysis, we used Repbase (version 20120418) to perform a TE search with RepeatMasker (4.0.5) and the WuBlast search engine. For the de novo prediction analysis, we used RepeatModeler to construct a TE library. Elements within the library were then classified by homology to Repbase sequences (see subsection below “Preparation of repeat sequences for evidence-based gene prediction”).Protein-coding gene predictionS. leprosula protein-coding genes were predicted by AUGUSTUS v3.245. For ab initio gene prediction, we used a pre-trained A. thaliana metaparameter implemented in AUGUSTUS. For the evidence-based gene prediction, we used the information of exon, intron and repeat sequences of S. leprosula as hints for the AUGUSTUS gene prediction. The details of the preparation of the hints were described in the following subsections.Preparation of repeat sequences for evidence-based gene predictionWe used RepeatModeler to construct a de novo library of repeated sequences in the S. leprosula assembly. Then, using RepeatMasker, we generated a file containing the information of the positions of repeat sequences in the S. leprosula genome based on the RepeatModeler library. Elements within the library were then classified by homology to Repbase sequences. Finally, the hint file for repeat sequences in GFF format was prepared using the two scripts, “10_makeGffRm.pl” and “12_makeTeHints.pl”, stored in https://gitlab.com/rbrisk/ahalassembly.Preparation of the exon and intron information for evidence-based gene predictionTo obtain the exon and intron hints, we used the mapping data of RNA-seq obtained from seven organs of the sequenced S. leprosula individual as described above. First, we merged all the mapping data stored in different BAM files into a single BAM file using SAMtools. Then, we prepared the intron hint file in GFF format using the, “bam2hints” script of AUGUSTUS. The exon hint file was also generated from the merged BAM file using the two AUGUSTUS scripts, “bam2wig” and “wig2hints.pl”. To conduct evidence-based gene prediction with AUGUSTUS, the three hint files (repeat sequences, intron and exon) described above were merged into a single file in GFF format.BUSCO analysisGenome annotation completeness were assessed with BUSCO v3.1.044 using the Embryophyta odb9 dataset composed of 1440 universal Embryophyta single-copy genes. We referred to these 1440 genes as core genes in the main text.Comparison with the proteome of Theobroma cacao
T. cacao’s gene models18 were downloaded from Phytozome 11 (https://phytozome.jgi.doe.gov/pz/portal.html). Then, comparison was conducted with BLASTP96 using the T. cacao proteomes as the BLAST database (E-value cutoff: 1.0E-10). Only the best hit was stored for each gene. We considered these best hits of the T. cacao genes as orthologs of the S. leprosula genes. When the T. cacao orthologs were identified by the BLASTP search, the orthologs of A. thaliana were defined based on the T. cacao-A. thaliana orthologous information provided by Phytozome 11 (Supplementary Table 4). When the T. cacao orthologs were not identified, the orthologs of A. thaliana were searched by BLASTP (E-value cutoff: 1.0E-10) using the A. thaliana proteomes obtained from TAIR 10 (https://www.arabidopsis.org) as the BLAST database.Synteny analysisBased on the result of the above BLASTP searches, we assessed synteny between the S. leprosula scaffolds and the T. cacao chromosomes using MCScanX97. Genome information of T. cacao in GFF format was also obtained from Phytozome 11 as described above, which was used as an input file for MCScanX.Assessment of the genome assemblyPopulation data and other dipterocarp speciesTo assess whether the genome assembly could be used as a reference for the S. leprosula individuals from various populations, we checked mapping ratio, SNP positions, and admixture using the distribution-wide S. leprosula samples. Similarly, to assess whether the S. leprosula assembly could be used as a reference for aligning data from closely related species and determining their mapping ratios. For interspecific analysis, the following three Dipterocarpoideae species: S. platycarpa, D. aromatica, and N. heimii were used (Supplementary Table 7).Sample collection and DNA extractionLeaf samples of 19 S. leprosula individuals from different populations and three other dipterocarp species (S. platycarpa, D. aromatica, and N. heimii) were used as described in Supplementary Tables 6 and 7. Genomic DNA was extracted using the same method as described above.Library preparation and sequencingPaired-end genomic libraries (200 bp) were prepared using a TruSeq DNA Library Preparation kit (Illumina). DNA libraries were then sequenced (~16× coverage each) using Illumina HiSeq2000.Mapping and SNP callingAdapters and low-quality bases from resequencing reads were removed using Trimmomatic. All trimmed reads were then mapped and aligned to the S. leprosula assembly using BWA. Variants were called using GATK v3.598. Duplicated reads were marked using Picard 2.6.0. Within GATK, HaplotypeCaller was used to identify variants for each sample by generating an intermediate genomic variant call format (gVCF). Subsequently, gVCF files were merged using GenotypeGVCFs to produce a raw VCF file containing SNPs and INDELs. Low-quality variants were removed from the raw VCF file by applying the hard filters implemented in GATK. Variants with genotype quality (GQ)  More

Host and gut content samplingA total of 94 individuals (42 A. speciosus, 9 A. argenteus, and 43 M. rufocanus) were captured from four sites within the Kamikawa Chubu national forest in the central area on the island of Hokkaido, Japan (Supplementary Table S1), and a total of 280 gut content (from the small intestine, cecum, and colon) and fecal matter (from the rectum) samples were collected for microbiome analysis (Supplementary Table S2). Based on 16S rRNA amplicon sequencing using Illumina Miseq, a total of 12,286,171 paired-end reads were obtained after quality filtering and chimeric sequence removal. There was an average of 43,879 reads per sample, although it varied among species and gut region (Supplementary Table S3).Within host species/among gut region gut microbiota alpha diversityAlpha diversity of the gut microbiota in the small intestine was significantly lower than the rectum, colon, and cecum in all three host species based on Shannon diversity, Faith’s PD, evenness, and number of ASVs as expected (GLME: all p  0.05; Fig. 1, Supplementary Fig. S2, Supplementary Tables S4–S7). Males had significantly higher alpha diversity within all gut regions of A. speciosus while female A. argenteus had significantly higher alpha diversity as compared to males (GLME, all p  0.05; Supplementary Tables S4–S7) while age had no effect in any gut region of any rodent species (GLME: all p  > 0.05; Supplementary Tables S4–S7).Figure 1Alpha diversity within each gut region of each species based on (a) Shannon diversity and (b) Faith’s PD. Dashed lines separate host species.Full size imageAmong host species alpha diversityMyodes rufocanus had significantly higher alpha diversity in all four gut regions as compared to both A. speciosus and A. argenteus based on all four diversity measurements (GLME: all p  More

Isolation and characterization of bean dreg-degrading strainsA 1362-bp amplification fragment of 16S rDNA was obtained by PCR (GenBank accession number MW406939). This sequence was compared with others in the GenBank database, aligning the 16S rDNA sequences with several Geobacillus sp. strains and constructed a phylogenetic tree (Fig. 2a). The phylogenetic tree clearly showed that strain GT-02 belongs to the G.toebii branch and was similar to G.toebii R-32652, G.toebii NBRC 107807, and G.toebii SK-1 with 99.78%, 99.63% and 99.05% similarities, respectively. According to the study described previously, G.toebii was a gram-positive, aerobic rod and motile bacterial26. G.toebii could produce acid from inositol and gas from nitrate. G.toebii could hydrolysis casein and utilize n-alkanes as carbon source27.Figure 2(a) Phylogenetic tree based on 16S rDNA gene sequences from related species of the genus Geobacillus constructed using the neighbour-joining method with 1000 bootstrap replicates. Branch length is indicated at each node. (b) The growth curve of strain GT-02 with temperature. (c) The growth curve of strain GT-02 with pH.Full size imageThe growth characteristics of strain GT-02, such as temperature and pH values, were investigated. The bacterial strain could grow within a range of 40–75 °C and pH 6.50–9.50, and the optimum temperature and pH were 65 °C and 7.50, respectively (Fig. 2b,c). Compared to other G.toebii strains, the maximum growth temperature and pH of strains R-32652 and SK-1 were 70 °C and 9.0026,28, respectively. These results showed that strain GT-02 was more resistant to high temperature and alkalinity. Fermentation temperature above 70 °C could effectively inactivate harmful microorganisms in organic solid waste12. Therefore, the fermentation temperature was set at 70 °C in this study.Changes in the composition of bean dregs during fermentationChanges in GI, TOC and TN of bean dregs during fermentationThe GI is traditionally used to evaluate the phytotoxicity and maturity of organic fertilizer12. As shown in Fig. 3a, both groups of experiments reached the standard of maturity (GI ≥ 85.00%). Therefore, the fermentation was terminated in five days. In the initial stage of fermentation, the GI of CK dropped to 51.85% on day 2, and the GI of T1 dropped to 41.98% on day 1. Phytotoxicity, which is usually caused by various heavy metals and low-molecular-weight substances, such as NH3 and organic acids, can reduce seed germination and inhibit root development29. During fermentation, bean dregs might produce NH3, organic acids and other substances, which could trigger a decrease in the GI. The GI of T1 showed a clear decrease, which was likely due to the production of toxic organic acids and might also explain the decrease in pH observed in T1 (Fig. 3d). Due to the degradation of organic acids, the GI of T1 increased to 95.06% on the third day and continued to increase to more than 100.00%, whereas in CK, the GI only reached 86.42% at the end of the fermentation. These results revealed that the maturity of T1 on day 3 was markedly higher than that of CK on day 5 and thus suggest that G.toebii can significantly enhance the fermentation efficiency by accelerating the maturation process and thus reducing the thermophilic fermentation period from 5 to 3 days.Figure 3Profiles of GI (a), TOC (b), TN (c), pH (d) and EC (e) during the fermentation process of CK and T1. The data represent the means ± standard deviations from three measurements.Full size imageTOC is usually used as an energy source by microorganisms30. The TOC loss in both CK and T1 increased during fermentation (Fig. 3b). The reduction of TOC was mainly caused by the production of carbon dioxide from bacterial respiration. The rate of TOC loss in T1 was higher than that in CK. At the end of the fermentation, the TOC loss of T1 was 11.78% higher than that in CK. Because of the addition of G.toebii, bacterial metabolism in T1 was more active, and organic degradation was faster.The TN loss in both CK and T1 also showed an upward trend (Fig. 3c). The loss of TN was mainly caused by the volatilization of ammonia nitrogen31. The rate of TN loss in T1 increased more than that of CK group. After fermentation (day 5), the TN loss in T1 was 6.83% higher than that of CK. The mineralization in T1 was more active and thus ammonia nitrogen was more, which was easy to cause volatilization. However, the bean dregs in CK were mature on the 5th day, while those in T1 were on the 3rd day. At this time, the TN loss of mature bean dregs in T1 was 5.66% lower than that in CK, which indicated that the bean dregs lost less nitrogen source when they reached the standard of maturity after the addition of G.toebii.Changes in pH and EC of bean dregs during fermentationThe variation in pH observed during fermentation is due to the interaction between inorganic nitrogen and organic acids produced by the decomposition of organic matter32. As shown in Fig. 3d, the pH of CK gradually increased to 8.72 at the end of the fermentation. The ammonification process and the release of free NH3 during organic matter (OM) degradation lead to increases in pH33. The pH of T1 decreased to 5.73 on day 1, which was due to the formation of more organic acids than CK, and then increased to 8.76 on day 2, which was due to acid consumption and ammonia formation. Figure 2c showed that GT-02 could hardly grow when the pH was lower than 6.00, but the heterogeneity of solid fermentation provided a possible living environment for the growth of GT-02. Subsequently, the pH of T1 slowly decreased to 8.10 due to ammonia volatilization or ammonia conversion. These study findings showed that the pH value of the fermentation process was significantly affected by the addition of GT-02. G.toebii can produce abundant high-temperature enzymes, such as amylase, protease, cellulase, xylanase, and mannanase17, which explains why the ammonification process was faster in T1 than in CK and thus the higher pH was found in T1.The EC, which is a measure of the total ion concentration, describes changes in the levels of organic and inorganic ions such as SO42−, Na+, NH4+, K+, Cl−, and NO3− during the fermentation process34. As shown in Fig. 3e, the EC of the two groups increased significantly during fermentation process (P  More

1.Collins, M. et al. Long-term climate change: Projections, commitments and irreversibility. In Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change (Cambridge University Press, 2013).2.Kawahata, H. et al. Perspective of the response by marine calcifiers to global warming and ocean acidification –Behavior of corals and foraminifers in the high CO2 world in “hot house”. Prog. Earth Planet Sci. 6, 5 (2019).Article 

Google Scholar 
3.Kroeker, K. J., Kordas, R. L., Crim, R. N. & Singh, G. G. Meta-analysis reveals negative yet variable effects of ocean acidification on marine organisms. Ecol. Lett. 13, 1419–1434 (2010).Article 

Google Scholar 
4.Orr, J. C. et al. Anthropogenic ocean acidification over the twenty-first century and its impact on calcifying organisms. Nature 437, 681–686 (2005).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
5.Schiebel, R. Planktic foraminiferal sedimentation and the marine calcite budget. Glob. Biogeochem. Cycles 16, 1065 (2002).ADS 
Article 
CAS 

Google Scholar 
6.Keul, N., Langer, G., de Nooijer, L. J. & Bijma, J. Effect of ocean acidification on the benthic foraminifera Ammonia sp. is caused by a decrease in carbonate ion concentration. Biogeosciences 10, 6185–6198 (2013).ADS 
CAS 
Article 

Google Scholar 
7.Doo, S. S., Fujita, K., Byrne, M. & Uthicke, S. Fate of calcifying tropical symbiont-bearing large benthic Foraminifera: Living sands in a changing ocean. Biol. Bull. 226, 169–186 (2014).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
8.Prazeres, M., Uthicke, S. & Pandolfi, J. M. Ocean acidification induces biochemical and morphological changes in the calcification process of large benthic foraminifera. Proc. R. Soc. B 282, 20142782 (2015).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
9.Iwasaki, S. et al. Sensitivity of planktic foraminiferal test bulk density to ocean acidification. Sci. Rep. 9, 9803 (2019).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
10.Hohenegger, J., Kinoshita, S., Briguglio, A., Eder, W. & Wöger, J. Lunar cycles and rainy seasons drive growth and reproduction in nummulitid foraminifera, important producers of carbonate buildups. Sci. Rep. 9, 8286 (2019).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
11.Kinoshita, S. et al. Temperature effects on the shell growth of a larger benthic foraminifer (Sorites orbiculus): Results from culture experiments and micro X-ray computed tomography. Mar. Micropaleontol. 163, 101960 (2021).ADS 
Article 

Google Scholar 
12.Fujita, K. & Fujimura, H. Organic and inorganic carbon production by algal symbiont-bearing foraminifera on northwest Pacific coral-reef flat. J. Foraminifer. Res. 38, 117–126 (2008).Article 

Google Scholar 
13.Raja, R., Saraswati, P. K., Rogers, K. & Iwao, K. Magnesium and strontium compositions of recent symbiont-bearing benthic foraminifera. Mar. Micropaleontol. 58, 31–44 (2005).ADS 
Article 

Google Scholar 
14.Narayan, G. R. et al. Response of large benthic foraminifera to climate and local changes: Implications for future carbonate production. Sedimentology. 12858. https://doi.org/10.1111/sed.12858 (2021).
15.Morse, J. W., Andersson, A. J. & Mackenzie, F. T. Initial responses of carbonate-rich shelf sediments to rising atmospheric pCO2 and “ocean acidification”: Role of high Mg-calcites. Geochim. Cosmochim. Acta 70, 5814–5830 (2006).ADS 
CAS 
Article 

Google Scholar 
16.Fujita, K., Nishi, H. & Saito, T. Population dynamics of Marginopora kudakajimaensis Gudmundsson (Foraminifera: Soritidae) in the Ryukyu Islands, the tropical northwest Pacific. Mar. Micropaleontol. 38, 267–284 (2000).ADS 
Article 

Google Scholar 
17.Kuroyanagi, A., Kawahata, H., Suzuki, A., Fujita, K. & Irie, T. Impacts of ocean acidification on large benthic foraminifers: Results from laboratory experiments. Mar. Micropaleontol. 73, 190–195 (2009).ADS 
Article 

Google Scholar 
18.Barker, S. & Elderfield, H. Foraminiferal calcification response to glacial–interglacial changes in atmospheric CO2. Science 297, 833–836 (2002).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
19.Osborne, E. B. et al. Calcification of the planktonic foraminifera Globigerina bulloides and carbonate ion concentration: Results from the Santa Barbara Basin. Paleoceanography 31, 1083–1102 (2016).ADS 
Article 

Google Scholar 
20.Mollica, N. R. et al. Ocean acidification affects coral growth by reducing skeletal density. Proc. Natl. Acad. Sci. 115, 1754–1759 (2018).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
21.Schmidt, C., Kucera, M. & Uthicke, S. Combined effects of warming and ocean acidification on coral reef Foraminifera Marginopora vertebralis and Heterostegina depressa. Coral Reefs 33, 805–818 (2014).ADS 
Article 

Google Scholar 
22.Sinutok, S., Hill, R., Kühl, M., Doblin, M. & Ralph, P. Ocean acidification and warming alter photosynthesis and calcification of the symbiont-bearing foraminifera Marginopora vertebralis. Mar. Biol. 161, 2143–2154 (2014).CAS 
Article 

Google Scholar 
23.ter Kuile, B., Erez, J. & Padan, R. Mechanisms for the uptake of inorganic carbon by two species of symbiont-bearing foraminifera. Mar. Biol. 103, 241–251 (1989).Article 

Google Scholar 
24.Nijweide, P. J., Kawilarang-de Haas, E. W. & Wassenaar, A. M. Alkaline phosphatase and calcification, correlated or not?. Metab. Bone Dis. Relat. Res. 3, 61–66 (1981).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
25.Guo, M. K. & Messer, H. H. A comparison of Ca2+-, Mg2+-ATPase and alkaline phosphatase activities of rat incisor pulp. Calc. Tissue Res. 26, 33–38 (1978).CAS 
Article 

Google Scholar 
26.Vogel, N. & Uthicke, S. Calcification and photobiology in symbiont-bearing benthic foraminifera and responses to a high CO2 environment. J. Exp. Mar. Biol. Ecol. 424–425, 15–24 (2012).Article 
CAS 

Google Scholar 
27.Schiebel, R. & Hemleben, C. Planktic Foraminifers in the Modern Ocean (Springer, 2017).Book 

Google Scholar 
28.Bassinot, F. C., Mélières, F., Gehlen, M., Levi, C. & Labeyrie, L. Crystallinity of foraminifera shells: A proxy to reconstruct past botto m water CO3= changes?. Geochem. Geophys. Geosyst. 5, Q08D10 (2004).Article 

Google Scholar 
29.Broecker, W. & Clark, E. Shell weights from the South Atlantic. Geochem. Geophys. Geosyst. 5, Q03003 (2004).ADS 
Article 

Google Scholar 
30.Beer, C. J., Schiebel, R. & Wilson, P. A. Testing planktic foraminiferal shell weight as a surface water [CO32−] proxy using plankton net samples. Geology 38, 103–106 (2010).ADS 
CAS 
Article 

Google Scholar 
31.Naik, S. S., Naidu, P. D., Govil, P. & Godad, S. Relationship between weights of planktonic foraminifer shell and surface water CO3= concentration during the Holocene and Last Glacial Period. Mar. Geol. 275, 278–282 (2010).ADS 
CAS 
Article 

Google Scholar 
32.Moy, A. D., Howard, W. R., Bray, S. G. & Trull, T. W. Reduced calcification in modern Southern Ocean planktonic foraminifera. Nat. Geosci. 2, 276–280 (2009).ADS 
CAS 
Article 

Google Scholar 
33.Gonzalez-Mora, B., Sierro, F. J. & Flores, J. A. Controls of shell calcification in planktonic foraminifers. Quat. Sci. Rev. 27, 956–961 (2008).ADS 
Article 

Google Scholar 
34.Marr, J. P. et al. Ecological and temperature controls on Mg/Ca ratios of Globigerina bulloides from the southwest Pacific Ocean. Paleoceanography 26, PA2209 (2011).ADS 
Article 

Google Scholar 
35.de Villiers, S. A 425 ka record of foraminiferal shell weight variability in the western Equatorial Pacific. Paleoceanography 18, 1080 (2003).ADS 

Google Scholar 
36.de Villiers, S. Occupation of an ecological niche as the fundamental control on the shell-weight of calcifying planktonic foraminifera. Mar. Biol. 144, 45–50 (2004).Article 

Google Scholar 
37.Reymond, C. E., Lloyd, A., Kline, D. I., Dove, S. G. & Pandolfi, J. M. Decline in growth of foraminifer Marginopora rossi under eutrophication and ocean acidification scenarios. Glob. Change Biol. 19, 291–302 (2013).ADS 
Article 

Google Scholar 
38.Weinkauf, M. F. G., Moller, T., Koch, M. C. & Kucera, M. Calcification intensity in planktic foraminifera reflects ambient conditions irrespective of environmental stress. Biogeosciences 10, 6639–6655 (2013).ADS 
CAS 
Article 

Google Scholar 
39.Doo, S. S. et al. Amelioration of ocean acidification and warming effects through physiological buffering of a macroalgae. Ecol. Evol. 10, 8465–8475 (2020).PubMed 
PubMed Central 
Article 

Google Scholar 
40.Hikami, M. et al. Contrasting calcification responses to ocean acidification between two reef foraminifers harboring different algal symbionts. Geophys. Res. Lett. 38, L19601 (2011).ADS 
Article 
CAS 

Google Scholar 
41.Sanyal, A. et al. Oceanic pH control on the boron isotopic composition of foraminifera: Evidence from culture experiments. Paleoceanography 11, 513–517 (1996).ADS 
Article 

Google Scholar 
42.Anagnostou, E. et al. Changing atmospheric CO2 concentration was the primary driver of early Cenozoic climate. Nature 533, 380–384 (2016).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
43.Foster, G. L. & Rae, J. W. B. Reconstructing ocean pH with boron isotopes in foraminifera. Annu. Rev. Earth Planet. Sci. 44, 207–237 (2016).ADS 
CAS 
Article 

Google Scholar 
44.Kroeker, K. J. et al. Impacts of ocean acidification on marine organisms: Quantifying sensitivities and interaction with warming. Glob. Change Biol. 19, 1884–1896 (2013).ADS 
Article 

Google Scholar 
45.Dove, S. G. et al. Future reef decalcification under a business-as-usual CO2 emission scenario. Proc. Nat. Acad. Sci. 110, 15342–15347 (2013).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
46.Cornwall, C. E. et al. Global declines in coral reef calcium carbonate production under ocean acidification and warming. Proc. Nat. Acad. Sci. 118, 2015265118 (2021).Article 
CAS 

Google Scholar 
47.Langer, M. R., Silk, M. T. & Lipps, J. H. Global ocean carbonate and carbon dioxide production: the role of reef foraminifera. J. Foraminifer. Res 27, 271–277 (1997).Article 

Google Scholar 
48.Pierrot, D., Lewis E. D. & Wallace, D.W. MS EXCEL Program Developed for CO2 System Calculations. ORNL/CDIAC-105a. (Carbon Dioxide Information Analysis Center, Oak Ridge National Laboratory, U.S. Department of Energy, 2006). https://doi.org/10.3334/cdiac/otg.co2sys_xls_cdiac105a.49.Shapiro, S. S. & Wilk, M. B. An analysis of variance test for normality (complete samples). Biometrika 52, 591–611 (1965).MathSciNet 
MATH 
Article 

Google Scholar 
50.Bartlett, M. S. Properties of sufficiency and statistical test. Proc. R. Soc. A 160, 268–282 (1937).ADS 
MATH 

Google Scholar  More

Infection-age structured dynamicsFor the description of the dynamics, we follow the customary infection-age structured approach (for details see for instance Refs.4,10,11,12). Explicitly, we consider the infection-age structured dynamics of the number of individuals ({u}_{I}left(t,tau right)) at time (t) who were infected at time (t-tau) given by$$begin{array}{c}frac{partial }{partial t}{u}_{I}left(t,tau right)+frac{partial }{partial tau }{u}_{I}left(t,tau right)=0end{array}$$
(7)
with boundary condition$$begin{array}{c}{u}_{I}left(t,0right)=jleft(tright).end{array}$$
(8)
Here, (tau) is the time elapsed after infection, referred to as infection age, and (jleft(tright)={int }_{0}^{infty }{k}_{I}(t,tau ){u}_{I}left(t,tau right)dtau) is the incidence, with ({k}_{I}(t,tau )) being the rate of secondary transmissions per single primary case.The solution is obtained through the method of characteristics32 as$$begin{array}{c}{u}_{I}left(t,tau right)=jleft(t-tau right)end{array}$$
(9)
for (tge tau) and ({u}_{I}left(t,tau right)=0) for (t1 for countries and for US locations.The daily death counts (Delta {n}_{W}left(tright)={n}_{W}left(tright)-{n}_{W}left(t-1right)) are considered to contain reporting artifacts if they are negative or if they are unrealistically large. This last condition is defined explicitly as larger than 4 times its previous 14-day average value plus 10 deaths, (Delta {n}_{W}left(tright) >10+4times frac{1}{14}left({n}_{W}left(tright)-{n}_{W}left(t-14right)right)), from a non-sparse reporting schedule with at least 2 consecutive non-zero values before and after the time (t), (Delta {n}_{W}left(tright)ne frac{1}{5}left({n}_{W}left(t+2right)-{n}_{W}left(t-3right)right)).Reporting artifacts identified at time (t) are considered to be the result of previous miscounting. The excess or lack of deaths are imputed proportionally to previous death counts. Explicitly, death counts are updated as$$begin{array}{c}{n}_{W}left(t-1-iright)leftarrow {n}_{W}left(t-1-iright)frac{{n}_{W}{left(t-1right)}_{estimated}}{{n}_{W}left(t-1right)}end{array}$$
(32)
with ({n}_{W}{left(t-1right)}_{estimated}={n}_{W}left(tright)-frac{1}{7}left({n}_{W}left(t-1right)-{n}_{W}left(t-8right)right)) for all (ige 0). In this way, (Delta {n}_{W}left(tright)) is assigned its previous seven-day average value.The expected daily deaths, (Delta {n}_{D}(t)), are obtained through a density estimation multiscale functional, ({f}_{de}), as (Delta {n}_{D}(t)={f}_{de}left(Delta {n}_{W}left(tright)right)), which leads to the estimation of the expected cumulative deaths at time (t) as ({n}_{D}left(tright)={n}_{W}left({t}_{0}right)+{sum }_{s={t}_{0}+1}^{t}Delta {n}_{D}(s)). Specifically,$$begin{array}{c}{f}_{de}left(Delta {n}_{W}left(tright)right)=left(1-{r}_{1}right)d{d}_{0}+{r}_{1}left(left(1-{r}_{2}right)d{d}_{1}+{r}_{2}d{d}_{2}right)end{array}$$
(33)
with$$begin{array}{c}{r}_{1} = {e}^{-0.3d{d}_{1}},end{array}$$
(34)
$$begin{array}{c}{r}_{2} = {e}^{-3d{d}_{2}},end{array}$$
(35)
$$begin{array}{c}d{d}_{0}={ma}_{14}left({ma}_{14}left(Delta {n}_{W}left(tright)right)right),end{array}$$
(36)
$$begin{array}{c}d{d}_{1}={rg}_{12}left({ma}_{14}left(Delta {n}_{W}left(tright)right)right),end{array}$$
(37)
$$begin{array}{c}d{d}_{2}={rg}_{48}left({ma}_{14}left(Delta {n}_{W}left(tright)right)right),end{array}$$
(38)
where ({ma}_{14}left(cdot right)) is a centered moving average with window size of 14 days and ({rg}_{sigma }left(cdot right)) is a centered rolling average through a Gaussian window with standard deviation (sigma). The specific value of the window size has been chosen to mitigate weekly reporting effects. The values of the standard deviations of the Gaussian windows have been selected to achieve a smooth representation of the expected death estimation for each country as shown in the bottom panels of Supplementary Fig. S1.Reporting delaysWe consider an average delay of two days between reporting a death and its occurrence. This value is obtained by comparing the daily death counts reported for Spain1 and their actual values33 from February 15 to March 31, 2020. The values of the root-mean-squared deviation between reported and actual deaths shifted by 0, 1, 2, 3, and 4 days are 77.9, 58.4, 38.5, 58.7, and 88.6 deaths respectively.Infection fatality rate ((IFR))The infection fatality rate is computed assuming homogeneous attack rate as$$begin{array}{c}IFR=frac{1}{{sum }_{a}{g}_{a}}{sum }_{a}{IFR}_{a}{g}_{a} ,end{array}$$
(39)
where ({mathrm{IFR}}_{a}) is the previously estimated (IFR) for the age group (a)5 and ({g}_{a}) is the population in the age group (a) as reported by the United Nations for countries18 and the US Census for states19.Clinical parametersWe obtained the values of the average ({tau }_{G}) and standard deviation ({sigma }_{G}) of the generation time from Ref.13, of the averages of the incubation ({tau }_{I}) and symptom onset-to-death ({tau }_{OD}) times from Refs.5,14, and of the average number of days (Delta {t}_{TP}) of positive testing by an infected individual from Refs.15,17. The average time at which an individual tested positive after infection ({tau }_{TP}) was computed as ({tau }_{TP}={tau }_{I}-2+Delta {t}_{TP}/2), where we have assumed that on average an individual started to test positive 2 days before symptom onset. The average seroconversion time after infection ({tau }_{SP}) was estimated as ({tau }_{I}) plus the 7 days of 50% seroconversion after symptom onset reported in Ref.16.Dynamical constraints implementation with discrete timeWe implemented the dynamical constraints to compute the infectious and infected population as outlined in the main text and as detailed in the previous section of this document, using days as time units. Time delays were rounded to days to assign daily values.The first derivative of the cumulative number of deaths is computed as$$begin{array}{c}frac{d{n}_{D}left(tright)}{dt}=Delta {n}_{D}left(tright),end{array}$$
(40)
with (Delta {n}_{D}left(tright)={n}_{D}left(tright)-{n}_{D}(t-1)).The growth rate was computed explicitly from the discrete time series as the centered 7-day difference$$begin{array}{c}{k}_{G}left(tright)=frac{1}{7}left({mathrm{ln}}left(Delta {n}_{D}left(t+4right)+Delta {n}_{D}left(t+3right)right)-{mathrm{ln}}left(Delta {n}_{D}left(t-3right)+Delta {n}_{D}left(t-4right)right)right).end{array}$$
(41)
The 7-day value was chosen to mitigate reporting artifacts.Confidence and credibility intervalsConfidence intervals associated with death counts were computed using bootstrapping with 10,000 realizations34. These confidence intervals were combined with the credibility intervals of the (IFR) in infectious and infected populations assuming independence and additivity on a logarithmic scale.Fold accuracyThe fold accuracy, ({F}_{A}), is explicitly computed as$$begin{array}{c}{mathrm{log}}{F}_{A}=frac{1}{N}{sum }_{i=1}^{N}left|{mathrm{log}}{x}_{i}^{obs}-{mathrm{log}}{x}_{i}^{est}right|,end{array}$$
(42)
where (left|cdot right|) is the absolute value function, ({x}_{i}^{obs}) is the ({i}^{th}) observation, ({x}_{i}^{est}) is its corresponding estimation, and (N) is the total number of observations.Inference and extrapolationBecause of the delay between infections and deaths, inference for the values of the growth rate and infectious populations ends on December 30, 2020 and for the values of the infected populations ends on December 26, 2020. Extrapolation to the current time (January 21, 2021) is carried out assuming the last growth rate computed.Reproduction numberThe quantities ({R}_{t}) and ({k}_{G}left(tright)) are related to each other through the Euler–Lotka equation, ({R}_{t}^{-1}={int }_{0}^{infty }{f}_{GT}left(tau right){e}^{-{k}_{G}left(tright)tau }dtau ,) which considers (jleft(t-tau right)simeq {e}^{-{k}_{G}left(tright)tau }jleft(tright)) in the renewal equation (jleft(tright)={int }_{0}^{infty }{k}_{I}left(t,tau right)jleft(t-tau right)dtau). Generation times can generally be described through a gamma distribution ({f}_{GT}left(tau right)=frac{{beta }^{alpha }}{Gamma left(alpha right)}{tau }^{alpha -1}{e}^{-beta tau }) with (alpha ={tau }_{G}^{2}/{sigma }_{G}^{2}) and (beta ={tau }_{G}/{sigma }_{G}^{2}), which leads to ({R}_{t}={left(1+{k}_{G}(t)/beta right)}^{alpha }) for ({k}_{G}(t) >-beta) and ({R}_{t}=0) for ({k}_{G}left(tright)le -beta). In the case of the exponentially distributed limit ((alpha simeq 1)) or small values of ({k}_{G}(t)/beta), it simplifies to ({R}_{t}=1+{k}_{G}left(tright){tau }_{G}) for ({k}_{G}left(tright) >-1/{tau }_{G}) and ({R}_{t}=0) for ({k}_{G}left(tright)le -1/{tau }_{G}). Global prevalence data were obtained from multiple data sources35,36,37,38,39,40,41,42, as described in Supplementary Table S1. More

1.Oerke, E.-C. Crop losses to pests. J. Agric. Sci. 144, 31–43 (2006).Article 

Google Scholar 
2.R4P Network. Trends and challenges in pesticide resistance detection. Trends Plant Sci. 21, 834–853 (2016).3.Heap, I. M. The international herbicide-resistant weed database. http://www.weedscience.org/Home.aspx (2021).4.Délye, C., Jasieniuk, M. & Le Corre, V. Deciphering the evolution of herbicide resistance in weeds. Trends Genet. 29, 649–658 (2013).PubMed 
Article 
CAS 

Google Scholar 
5.Gaines, T. A. et al. Mechanisms of evolved herbicide resistance. J. Biol. Chem. 295, 10307–10330 (2020).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
6.Murphy, B. P. & Tranel, P. J. Target-site mutations conferring herbicide resistance. Plants 8, 382 (2019).CAS 
PubMed Central 
Article 
PubMed 

Google Scholar 
7.Beckie, H. J. & Tardif, F. J. Herbicide cross resistance in weeds. Crop Prot. 35, 15–28 (2012).CAS 
Article 

Google Scholar 
8.Han, H. et al. Cytochrome P450 CYP81A10v7 in Lolium rigidum confers metabolic resistance to herbicides across at least five modes of action. Plant J. 105, 79–92 (2021).CAS 
PubMed 
Article 

Google Scholar 
9.Kreiner, J. M. et al. Multiple modes of convergent adaptation in the spread of glyphosate-resistant Amaranthus tuberculatus. Proc. Natl. Acad. Sci. 116, 21076–21084 (2019).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
10.Milani, A. et al. Population structure and evolution of resistance to acetolactate synthase (ALS)-inhibitors in Amaranthus tuberculatus in Italy. Pest Manag. Sci. 77, 2971–2980 (2021).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
11.Clements, D. R. et al. Adaptability of plants invading North American cropland. Agric. Ecosyst. Environ. 104, 379–398 (2004).Article 

Google Scholar 
12.Essl, F. et al. Biological flora of the British Isles: Ambrosia artemisiifolia. J. Ecol. 103, 1069–1098 (2015).Article 

Google Scholar 
13.Cowbrough, M. J., Brown, R. B. & Tardif, F. J. Impact of common ragweed (Ambrosia artemisiifolia) aggregation on economic thresholds in soybean. Weed Sci. 51, 947–954 (2003).CAS 
Article 

Google Scholar 
14.Swinton, S. M., Buhler, D. D., Forcella, F., Gunsolus, J. L. & King, R. P. Estimation of crop yield loss due to interference by multiple weed species. Weed Sci. 42, 103–109 (1994).Article 

Google Scholar 
15.Bassett, I. J. & Crompton, C. W. The biology of Canadian weeds: Ambrosia artemisiifolia L. and A. psilostachya DC. Can. J. Plant Sci. 55, 463–476 (1975).16.Chauvel, B., Dessaint, F., Cardinal-Legrand, C. & Bretagnolle, F. The historical spread of Ambrosia artemisiifolia L. France from herbarium records. J. Biogeogr. 33, 665–673 (2006).Article 

Google Scholar 
17.Sala, C. A., Bulos, M., Altieri, E. & Ramos, M. L. Genetics and breeding of herbicide tolerance in sunflower. Helia 35, 57–69 (2012).Article 

Google Scholar 
18.Yu, Q. & Powles, S. B. Resistance to AHAS inhibitor herbicides: Current understanding. Pest Manag. Sci. 70, 1340–1350 (2014).CAS 
PubMed 
Article 

Google Scholar 
19.Tranel, P. J., Wright, T. R. & Heap, I. M. ALS mutations from resistant weeds. http://www.weedscience.com (2021).20.Patzoldt, W. L., Tranel, P. J., Alexander, A. L. & Schmitzer, P. R. A common ragweed population resistant to cloransulam-methyl. Weed Sci. 49, 485–490 (2001).CAS 
Article 

Google Scholar 
21.Rousonelos, S. L., Lee, R. M., Moreira, M. S., VanGessel, M. J. & Tranel, P. J. Characterization of a common ragweed (Ambrosia artemisiifolia) population resistant to ALS- and PPO-inhibiting herbicides. Weed Sci. 60, 335–344 (2012).CAS 
Article 

Google Scholar 
22.Zheng, D., Patzoldt, W. L. & Tranel, P. J. Association of the W574L ALS substitution with resistance to cloransulam and imazamox in common ragweed (Ambrosia artemisiifolia). Weed Sci. 53, 424–430 (2005).CAS 
Article 

Google Scholar 
23.Van Wely, A. C. et al. Glyphosate and acetolactate synthase inhibitor resistant common ragweed (Ambrosia artemisiifolia L.) in southwestern Ontario. Can. J. Plant Sci. 95, 335–338 (2015)24.Marsan-Pelletier, F., Vanasse, A., Simard, M.-J. & Cuerrier, M.-E. Survey of imazethapyr-resistant common ragweed (Ambrosia artemisiifolia L.) in Quebec. Phytoprotection 99, 36–44 (2019).25.Owen, M. D. & Zelaya, I. A. Herbicide-resistant crops and weed resistance to herbicides. Pest Manag. Sci. 61, 301–311 (2005).CAS 
PubMed 
Article 

Google Scholar 
26.Duke, S. O. & Powles, S. B. Glyphosate: A once-in-a-century herbicide. Pest Manag. Sci. 64, 319–325 (2008).CAS 
PubMed 
Article 

Google Scholar 
27.Barnes, E. R., Knezevic, S. Z., Sikkema, P. H., Lindquist, J. L. & Jhala, A. J. Control of glyphosate-resistant common ragweed (Ambrosia artemisiifolia L.) in glufosinate-resistant soybean [Glycine max (L.) Merr]. Front. Plant Sci. 8, 1455 (2017).28.Tranel, P. J. & Wright, T. R. Resistance of weeds to ALS-inhibiting herbicides: What have we learned?. Weed Sci. 50, 700–712 (2002).CAS 
Article 

Google Scholar 
29.Li, J., Li, M., Gao, X. & Fang, F. A novel amino acid substitution Trp574Arg in acetolactate synthase (ALS) confers broad resistance to ALS-inhibiting herbicides in crabgrass (Digitaria sanguinalis). Pest Manag. Sci. 73, 2538–2543 (2017).CAS 
PubMed 
Article 

Google Scholar 
30.Duggleby, R. G., Pang, S. S., Yu, H. & Guddat, L. W. Systematic characterization of mutations in yeast acetohydroxyacid synthase. Interpretation of herbicide-resistance data. Eur. J. Biochem. 270, 2895–2904 (2003).31.Jung, S.-M. et al. Amino acid residues conferring herbicide resistance in tobacco acetohydroxyacid synthase. Biochem. J. 383, 53–61 (2004).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
32.Owen, M. J., Walsh, M. J., Llewellyn, R. S. & Powles, S. B. Widespread occurrence of multiple herbicide resistance in Western Australian annual ryegrass (Lolium rigidum) populations. Aust. J. Agric. Res. 58, 711–718 (2007).CAS 
Article 

Google Scholar 
33.Owen, M. J., Martinez, N. J. & Powles, S. B. Multiple herbicide-resistant Lolium rigidum (annual ryegrass) now dominates across the Western Australian grain belt. Weed Res. 54, 314–324 (2014).CAS 
Article 

Google Scholar 
34.Délye, C. Nucleotide variability at the acetyl coenzyme A carboxylase gene and the signature of herbicide selection in the grass weed Alopecurus myosuroides (Huds.). Mol. Biol. Evol. 21, 884–892 (2004).35.Délye, C., Clément, J. A. J., Pernin, F., Chauvel, B. & Le Corre, V. High gene flow promotes the genetic homogeneity of arable weed populations at the landscape level. Basic Appl. Ecol. 11, 504–512 (2010).Article 

Google Scholar 
36.Délye, C., Pernin, F. & Scarabel, L. Evolution and diversity of the mechanisms endowing resistance to herbicides inhibiting acetolactate-synthase (ALS) in corn poppy (Papaver rhoeas L.). Plant Sci. 180, 333–342 (2011).37.Sudheesh, M. An analysis of polygenic herbicide resistance evolution and its management based on a population genetics approach. Basic Appl. Ecol. 16, 104–111 (2015).Article 

Google Scholar 
38.Bullock, J. M. Assessing and controlling the spread and the effects of common ragweed in Europe. Report, Contractor: Natural environment research Council UK (2012).39.Yu, Q., Nelson, J. K., Zheng, M. Q., Jackson, J. & Powles, S. B. Molecular characterisation of resistance to ALS-inhibiting herbicides in Hordeum leporinum biotypes. Pest Manag. Sci. 63, 918–927 (2007).CAS 
PubMed 
Article 

Google Scholar 
40.Simard, M.-J., Laforest, M., Soufiane, B., Benoit, D. L. & Tardif, F. Linuron resistant common ragweed (Ambrosia artemisiifolia) populations in Quebec carrot fields: presence and distribution of target-site and non-target site resistant biotypes. Can. J. Plant Sci. 98, 345–352 (2017).
Google Scholar 
41.Ganie, Z., Jugulam, M., Varanasi, V. & Jhala, A. J. Investigating mechanism of glyphosate resistance in a common ragweed (Ambrosia artemisiifolia L.) biotype from Nebraska. Can. J. Plant Sci. (2017). https://doi.org/10.1139/CJPS-2017-0036.42.Duhoux, A., Carrère, S., Duhoux, A. & Délye, C. Transcriptional markers enable identification of rye-grass (Lolium sp.) plants with non-target-site-based resistance to herbicides inhibiting acetolactate-synthase. Plant Sci. 257, 22–36 (2017).43.Gardin, J. A. C., Gouzy, J., Carrère, S. & Délye, C. ALOMYbase, a resource to investigate non-target-site-based resistance to herbicides inhibiting acetolactate-synthase (ALS) in the major grass weed Alopecurus myosuroides (black-grass). BMC Genomics 16, 590 (2015).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
44.Torra, J. et al. Target-site and non-target-site resistance mechanisms confer multiple and cross- resistance to ALS and ACCase inhibiting herbicides in Lolium rigidum from Spain. Front. Plant Sci. 12, 625138 (2021).45.Manley, B. S., Hatzios, K. K. & Wilson, H. P. Absorption, translocation, and metabolism of chlorimuron and nicosulfuron in imidazolinone-resistant and susceptible smooth pigweed (Amaranthus hybridus). Weed Technol. 13, 759–764 (1999).CAS 
Article 

Google Scholar 
46.Jeffers, G. M., O’Donovan, J. T. & Hall, L. M. Wild mustard (Brassica kaber) resistance to ethametsulfuron but not to other herbicides. Weed Technol. 10, 847–850 (1996).CAS 
Article 

Google Scholar 
47.Veldhuis, L. J., Hall, L. M., O’Donovan, J. T., Dyer, W. & Hall, J. C. Metabolism-based resistance of a wild mustard (Sinapis arvensis L.) biotype to ethametsulfuron-methyl. J. Agric. Food Chem. 48, 2986–2990 (2000).48.Scarabel, L., Pernin, F. & Délye, C. Occurrence, genetic control and evolution of non-target-site based resistance to herbicides inhibiting acetolactate synthase (ALS) in the dicot weed Papaver rhoeas. Plant Sci. 238, 158–169 (2015).CAS 
PubMed 
Article 

Google Scholar 
49.Nakka, S., Thompson, C. R., Peterson, D. E. & Jugulam, M. Target site-based and non-target site based resistance to ALS Inhibitors in Palmer Amaranth (Amaranthus palmeri). Weed Sci. 65, 681–689 (2017).Article 

Google Scholar 
50.Meyer, L. et al. New gSSR and EST-SSR markers reveal high genetic diversity in the invasive plant Ambrosia artemisiifolia L. and can be transferred to other invasive Ambrosia species. PLOS ONE 12, e0176197 (2017).51.Van Boheemen, L. A. et al. Multiple introductions, admixture and bridgehead invasion characterize the introduction history of Ambrosia artemisiifolia in Europe and Australia. Mol. Ecol. 26, 5421–5434 (2017).PubMed 
Article 

Google Scholar 
52.Délye, C. et al. Harnessing the power of next-generation sequencing technologies to the purpose of high-throughput pesticide resistance diagnosis. Pest Manag. Sci. 76, 543–552 (2020).PubMed 
Article 
CAS 

Google Scholar 
53.Délye, C., Matéjicek, A. & Gasquez, J. PCR-based detection of resistance to acetyl-CoA carboxylase-inhibiting herbicides in black-grass (Alopecurus myosuroides Huds) and ryegrass (Lolium rigidum Gaud). Pest Manag. Sci. 58, 474–478 (2002).PubMed 
Article 
CAS 

Google Scholar 
54.Duggleby, R. G., McCourt, J. A. & Guddat, L. W. Structure and mechanism of inhibition of plant acetohydroxyacid synthase. Plant Physiol. Biochem. 46, 309–324 (2008).CAS 
PubMed 
Article 

Google Scholar 
55.Leigh, J. W. & Bryant, D. POPART: full-feature software for haplotype network construction. Methods Ecol. Evol. 6, 1110–1116 (2015).Article 

Google Scholar 
56.Neff, M. M., Neff, J. D., Chory, J. & Pepper, A. E. dCAPS, a simple technique for the genetic analysis of single nucleotide polymorphisms: Experimental applications in Arabidopsis thaliana genetics. Plant J. 14, 387–392 (1998).CAS 
PubMed 
Article 

Google Scholar 
57.Délye, C. & Boucansaud, K. A molecular assay for the proactive detection of target site-based resistance to herbicides inhibiting acetolactate synthase in Alopecurus myosuroides. Weed Res. 48, 97–101 (2008).Article 

Google Scholar 
58.Livak, K. J. & Schmittgen, T. D. Analysis of relative gene expression data using real-time quantitative PCR and the 2ddCT method. Methods 25, 402–408 (2001).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
59.Bustin, S. A. et al. The MIQE guidelines: minimum information for publication of quantitative real-time PCR experiments. Clin. Chem. 55, 611–622 (2009).CAS 
PubMed 
Article 

Google Scholar  More

In this study, the HFLUX model was coupled with the SCEM-UA algorithm for analyzing the uncertainties of the model inputs. The specific procedures started with selecting the inputs of the HFLUX model. With the linked HFLUX and SCEM-UA model and implementation of an iteration scheme, the uncertainty of each of the selected inputs was obtained based on the ranges (minimum and maximum values) of the input data/parameters and the Latin hypercube sampling. The simulations were then compared against the observed data to evaluate the performance of the SCEM-UA algorithm. These steps are depicted in Fig. 1.Figure 1Flowchart for the uncertainty analysis.Full size imageRiver water temperatures simulated by the HFLUX modelRiver water temperature affects the water quality and the ecosystem health, and hence control of river water temperature is important to mitigation of its adverse effects1. The HFLUX model was used to simulate the streamflow temperatures at different locations and times. The model is highly flexible in terms of choosing the solution methods for solving the governing equations and selecting the energy budget terms such as shortwave solar radiation, latent heat flux, and sensible heat transfer flux. The model input data include the initial spatial and temporal temperature conditions, stream geometry data, discharge data, and meteorological data8. The water balance and energy balance equations are respectively given by8:$$frac{partial A}{{partial t}} + frac{partial Q}{{partial x}} = mathop qnolimits_{L}$$
(1)
$$frac{{partial left( {Amathop Tnolimits_{w} } right)}}{partial t} + frac{{partial left( {Qmathop Tnolimits_{w} } right)}}{partial x} = mathop qnolimits_{L} mathop Tnolimits_{L} + R$$
(2)
$$R = frac{{Bmathop varphi nolimits_{total} }}{{mathop rho nolimits_{w} mathop Cnolimits_{w} }}$$
(3)
where A is the cross section area of the stream (m2), x is the distance along the stream (m), t is the time (s), Q is the discharge of the stream (m3/s), qL is the lateral inflow per unit stream length (m2/s), Tw is the stream temperature ((^circ C)), TL is the temperature of the lateral inflow ((^circ C)), R is the energy flux (source or sink) per unit stream length ((^circ C) m2/s), B is the width of the stream (m), (mathop varphi nolimits_{total}) is the total energy flux to the stream per surface area (W/m2), (mathop rho nolimits_{w}) is the density of water (kg/m3), and (mathop Cnolimits_{w}) is the specific heat of water (J/kg (^circ C)). Equation (3) is based on a thermal datum of 0 (^circ C) and the impact on the absolute value of the advective heat flux term. In Eq. (2), if qL is negative, the first term on the right-hand side of the equation becomes a loss of qLTw. Also, dispersive heat transport that is omitted in Eq. 2 is negligible when the longitudinal change in water temperature is small in comparison to the temporal changes8.SCEM-UA algorithmThe SCEM-UA algorithm provides posterior distribution functions for the model parameters and input data by generating an initial sample from the parameter space. First, the indicators of n, q, and s that are respectively dimension (the number of investigate inputs), number of complexes (the population to be divided), and population (the number of sample points) are determined for the algorithm. Then, the algorithm searches the sampling points in the feasible space and sorts the points according to the density. The algorithm determines the sequence and complexes based on those points. The sequence is the first q points of the population and complexes are a collection of m points from the population. Note that m = s/q. In the next step, the points of each complex are sorted based on the density, which can be mathematically expressed as20:$$left{ {begin{array}{*{20}c} {mathop alpha nolimits^{k} le T,,,,,,,,,mathop theta nolimits^{t + 1} = Nleft( {mathop theta nolimits^{t} ,,mathop Cnolimits_{n}^{2} mathop Sigma nolimits^{k} } right)} \ {mathop alpha nolimits^{k} > T,,,,,,,,mathop theta nolimits^{t + 1} = Nleft( {mathop mu nolimits^{k} ,,mathop Cnolimits_{n}^{2} mathop Sigma nolimits^{k} } right)} \ end{array} } right.$$
(4)
where k = 1,2,…,q, α is the ratio of the mean posterior density of the m points of complexes to the mean posterior density of the last m generated points of sequences, (theta) is the points of complexes, ({c}_{n}=frac{2.4}{sqrt{n}}) , (T={10}^{6}), (mu) is the mean, and ∑ denotes the covariance. To investigate the new points created by the algorithm, the points of complexes are replaced by20:$$left{ {begin{array}{*{20}l} {Omega ge Zquad replace,best,member,of,mathop Cnolimits^{k} ,with,mathop theta nolimits^{t + 1} } \ {Omega < Zquad mathop theta nolimits^{t + 1} = mathop theta nolimits^{t} ,,,,,,,,,,,,,,,,,,,,,} \ end{array} } right.$$ (5) where (mathop Cnolimits^{k}) is the Kth complex, Z is drawn from the uniform distribution in the range of 0–1, and Ω is calculated by20:$$Omega = frac{{Pleft( {left. {mathop theta nolimits^{t + 1} } right|y} right)}}{{Pleft( {left. {mathop theta nolimits^{t} } right|y} right)}}$$ (6) where (Pleft( {left. {mathop theta nolimits^{t + 1} } right|y} right)) and (Pleft( {left. {mathop theta nolimits^{t} } right|y} right)) are the posterior probability distributions for (mathop theta nolimits^{t + 1}) and (mathop theta nolimits^{t}), respectively. Then, the algorithm examines the following condition for each complex. If it is rejected, the algorithm replaces the worst member ({c}^{k})(the point with the lowest density) with ({theta }^{t+1}) 20.$$mathop Gamma nolimits^{k} le T,,and,,Pleft( {{{mathop theta nolimits^{t + 1} } mathord{left/ {vphantom {{mathop theta nolimits^{t + 1} } y}} right. kern-nulldelimiterspace} y}} right) < ,Pleft( {{{mathop Cnolimits_{m}^{k} } mathord{left/ {vphantom {{mathop Cnolimits_{m}^{k} } y}} right. kern-nulldelimiterspace} y}} right)$$ (7) where ({Gamma }^{k}) is the ratio of the posterior density of the best (the point with the highest density) to the posterior density of the worst member of ({c}^{k}). The last step is to examine (beta) and L. Note that (beta) = 1 and L = m/10. If (beta < L), (beta = beta + 1) and the algorithm returns to sort complex points. Otherwise, the algorithm examines the Gelman and Rubin convergence6, and eventually provides the posterior distribution functions20. The value of the Gelman and Rubin convergence should be less than 1.2. The Gelman and Rubin convergence is examined by:$$R = sqrt {frac{g - 1}{g} + frac{q + 1}{{q.g}}frac{B}{W}}$$ (8) where g is the number of iterations within each sequence, B is the variance between the q sequence means, and W is the average of the q within-sequence variances for the parameter under consideration20.Study AREAMeadowbrook Creek was selected to test the methods proposed in this study8. The creek flows through the City of Syracuse in New York. Thus, this catchment consists of high residential and industrial land covers, which contribute runoff to the main channel. The creek is about 4 km long. A portion of this creek (475 m long) was selected for the modeling for a period of June 13–19, 2012 in this study. The upstream boundary condition in the HFLUX model was set based on the water temperature of the creek observed at the upstream station8. The uncertainty of the model inputs was examined at three selected points as shown in Fig. 2. Note that the input values at these three points had greater relative changes than the changes at other locations, which provided the possibility to improve the evaluation of the algorithm performance. In addition, these three locations had the same sampling of the selected input data. During the simulation period, the streamflow velocity varied within a range of 0.06–0.63 (m/s). The daily temperature changed between 8.9 and 28.2 °C. The relative humidity, used to calculate the total energy flux to the stream per surface area, changed from 36 to 93%. The creek bed mainly consisted of clay, cobbles, sand, and gravel materials. The basic statistics of the data/variables used in the HFLUX model are presented in Table 1. Figure 2 shows the study area, the creek, and the three selected points for analysis.Figure 2Study area and the locations of three evaluation sections (the gray enlarged map shows the State of New York), the map in this Figure is created by Google Earth 7.0.2.8415 (https://google.com/earth/versions).Full size imageTable 1 Basic statistics of the data/variables used in the HFLUX model.Full size tableEthical approvalAll authors accept all ethical approvals.Consent to participateAll authors consent to participate.Consent to publishAll authors consent to publish. More