Philippa Kaur

Data sourcesData for the trophic parameters and data categories listed in Tables 1 and 2 were gathered from peer-reviewed scientific publications, grey literature (e.g., agency reports, theses, and dissertations) and unpublished data by the authors of this paper. Data compilation on stomach contents, stable isotopes, FATM, and trophic positions, focussed on mesopelagic organisms, their potential prey and predators. For major and trace elements, energy density and estimates of diet proportions, our search concentrated on mesopelagic taxa. Nevertheless, we also gathered information from small or intermediate-sized epi-, bathy- or benthopelagic species found in the compiled data sources. These species were included because they play key roles in most marine ecosystems, both as important consumers of phytoplankton and zooplankton, and prey for many top predators, and can represent alternative energy pathways to mesopelagic organisms. However, we stress that the data coverage for these species in the current version of the database is very incomplete. Our main interest was on data from the central and eastern North Atlantic, and the Mediterranean, corresponding to the study regions of the SUMMER project. When we could not find suitable data within this region, we extended the geographic scope of our literature search to the western North Atlantic. We did not search for datasets in open access repositories since those data can be easily accessed and extracted. However, some of the data provided by the authors of this paper have been previously deposited in PANGAEA.DNA sequencing-based methods, such as metabarcoding and direct shotgun sequencing, are emerging as promising tools in dietary analyses due to the high resolution in taxonomic identification of many prey simultaneously, and the potential to provide quantitative diet estimates from relative read abundance29. However, recent studies have shown that various methodological and biological factors can break the correlation between the number and abundance of ingested prey and the prey DNA present in the sample, and lead to biased estimates of taxonomic diversity and composition of diet29,30. Given the uncertainties remaining in the interpretation of DNA sequencing-based diet data, we decided not to include these data in MesopTroph until additional research demonstrates that these techniques can be confidently applied for quantitative diet assessment.We identified available data sources in the literature through systematic searches on Web of Science, Google Scholar, ResearchGate, and the Google search engine. We used multiple combinations of terms related to specific data categories (Table 3), in conjunction with the common or scientific taxon names (from genus to order), and the ocean basin. For example, the search for stomach content data of fishes belonging to the family Myctophidae was undertaken using the following terms: “stomach content” OR “gut content” OR “prey composition” OR “diet composition”, AND “mesopelagic fish” OR “myctophid” OR “Myctophiformes” OR “Myctophidae”, AND “Atlantic” or “Mediterranean”. For the mesopelagic and predator species known to be numerically abundant in the SUMMER study regions, we performed a second literature search using the common or scientific name of the species, along with the terms “diet”, “feeding habits”, “trophic ecology”, “trophic markers”, or “food web”. We also examined the literature cited within each collected publication to locate additional data sources.Table 3 Terms used in the literature search for each data category.Full size tableWe next screened the full text of the compiled studies and retained data sources that: (1) were collected within the region of interest, (2) reported quantitative data for the trophic parameters of interest, (3) reported the number of samples for pooled or aggregated data, and (4) provided sufficient details on the methodology to enable a quality check. In the case of stable isotope data, we only included data sources reporting both δ13C and δ15N measurements.Data extraction, cleaning, and formattingWe created a template table for each data category in Microsoft Excel to assemble all datasets into a single file, and to facilitate cleaning and standardization of data records. We added a large number of metadata fields to the tables to annotate details about the sampling (e.g., location, date, methods), sampled specimen(s) (e.g., taxonomy, number and size of individuals, number of replicates, tissue analysed), and data source (e.g., full reference, DOI) for every record.Data contributors formatted and incorporated their datasets directly into the tables. For published sources, the data and associated metadata were extracted manually or digitized from the article text, tables, or supplementary material into the tables. Extraneous or hidden characters, and values such as “NA” (Not Available) or “ND” (Not Determined), were deleted from the parameter and metadata fields. Measurements of trophic parameters were standardized to the same units (see Tables 1 and 2). Parameter values that were clearly incorrect (e.g., δ15N  > 20, or the frequency of occurrence of a prey higher than the number of stomachs sampled) were corrected by searching for the value within the data source. When values could not be corrected, we deleted that data record.When available, we extracted information at the individual level. However, most studies reported data obtained from pooled samples of the same species. In some cases (e.g., small specimens such as planktonic organisms), a minimum and maximum number of individuals in the sample was provided instead of the actual number of individuals sampled. We added two columns to the tables presenting the minimum and maximum number of individuals in the sample. By filtering the column “Ind No (maximum per sample)” for values >1, users can easily identify records with aggregated data and differentiate them from records where information was drawn from a single individual (i.e., where “Ind No (maximum per sample)” =1). In addition, the tables Stomach contents and Estimates of diet proportions include a field “Sample ID” with a unique identifier of the sample. If data are reported at the individual level (i.e., “Ind No (maximum per sample)” =1) then Sample ID is the individual animal ID. If the data are from a group of individuals (i.e., “Ind No (maximum per sample)” >1), then Sample ID identifies that group.We standardized the taxonomic classification and nomenclature of fishes and elasmobranchs following the Eschmeyer’s Catalog of Fishes (http://researcharchive.calacademy.org/research/ichthyology/catalog/fishcatmain.asp)31,32. For the remaining taxa, we used the World Register of Marine Species (http://www.marinespecies.org/)33. Unaccepted or alternate taxon names were replaced by the most up-to-date valid name. When the identification of a taxon was uncertain, the taxonomic level of identification was decreased to a satisfactory level. For example, prey reported as “Cephalopods” were changed to “Cephalopoda”, “Sepiolids” to “Sepiolidae”, and “Myctophum punctatum?” to the genus “Myctophum”.Stomach contentsStomach contents analysis is a standard dietary assessment method that potentially enables quantifying diet components with high taxonomic resolution34. Three parameters are typically used to describe diet composition from stomach contents: the number of individuals of a prey type as a proportion of the total number of prey items (%N), the proportion of a prey item by weight or volume (%W), and the proportion of stomachs containing a particular prey item (i.e., percent frequency of occurrence, %F)35. When available, we collected data on the three parameters, as well as on the absolute number, weight, and frequency of occurrence of each prey type in the stomachs of each sampled individual or group of individuals. If stated in the data source, we indicate if prey weights were directly measured or reconstructed from hard remains (fish otoliths and vertebrae, cephalopod beaks), and if they represent dry or wet weight. Some datasets contained records of prey items without corresponding weights or numbers. As a result, the cumulative percent of all prey items did not sum to 100%. This occurred in 11 data records for the cumulative %W, and nine for the cumulative %N. While we checked the accuracy of percentage values and adjusted rounding errors, we did not attempt to fill in missing values nor did we remove records with missing values. When prey values were reported by an upper bound (e.g., “ More

Hart, M. B. et al. The search for the origin of the planktic foraminifera. J. Geol. Soc. Lond. 160, 341–343 (2003).Article 

Google Scholar 
Aze, T. et al. A phylogeny of Cenozoic macroperforate planktonic foraminifera from fossil data. Biol. Rev. 86, 900–927 (2011).Article 
PubMed 

Google Scholar 
Gradstein, F., Waskowska, A. & Glinskikh, L. The first 40 million years of planktonic foraminifera. Geosci 11, 1–25 (2021).Article 

Google Scholar 
Ujiié, Y., Kimoto, K. & Pawlowski, J. Molecular evidence for an independent origin of modern triserial planktonic foraminifera from benthic ancestors. Mar. Micropaleontol. 69, 334–340 (2008).Article 
ADS 

Google Scholar 
Darling, K. F. et al. Surviving mass extinction by bridging the benthic/planktic divide. Proc. Natl Acad. Sci. USA 106, 12629–33 (2009).Article 
ADS 
PubMed 
PubMed Central 
CAS 

Google Scholar 
Kucera, M. et al. Caught in the act: anatomy of an ongoing benthic–planktonic transition in a marine protist. J. Plankton Res. 39, 436–449 (2017).
Google Scholar 
Ezard, T. H. G., Aze, T., Pearson, P. N. & Purvis, A. Interplay between changing climate and species’ ecology drives macroevolutionary dynamics. Science 332, 349–352 (2011).Article 
ADS 
PubMed 
CAS 

Google Scholar 
Lowery, C. M., Bown, P. R., Fraass, A. J. & Hull, P. M. Ecological response of plankton to environmental change: thresholds for extinction. Annu. Rev. Earth Planet. Sci. 48, 403–429 (2020).Article 
ADS 
CAS 

Google Scholar 
Pawlowski, J., Holzmann, M. & Tyszka, J. New supraordinal classification of foraminifera: molecules meet morphology. Mar. Micropaleontol. 100, 1–10 (2013).Article 
ADS 

Google Scholar 
Lecroq, B. et al. Ultra-deep sequencing of foraminiferal microbarcodes unveils hidden richness of early monothalamous lineages in deep-sea sediments. Proc. Natl Acad. Sci. USA 108, 13177–13182 (2011).Article 
ADS 
PubMed 
PubMed Central 
CAS 

Google Scholar 
Pawlowski, J. et al. The evolution of early foraminifera. Proc. Natl Acad. Sci. USA 100, 11494–8 (2003).Article 
ADS 
PubMed 
PubMed Central 
CAS 

Google Scholar 
Vachard, D. Macroevolution and biostratigraphy of paleozoic foraminifers. in Stratigraphy and Timescales (Ed. Montenari, M.) Vol. 1, 257–323 (Academic Press, 2016).Ibarbalz, F. M. et al. Global trends in marine plankton diversity across kingdoms of life. Cell 179, 1084–1097.e21 (2019).Article 
PubMed 
PubMed Central 
CAS 

Google Scholar 
Guillou, L. et al. The Protist Ribosomal Reference database (PR2): a catalog of unicellular eukaryote Small Sub-Unit rRNA sequences with curated taxonomy. Nucleic Acids Res. 41, D597–D604 (2013).Article 
PubMed 
CAS 

Google Scholar 
Holzmann, M. & Pawlowski, J. An updated classification of rotaliid foraminifera based on ribosomal DNA phylogeny. Mar. Micropaleontol. 132, 18–34 (2017).Article 
ADS 

Google Scholar 
John, A. W. G. The regular occurrence of Reophax Scottie Chaster, a benthic foraminiferan, in plankton samples from the North Sea. J. Micropalaeontol. 6, 61–63 (1987).Article 

Google Scholar 
Kucera, M. et al. Caught in the act: anatomy of an ongoing benthic-planktonic transition in a marine protist. J. Plankton Res. 39, 436–449 (2017).Darling, K. F., Wade, C. M., Kroon, D. & Brown, A. J. L. Planktic foraminiferal molecular evolution and their polyphyletic origins from benthic taxa. Mar. Micropaleontol. 30, 251–266 (1997).Article 
ADS 

Google Scholar 
Church, S. H., Ryan, J. F. & Dunn, C. W. Automation and evaluation of the SOWH test with SOWHAT. Syst. Biol. 64, 1048–1058 (2015).Article 
PubMed 
PubMed Central 

Google Scholar 
Shimodaira, H. An approximately unbiased test of phylogenetic tree selection. Syst. Biol. 51, 492–508 (2002).Article 
PubMed 

Google Scholar 
Pawlowski, J. et al. Extreme differences in rates of molecular evolution of foraminifera revealed by comparison of ribosomal DNA sequences and the fossil record. Mol. Biol. Evol. 14, 498–505 (1997).Article 
PubMed 
CAS 

Google Scholar 
Peijnenburg, K. T. C. A. et al. The origin and diversification of pteropods precede past perturbations in the Earth’s carbon cycle. Proc. Natl Acad. Sci. USA 117, 25609–25617 (2020).Article 
ADS 
PubMed 
PubMed Central 
CAS 

Google Scholar 
O’Brien, C. L. et al. Cretaceous sea-surface temperature evolution: constraints from TEX86 and planktonic foraminiferal oxygen isotopes. Earth-Sci. Rev. 172, 224–247 (2017).Article 
ADS 

Google Scholar 
Olsson, R. K., Berggren, W. A., Hemleben, C. & Huber, B. T. Atlas of Paleocene planktonic foraminifera. Smithson. Contrib. Paleobiol. 1–252 https://doi.org/10.5479/si.00810266.85.1 (1999).Arenillas, I. & Arz, J. A. Benthic origin and earliest evolution of the first planktonic foraminifera after the Cretaceous/Palaeogene boundary mass extinction. Hist. Biol. 29, 25–42 (2017).Article 

Google Scholar 
Huber, B. T., Petrizzo, M. R. & MacLeod, K. G. Planktonic foraminiferal endemism at southern high latitudes following the terminal cretaceous extinction. J. Foraminifer. Res. 50, 382–402 (2020).Article 

Google Scholar 
Arenillas, I., Arz, J. A. & Gilabert, V. An updated suprageneric classification of planktic foraminifera after growing evidence of multiple benthic-planktic transitions. Spanish J. Palaeontol. https://doi.org/10.7203/sjp.22189 (2022).Culver, S. J. Benthic foraminifera across the Cretaceous–Tertiary (K–T) boundary: a review. Mar. Micropaleontol. 47, 177–226 (2003).Article 
ADS 

Google Scholar 
Widmark, J. G. V. & Malmgren, B. A. Benthic foraminiferal changes across the Cretaceous/Tertiary boundary in the deep sea; DSDP sites 525, 527, and 465. J. Foraminifer. Res. 22, 81–113 (1992).Article 

Google Scholar 
Rigaud, S., Martini, R. & Vachard, D. Early evolution and new classification of the order Robertinida (foraminifera). J. Foraminifer. Res. 45, 3–28 (2015).Article 

Google Scholar 
Rigaud, S., Granier, B. & Masse, J. P. Aragonitic foraminifers: an unsuspected wall diversity. J. Syst. Palaeontol. 19, 461–488 (2021).Article 

Google Scholar 
Hull, P. M. et al. On impact and volcanism across the Cretaceous-Paleogene boundary. Science 367, 266–272 (2020).Article 
ADS 
PubMed 
CAS 

Google Scholar 
Morard, R. et al. PFR2: a curated database of planktonic foraminifera 18S ribosomal DNA as a resource for studies of plankton ecology, biogeography and evolution. Mol. Ecol. Resour. 15, 1472–1485 (2015).Article 
PubMed 
CAS 

Google Scholar 
Morard, R. et al. Genetic and morphological divergence in the warm-water planktonic foraminifera genus Globigerinoides. PLoS ONE 14, 1–30 (2019).Article 

Google Scholar 
Morard, R., Vollmar, N. M., Greco, M. & Kucera, M. Unassigned diversity of planktonic foraminifera from environmental sequencing revealed as known but neglected species. PLoS ONE 14, e0213936 (2019).Article 
PubMed 
PubMed Central 
CAS 

Google Scholar 
Camacho, C. et al. BLAST+: Architecture and applications. BMC Bioinforma. 10, 1–9 (2009).Article 

Google Scholar 
R Development Core Team. R: A Language and Environment for Statistical Computing (R Foundation for Statistical Computing, 2014).Liaw, A. & Wiener, M. Classification and Regression by randomForest. R. N. 2, 18–22 (2002).
Google Scholar 
Lang, M. et al. mlr3: a modern object-oriented machine learning framework in R. J. Open Source Softw. 4, 1903 (2019).Article 
ADS 

Google Scholar 
Katoh, K. & Standley, D. M. MAFFT multiple sequence alignment software version 7: improvements in performance and usability. Mol. Biol. Evol. 30, 772–780 (2013).Article 
PubMed 
PubMed Central 
CAS 

Google Scholar 
Darriba, D. et al. ModelTest-NG: a new and scalable tool for the selection of DNA and protein evolutionary models. Mol. Biol. Evol. 37, 291–294 (2020).Article 
MathSciNet 
PubMed 
CAS 

Google Scholar 
Kozlov, A. M. et al. RAxML-NG: a fast, scalable and user-friendly tool for maximum likelihood phylogenetic inference. Bioinformatics 35, 4453–4455 (2019).Article 
PubMed 
PubMed Central 
CAS 

Google Scholar 
Barbera, P. et al. EPA-ng: massively parallel evolutionary placement of genetic sequences. Syst. Biol. 68, 365–369 (2019).Article 
MathSciNet 
PubMed 

Google Scholar 
Letunic, I. & Bork, P. Interactive tree of life (iTOL) v5: an online tool for phylogenetic tree display and annotation. Nucleic Acids Res. 49, 293–296 (2021).Article 

Google Scholar 
Löytynoja, A. & Goldman, N. WebPRANK: a phylogeny-aware multiple sequence aligner with interactive alignment browser. BMC Bioinform. 11, 1–7 (2010).Ronquist, F. et al. MrBayes 3. 2: efficient Bayesian phylogenetic inference and model choice across a large model space. Syst. Biol. 61, 539–542 (2012).Article 
PubMed 
PubMed Central 

Google Scholar 
Minh, B. Q. et al. IQ-TREE 2: new models and efficient methods for phylogenetic inference in the genomic era. Mol. Biol. Evol. 37, 1530–1534 (2020).Article 
PubMed 
PubMed Central 
CAS 

Google Scholar 
Dos Reis, M., Donoghue, P. C. J. & Yang, Z. Bayesian molecular clock dating of species divergences in the genomics era. Nat. Rev. Genet. 17, 71–80 (2016).Article 
PubMed 

Google Scholar 
Song, H., Tong, J. & Chen, Z. Q. Evolutionary dynamics of the Permian-Triassic foraminifer size: Evidence for Lilliput effect in the end-Permian mass extinction and its aftermath. Palaeogeogr. Palaeoclimatol. Palaeoecol. 308, 98–110 (2011).Article 

Google Scholar 
Copestake, P. & Johnson, B. Lower Jurassic Foraminifera from the Llanbedr (Mochras Farm) Borehole, North Wales, UK. Monogr. Palaeontogr. Soc. 167, 1–403 (2013).Article 

Google Scholar 
Rigaud, S. & Blau, J. New Robertinid Foraminifers from the Early Jurassic of Adnet, Austria and Their Evolutionary Importance. Acta Palaeontol. Pol. 61, 721–734 (2016).Article 

Google Scholar 
Boudagher-fadel, M. K. Evolution and Geological Significance of Larger Benthic Foraminifera. Evolution and Geological Significance of Larger Benthic Foraminifera (UCL Press, 2018).Piuz, A. & Meister, C. Cenomanian rotaliids (Foraminiferida) from Oman and Morocco. Swiss J. Palaeontol. 132, 81–97 (2013).Article 

Google Scholar 
Kucera, M. & Schönfeld, J. The origin of modern oceanic foraminiferal faunas and Neogene climate change. in Deep-Time Perspectives on Climate Change: Marrying the Signal from Computer Models and Biological Proxies. (ed. The Micropalaeontological Society, S. P.) 409–425 (The Geological Society, 2007).Drummond, A. J. & Suchard, M. A. Bayesian random local clocks, or one rate to rule them all. BMC Biol. 8, 114 (2010).Article 
PubMed 
PubMed Central 

Google Scholar 
Rambaut, A. FigTree version 1.3.1. http://tree.bio.ed.ac.uk (2009).Groussin, M., Pawlowski, J. & Yang, Z. Bayesian relaxed clock estimation of divergence times in foraminifera. Mol. Phylogenet. Evol. 61, 157–166 (2011).Article 
PubMed 

Google Scholar 
Loeblich Jr, A. R. & Tappan, H. Foraminiferal Genera and Their Classification (Springer, 1988). More

DNA sample collectionTo analyze the population structure of wild medaka populations, we selected samples from the DNA collection of Takehana et al.29, deposited in University of Shizuoka. The original DNA collection had been made throughout 1980s and 2000s. The selected samples covered the major mitotypes and contained more than three individuals of each population (Table S11, Fig. 3), which were collected from three collection sites for O. sakaizumii and 12 collection sites for O. latipes. We also examined several artificial strains: HNI and Hd-rR, which are inbred strains derived from O. sakaizumii and O. latipes, respectively, and four Himedaka individuals from commercial stock (Uruma city, Okinawa Prefecture, Japan).In addition, samples were newly collected at Kunigami Village, Okinawa Prefecture. Live fish were anesthetized with MS-222 (aminobenzene methanesulfonate, FUJIFILM Wako Pure Chemical Corporation, Osaka, Japan) and then fixed in 99% ethanol. Genomic DNA was extracted using a DNeasy kit (Qiagen Inc., Hilden, Germany) from ethanol-fixed pectoral fin samples according to the manufacturer’s protocol. The DNA concentration was measured using a spectrophotometer (Nanodrop 1000, Thermo Fisher Scientific, Waltham, Massachusetts, USA), and the DNA was diluted with PCR-grade water to a concentration of c.a. 10 ng/µl (UltraPure™ DNase/RNase-Free Distilled Water, Thermo Fisher Scientific).Ethic statementAll methods were carried out in accordance with the Regulation for Animal Experiments at University of the Ryukyus for handling live fish. All experiments were approved by the Animal Care Ethics Committee of University of the Ryukyus (R2019035). All experimental methods are reported in accordance with ARRIVE guidelines.PCR primer designThe following steps were used to select primers for MAAS (Fig. 1). (1) All possible 10-mer sequence combinations (i.e., 410 = 1,048,576 sequences) were generated in silico. (2) The sequences containing simple sequence repeats, some of which had been used in the MIG-seq method17, were excluded. (3) Sequences containing a functional motif, such as a transcription factor-binding site, were also excluded because they may not be suitable for examining neutral genetic markers. We obtained a catalog of motifs from the JASPAR CORE40 (http://jaspar.genereg.net). (4) To avoid taxon-dependency in primer performance, we used information about the k-mer (k = 10) frequency of reference genomes from multiple phyla. Sequences that showed marked differences in frequency among taxa were excluded. The frequencies of each 10-mer sequence in the reference genomes of 17 species belonging to 12 phyla of metazoa were counted (Table S12) using the “oligonucleotideFrequency” function in the “Biostrings” package ver. 2.441. In each of these taxa, the frequencies of sequences were stratified into three grades ( 103). We then selected the sequences that showed the same grade in more than 80% (14/17) of the species. (5) To avoid synthesizing primer dimers, self-complementary sequences were excluded, taking Illumina adapter sequences (5′-CGCTCTTCCGATCT-3′ and 5′-TGCTCTTCCGATCT-3′) into account. Self-complementation of two bases at the 3′-end or every three continuous bases in primer sequences was then evaluated using a custom script in R ver. 3.5.0 (R Development Core Team, http://cran.r-project.org). Based on the selected 10-mer sequences (i.e., 129 sequences, Fig. 1), 7-mer primer sequences were designed by removing the 3 bases at the 3′ end. Finally, we selected 24 candidate sequences for both 10-mer and 7-mer primers for the subsequent step (Table S1).The primer sequence consisted of three parts17: partial sequence of the Illimina adapter, 7 N bases, and a short priming sequence, e.g., 5′-CGCTCTTCCGATCTNNNNNNNGTCGCCC-3′. PCR amplification was performed using the candidate primers using the first PCR protocol described below (Table S1). Banding patterns were observed by electrophoresis on 1% agarose gels (agarose S; TaKaRa, Japan). Of the candidate primers, we selected four 7-mer primers and four 10-mer primers that each gave a smeared banding pattern with amplification products ranging from 500 to 2000 bp, indicating uniform amplification of multiple target sequences (Table S1).Library construction and sequencingThe library was constructed by a two-step PCR approach using a modification of a MIG-seq protocol14. In the first PCR step, multiple regions of genomic DNA were amplified using a cocktail of primers with a Multiplex PCR Assay Kit Ver.2 (TaKaRa) (Table 1). The volume of the PCR reaction mixture was 7 μl, containing 1 μl of template DNA, 2 μM of each PCR primer, 3.5 μl of 2 × Multiplex PCR Buffer, and 0.035 μl of Multiplex PCR Enzyme Mix. PCR was performed under the following conditions: denaturation at 94 °C for 1 min; 25 cycles of 94 °C for 30 s, 38 °C for 1 min, and 72 °C for 1 min, followed by a final extension step at 72 °C for 10 min.The primers in the second PCR step contained the Illumina sequencing adapter and an index sequence to identify each sample. Following the Truseq indexes, we used the combinations of eight forward indexes (i5) and 12 reverse indexes (i7), which resulted in a total of 96 combinations. To be used as a template for the second PCR, the first PCR product from each sample was diluted 50 times with PCR-grade water. The second PCR was performed in a 15-μl reaction mixture containing, 3 μl of diluted first PCR product, 3 μl of 5 × PrimeSTAR GXL Buffer, 200 μM of each dNTP, 0.2 μM of forward index primer and reverse index primer, 0.375 U of PrimeSTAR GXL DNA Polymerase (TaKaRa). The PCR conditions were as follows: 12 cycles at 98 °C for 10 s, 54 °C for 15 s, and 68 °C for 30 s.The second PCR product of each sample was pooled by equal volume and size-selected from 600 to 1000 bp using solid phase reversible immobilization (SPRI) select beads (Beckman Coulter Inc, Brea, California, USA) according to the manufacturer’s protocol. The DNA concentration of the pooled library was measured using a Qubit fluorometer (Thermo Fisher Scientific). We sequenced the libraries using two NGS platforms, MiSeq (Illumina, MiSeq Reagent Kit v2 Micro, Paired-End (PE), 150 bp) and HiSeq X (Illumina, PE, 150 bp). Sequencing using the HiSeq X platform was performed by Macrogen Japan (Tokyo, Japan).To compare primer performance, the DNA libraries constructed using the 7-mer and 10-mer primers for one individual were sequenced using MiSeq. Then, a 7-mer primer cocktail containing four sets of mixed primers was used for the subsequent analyses (Table 1). We also constructed DNA libraries using 7-mer and MIG-seq primer cocktails for three individuals and sequenced them using the HiSeq X platform. Finally, we constructed DNA libraries using 7-mer primer cocktails for 67 wild individuals and six artificial strain individuals for population genetics analyses (Table S11, Fig. 3).Mapping and SNV callingGenotyping was conducted using the following BWA-GATK best-practices pipeline for each sample42. Primer sequences were removed using cutadapt with the –b option selected43. The Illumina adapter sequences were also removed and quality filtering was performed using fastp ver. 0.20.0 with the “–detect_adapter_for_pe, –cut_front” option selected44. The remaining reads were mapped on the reference genome of medaka, Hd-rR strain, GCA_002234675.1; ASM223467v127 using Burrows-Wheeler Alignment tool, BWA mem ver. 0.7.1745. After mapping, output files were converted to Binary Alignment/Map (BAM) format using SAMtools ver. 1.746. SNVs and InDels in the sample were determined following the best practice guidelines set out in the Genome Analysis Tool Kit (GATK ver. 3.8.0)42. We then filtered out SNVs and InDels based on the following criteria: “QD  60.0 || MQ  More

Prosser JI, Hink L, Gubry-Rangin C, Nicol GW. Nitrous oxide production by ammonia oxidizers: Physiological diversity, niche differentiation and potential mitigation strategies. Glob Chang Biol. 2020;26:103–18.Article 
PubMed 

Google Scholar 
Huang L, Chakrabarti S, Cooper J, Perez A, John SM, Daroub SH, et al. Ammonia-oxidizing archaea are integral to nitrogen cycling in a highly fertile agricultural soil. ISME Commun. 2021;1:19.Article 

Google Scholar 
Hink L, Gubry-Rangin C, Nicol GW, Prosser JI. The consequences of niche and physiological differentiation of archaeal and bacterial ammonia oxidisers for nitrous oxide emissions. ISME J. 2018;12:1084–93.Article 
PubMed 
PubMed Central 
CAS 

Google Scholar 
Li Y, Chapman SJ, Nicol GW, Yao H. Nitrification and nitrifiers in acidic soils. Soil Biol Biochem. 2018;116:290–301.Article 
CAS 

Google Scholar 
Ahlgren NA, Fuchsman CA, Rocap G, Fuhrman JA. Discovery of several novel, widespread, and ecologically distinct marine Thaumarchaeota viruses that encode amoC nitrification genes. ISME J 2019;13:618–31.Article 
PubMed 
CAS 

Google Scholar 
Kim J-G, Kim S-J, Cvirkaite-Krupovic V, Yu W-J, Gwak J-H, López-Pérez M, et al. Spindle-shaped viruses infect marine ammonia-oxidizing thaumarchaea. Proc Natl Acad Sci USA. 2019;116:15645–50.Article 
PubMed 
PubMed Central 
CAS 

Google Scholar 
Emerson JB. Soil Viruses: A New Hope. mSystems 2019;4:e00120–19.Article 
PubMed 
PubMed Central 
CAS 

Google Scholar 
Santos-Medellín C, Estera-Molina K, Yuan M, Pett-Ridge J, Firestone MK, Emerson JB. Spatial turnover of soil viral populations and genotypes overlain by cohesive responses to moisture in grasslands. Proc Natl Acad Sci USA. 2022;119:e2209132119.Article 
PubMed 
PubMed Central 

Google Scholar 
Wu R, Davison MR, Gao Y, Nicora CD, Mcdermott JE, Burnum-Johnson KE, et al. Moisture modulates soil reservoirs of active DNA and RNA viruses. Commun Biol. 2021;4:992.Article 
PubMed 
PubMed Central 
CAS 

Google Scholar 
Braga LPP, Spor A, Kot W, Breuil M-C, Hansen LH, Setubal JC, et al. Impact of phages on soil bacterial communities and nitrogen availability under different assembly scenarios. Microbiome 2020;8:52.Article 
PubMed 
PubMed Central 
CAS 

Google Scholar 
Albright MBN, Gallegos-Graves LV, Feeser KL, Montoya K, Emerson JB, Shakya M, et al. Experimental evidence for the impact of soil viruses on carbon cycling during surface plant litter decomposition. ISME Commun. 2022;2:24.Article 

Google Scholar 
Starr EP, Shi S, Blazewicz SJ, Koch BJ, Probst AJ, Hungate BA, et al. Stable-isotope-informed, genome-resolved metagenomics uncovers potential cross-kingdom interactions in rhizosphere soil. mSphere 2021;6:e0008521.Article 
PubMed 

Google Scholar 
Trubl G, Kimbrel JA, Liquet-Gonzalez J, Nuccio EE, Weber PK, Pett-Ridge J, et al. Active virus-host interactions at sub-freezing temperatures in Arctic peat soil. Microbiome 2021;9:208.Article 
PubMed 
PubMed Central 
CAS 

Google Scholar 
Lee S, Sieradzki ET, Nicolas AM, Walker RL, Firestone MK, Hazard C, et al. Methane-derived carbon flows into host–virus networks at different trophic levels in soil. Proc Natl Acad Sci USA. 2021;118:e2105124118.Article 
PubMed 
PubMed Central 
CAS 

Google Scholar 
Nicol GW, Leininger S, Schleper C, Prosser JI. The influence of soil pH on the diversity, abundance and transcriptional activity of ammonia oxidizing archaea and bacteria. Environ Microbiol. 2008;10:2966–78.Article 
PubMed 
CAS 

Google Scholar 
Chaumeil P-A, Mussig AJ, Hugenholtz P, Parks DH. GTDB-Tk: a toolkit to classify genomes with the Genome Taxonomy Database. Bioinformatics 2020;36:1925–7.CAS 

Google Scholar 
Alves RJE, Minh BQ, Urich T, von Haeseler A, Schleper C. Unifying the global phylogeny and environmental distribution of ammonia-oxidising archaea based on amoA genes. Nat Commun. 2018;9:1517.Article 
PubMed 
PubMed Central 

Google Scholar 
Cardinale DJ, Duffy S. Single-stranded genomic architecture constrains optimal codon usage. Bacteriophage 2011;1:219–24.Article 
PubMed 
PubMed Central 

Google Scholar 
Lee S, Sorensen JW, Walker RL, Emerson JB, Nicol GW, Hazard C. Soil pH influences the structure of virus communities at local and global scales. Soil Biol Biochem. 2022;166:108569.Article 
CAS 

Google Scholar 
Jang HB, Bolduc B, Zablocki O, Kuhn JH, Roux S, Adriaenssens EM, et al. Taxonomic assignment of uncultivated prokaryotic virus genomes is enabled by gene-sharing networks. Nat Biotechnol. 2019;37:632–9.Article 

Google Scholar 
Nishimura Y, Yoshida T, Kuronishi M, Uehara H, Ogata H, Goto S. ViPTree: the viral proteomic tree server. Bioinformatics 2017;33:2379–80.Article 
PubMed 
CAS 

Google Scholar 
Kerou M, Offre P, Valledor L, Abby SS, Melcher M, Nagler M, et al. Proteomics and comparative genomics of Nitrososphaera viennensis reveal the core genome and adaptations of archaeal ammonia oxidizers. Proc Natl Acad Sci Usa 2016;113:7937–46.Article 

Google Scholar 
Reyes C, Hodgskiss LH, Kerou M, Pribasnig T, Abby SS, Bayer B, et al. Genome wide transcriptomic analysis of the soil ammonia oxidizing archaeon Nitrososphaera viennensis upon exposure to copper limitation. ISME J 2020;14:2659–74.Article 
PubMed 
PubMed Central 
CAS 

Google Scholar 
Sieradzki ET, Greenlon A, Nicolas AM, Firestone MK, Pett-Ridge J, Blazewicz SJ, et al. Functional succession of actively growing soil microorganisms during rewetting is shaped by precipitation history. bioRxiv. 2022; 2022.06.28.498032. More

We incorporated a macromolecular model of phytoplankton (CFM-Phyto) into the global ocean model (MITgcm). This combined model predicts cellular growth rate based on the macromolecular allocation, which in turn is used to determine the elemental stoichiometry of phytoplankton for the next model time step.The phytoplankton component of the model is executed using the following algorithm, which is illustrated graphically in Extended Data Fig. 2: (1) relate the growth rate and elemental stoichiometry of phytoplankton based on the macromolecular allocation; (2) evaluate the possible growth rates under four different limiting nutrient assumptions and select the lowest rate: Liebig’s Law of the Minimum; (3) evaluate storage of non-limiting elements; (4) evaluate excess of non-limiting elements relative to maximum quotas; (5) based on that excess, evaluate effective nutrient uptake rate; and (6) evaluate the change in the elemental stoichiometry based on the balance between the growth rate and effective nutrient uptake rate. We describe the procedural details in the following text. Parameter values are listed in Extended Data Table 1. See ref. 21 for further details and justification of each equation in CFM-Phyto; here we repeat equations essential to explain the model used in the current study.Connecting the elemental stoichiometry and the growth rateThe first step of the algorithm is to obtain the relationship between the current elemental stoichiometry and the growth rate (μ). To do that, we use CFM-Phyto21 (Extended Data Fig. 1). The model is based on the assumption of pseudo-steady state with respect to macromolecular allocation; in other words, the cellular-scale acclimation occurs rapidly relative to environmental changes. Laboratory studies show that macromolecular re-allocation occurs on the timescale of hours to days19. This is fast relative to the rates of environmental change in our coarse-resolution ocean simulations and so steady state solutions21 are used to relate growth rate, macromolecular allocation and elemental stoichiometry, as described in detail below. We first describe the case of N quota (here defined as QN; moles cellular N per mole cellular C) in detail, and then we briefly explain the case of P and C quotas as the overall procedures are similar. After that, we describe the case with Fe quota, which extends the previously published model21 for this study.Relating N quota and growth rateCFM-Phyto describes the allocation of N quota as follows, focusing on the quantitatively major molecules:$$Q_{mathrm{N}} = Q_{mathrm{N}}^{{mathrm{Pro}}} + Q_{mathrm{N}}^{{mathrm{RNA}}} + Q_{mathrm{N}}^{{mathrm{DNA}}} + Q_{mathrm{N}}^{{mathrm{Chl}}} + Q_{mathrm{N}}^{{mathrm{Sto}}}$$
(2)
where QN is total N quota (per cellular C: mol N (mol C)−1), the terms on the right-hand side are the contributions from protein, RNA, DNA, chlorophyll and N storage. We use empirically determined fixed elemental stoichiometry of macromolecules21 (Extended Data Table 1) to connect the macromolecular contributions of different elements (here C and P):$$Q_{mathrm{N}} = Q_{mathrm{C}}^{{mathrm{Pro}}}Y_{{mathrm{Pro}}}^{{mathrm{N:C}}} + Q_{mathrm{P}}^{{mathrm{RNA}}}Y_{{mathrm{RNA}}}^{{mathrm{N:P}}} + Q_{mathrm{C}}^{{mathrm{DNA}}}Y_{{mathrm{DNA}}}^{{mathrm{N:C}}} + Q_{mathrm{C}}^{{mathrm{Chl}}}Y_{{mathrm{Chl}}}^{{mathrm{N:C}}} + Q_{mathrm{N}}^{{mathrm{Nsto}}}$$
(3)
Here (Y_l^{j:k}) represents the imposed elemental ratio (elements j and k) for each macromolecular pool (l). (Q_{mathrm{C}}^x) and (Q_{mathrm{P}}^x) describe the contributions of macromolecule x to the total C quota (mol C (mol C)−1) and P quota (mol P (mol C)−1), respectively.CFM-Phyto uses the following empirically supported relationship to describe (Q_{mathrm{P}}^{{mathrm{RNA}}}) (ref. 21):$$Q_{mathrm{P}}^{{mathrm{RNA}}} = A_{{mathrm{RNA}}}^{mathrm{P}}mu Q_{mathrm{C}}^{{mathrm{Pro}}} + Q_{{mathrm{P,min}}}^{{mathrm{RNA}}}$$
(4)
where (A_{{mathrm{RNA}}}^{mathrm{P}}) is constant (below, A values represent constant except (A_{{mathrm{Chl}}}); see below), μ is growth rate (d−1) and (Q_{{mathrm{P,min}}}^{{mathrm{RNA}}}) represents the minimum amount of RNA in phosphorus per cellular C (mol P (mol C)−1). Substituting this equation into equation (3) gives:$$begin{array}{l}Q_{mathrm{N}} = Q_{mathrm{C}}^{{mathrm{Pro}}}Y_{{mathrm{Pro}}}^{{mathrm{N:C}}} + left( {A_{{mathrm{RNA}}}^{mathrm{P}}mu Q_{mathrm{C}}^{{mathrm{Pro}}} + Q_{{mathrm{P,min}}}^{{mathrm{RNA}}}} right)\Y_{{mathrm{RNA}}}^{{mathrm{N:P}}} + Q_{mathrm{C}}^{{mathrm{DNA}}}Y_{{mathrm{DNA}}}^{{mathrm{N:C}}} + Q_{mathrm{C}}^{{mathrm{Chl}}}Y_{{mathrm{Chl}}}^{{mathrm{N:C}}} + Q_{mathrm{N}}^{{mathrm{Nsto}}}end{array}$$
(5)
In CFM-Phyto, we resolve three types of protein, photosynthetic, biosynthetic and other:$$Q_{mathrm{C}}^{{mathrm{Pro}}} = Q_{mathrm{C}}^{{mathrm{Pro}}_{mathrm{Pho}}} + Q_{mathrm{C}}^{{mathrm{Pro}}_{mathrm{Bio}}} + Q_{mathrm{C}}^{{mathrm{Pro}}_{mathrm{Other}}}$$
(6)
Photosynthetic proteins represent those in chloroplasts largely responsible for light harvesting and electron transport. The model assumes a constant composition of chloroplasts; thus, the amount of photosynthetic protein is proportional to the amount of chlorophyll21:$$Q_{mathrm{C}}^{{mathrm{Pro}}_{mathrm{Pho}}} = A_{{mathrm{Pho}}}Q_{mathrm{C}}^{{mathrm{Chl}}}$$
(7)
Biosynthetic proteins represent proteins related to producing new material such as proteins, carbohydrates, lipids, RNAs, DNAs and other molecules. The models use the following empirically derived relationship21:$$Q_{mathrm{C}}^{{mathrm{Pro}}_{mathrm{Bio}}} = A_{{mathrm{Bio}}}mu$$
(8)
Substituting equations (6)–(8) (in this order) into equation (5) leads to the following equation:$$begin{array}{l}Q_{mathrm{N}} = left( {A_{{mathrm{Pho}}}Q_{mathrm{C}}^{{mathrm{Chl}}} + A_{{mathrm{Bio}}}mu + Q_{mathrm{C}}^{{mathrm{Pro}}_{mathrm{Other}}}} right)Y_{{mathrm{Pro}}}^{{mathrm{N:C}}}\ + left( {A_{{mathrm{RNA}}}^{mathrm{P}}mu left( {A_{{mathrm{Pho}}}Q_{mathrm{C}}^{{mathrm{Chl}}} + A_{{mathrm{Bio}}}mu + Q_{mathrm{C}}^{{mathrm{Pro}}_{mathrm{Other}}}} right) + Q_{{mathrm{P,min}}}^{{mathrm{RNA}}}} right)Y_{{mathrm{RNA}}}^{{mathrm{N:P}}}\ + Q_{mathrm{C}}^{{mathrm{DNA}}}Y_{{mathrm{DNA}}}^{{mathrm{N:C}}} + Q_{mathrm{C}}^{{mathrm{Chl}}}Y_{{mathrm{Chl}}}^{{mathrm{N:C}}} + Q_{mathrm{N}}^{{mathrm{Sto}}}end{array}$$
(9)
Empirically, chlorophyll depends on the growth rate and equation (10) accurately describes the relationship between the growth-rate dependences of chlorophyll under different light intensities21:$$Q_{mathrm{C}}^{{mathrm{Chl}}} = A_{{mathrm{Chl}}}mu + B_{{mathrm{Chl}}}$$
(10)
with (A_{{mathrm{Chl}}} = left( {1 + E} right)/v_I) and (B_{Chl} = m/v_I) with E (dimensionless) as a constant representing growth-rate-dependent respiration, and m (d−1) describing maintenance respiration. vI (mol C (mol C in Chl)−1 d−1) represents chlorophyll-specific photosynthesis rate based on an established function of light intensity I (μmol m−2 s−1)21,57:$$v_I = v_I^{{mathrm{max}}}left( {1 – e^{A_II}} right)$$
(11)
where (v_I^{{mathrm{max}}}) is the maximum chlorophyll-specific photosynthesis rate, e is the natural base and AI is a combined coefficient for absorption cross-section and turnover time. Substitution of equation (10) into equation (9) leads to the following quadratic relationship between QN and μ:$$Q_{mathrm{N}} = a_{mathrm{N}}mu ^2 + b_{mathrm{N}}mu + c_{mathrm{N}} + Q_{mathrm{N}}^{{mathrm{Sto}}}$$
(12)
where$$begin{array}{l}a_{mathrm{N}} = A_{{mathrm{RNA}}}^{mathrm{P}}left( {A_{{mathrm{Pho}}}A_{{mathrm{Chl}}} + A_{{mathrm{Bio}}}} right)Y_{{mathrm{RNA}}}^{{mathrm{N:P}}}\ b_{mathrm{N}} = left( {A_{{mathrm{Pho}}}A_{{mathrm{Chl}}} + A_{{mathrm{Bio}}}} right)Y_{{mathrm{Pro}}}^{{mathrm{N:C}}} + A_{{mathrm{Chl}}}Y_{{mathrm{Chl}}}^{{mathrm{N:C}}} + A_{{mathrm{RNA}}}^{mathrm{P}}left( {A_{{mathrm{Pho}}}B_{{mathrm{Chl}}} + Q_{mathrm{C}}^{{mathrm{Pro}}_{mathrm{Other}}}} right)Y_{mathrm{{RNA}}}^{{mathrm{N:P}}}\ c_{mathrm{N}} = B_{{mathrm{Chl}}}Y_{{mathrm{Chl}}}^{{mathrm{N:C}}} + left( {A_{{mathrm{Pho}}}B_{{mathrm{Chl}}} + Q_{mathrm{C}}^{{mathrm{Pro}}_{mathrm{Other}}}} right)Y_{{mathrm{Pro}}}^{{mathrm{N:C}}}\ + Q_{{mathrm{P}},{mathrm{min}}}^{{mathrm{RNA}}}Y_{{mathrm{RNA}}}^{{mathrm{N:P}}} + Q_{mathrm{C}}^{{mathrm{DNA}}}Y_{{mathrm{DNA}}}^{{mathrm{N:C}}}end{array}$$Relating P quota and growth rateSimilarly, CFM-Phyto describes the relationship between the current P quota QP and μ. P is allocated to its major molecular reservoirs:$$Q_{mathrm{P}} = Q_{mathrm{P}}^{{mathrm{RNA}}} + Q_{mathrm{C}}^{{mathrm{DNA}}}Y_{{mathrm{DNA}}}^{{mathrm{P:C}}} + Q_{mathrm{P}}^{{mathrm{Thy}}} + Q_{mathrm{P}}^{{mathrm{Other}}} + Q_{mathrm{P}}^{{mathrm{Sto}}}$$
(13)
Similar to equation (7), with the assumption of the constant composition of photosynthetic apparatus, the model connects the amount of the chlorophyll to phosphorus in thylakoid membranes:$$Q_{mathrm{P}}^{{mathrm{Thy}}} = A_{{mathrm{Pho}}}^{{mathrm{P:Chl}}}Q_{mathrm{C}}^{{mathrm{Chl}}}$$
(14)
As for N allocation, substitution of equations (14), (4), (6), (7), (8) and (10) (in this order) into equation (13) leads to a quadratic relationship between QP and μ:$$Q_{mathrm{P}} = a_{mathrm{P}}mu ^2 + b_{mathrm{P}}mu + c_{mathrm{P}} + Q_{mathrm{P}}^{{mathrm{Sto}}}$$
(15)
where$$begin{array}{l}a_{mathrm{P}} = A_{{mathrm{RNA}}}^{mathrm{P}}left( {A_{{mathrm{Pho}}}A_{{mathrm{Chl}}} + A_{{mathrm{Bio}}}} right)\ b_{mathrm{P}} = A_{{mathrm{RNA}}}^{mathrm{P}}left( {A_{{mathrm{Pho}}}B_{{mathrm{Chl}}} + Q_{mathrm{C}}^{{mathrm{Pro}}_{mathrm{Other}}}} right)Y_{{mathrm{RNA}}}^{{mathrm{N:P}}} + A_{{mathrm{Pho}}}^{{mathrm{P:Chl}}}A_{{mathrm{Chl}}}\ c_{mathrm{P}} = Q_{{mathrm{P,min}}}^{{mathrm{RNA}}} + Q_{mathrm{C}}^{{mathrm{DNA}}}Y_{{mathrm{DNA}}}^{{mathrm{P:C}}} + A_{{mathrm{Pho}}}^{{mathrm{P:Chl}}}B_{{mathrm{Chl}}} + Q_{mathrm{P}}^{{mathrm{Other}}}end{array}$$Relating C quota and growth rateSimilarly, CFM-Phyto describes C allocation as follows:$$begin{array}{l}Q_{mathrm{C}} = 1 = Q_{mathrm{C}}^{{mathrm{Pro}}} + Q_{mathrm{C}}^{{mathrm{RNA}}} + Q_{mathrm{C}}^{{mathrm{DNA}}} + Q_{mathrm{C}}^{{mathrm{Other}}} + Q_{mathrm{C}}^{{mathrm{Plip}} – {mathrm{Thy}}}\qquad + Q_{mathrm{C}}^{{mathrm{Csto}}} + Q_{mathrm{C}}^{{mathrm{Nsto}}}end{array}$$
(16)
where Plip−Thy indicates P lipid in thylakoid membranes. The equation represents the allocation per total cellular C in mol C (mol C)−1, so the sum of the macromolecules in C (QC) becomes 1. Using the imposed elemental ratios of macromolecular pools ((Y_l^{j:k})) we relate the elemental contributions:$$Q_{mathrm{C}} = Q_{mathrm{C}}^{{mathrm{Pro}}} + Q_{mathrm{P}}^{{mathrm{RNA}}}Y_{{mathrm{RNA}}}^{{mathrm{C:P}}} + Q_{mathrm{C}}^{{mathrm{DNA}}} + Q_{mathrm{C}}^{{mathrm{Other}}} + Q_{mathrm{P}}^{{mathrm{Thy}}}Y_{{mathrm{Plip}}}^{{mathrm{C:P}}} + Q_{mathrm{C}}^{{mathrm{Sto}}} + Q_{mathrm{N}}^{{mathrm{Sto}}}Y_{{mathrm{Nsto}}}^{{mathrm{C:N}}}$$
(17)
Following the steps similar to those for the N and P allocations, substituting equations (14), (4), (6), (7), (8) and (10) (in this order) into equation (17) leads to the following quadratic relationship between cellular C quota QC (=1 mol C (mol C)−1) and μ:$$Q_{mathrm{C}} = a_{mathrm{C}}mu ^2 + b_{mathrm{C}}mu + c_{mathrm{C}} + Q_{mathrm{C}}^{{mathrm{Sto}}} + Q_{mathrm{N}}^{{mathrm{Sto}}}Y_{{mathrm{Nsto}}}^{{mathrm{C:N}}}$$
(18)
where$$begin{array}{l}a_{mathrm{C}} = A_{{mathrm{RNA}}}^{mathrm{P}}left( {A_{{mathrm{Pho}}}A_{{mathrm{Chl}}} + A_{{mathrm{Bio}}}} right)Y_{{mathrm{RNA}}}^{{mathrm{C:P}}}\ b_{mathrm{C}} = A_{{mathrm{Chl}}}left( {1 + A_{{mathrm{Pho}}} + A_{{mathrm{Pho}}}^{{mathrm{P:Chl}}}Y_{{mathrm{Plip}}}^{{mathrm{C:P}}}} right) + A_{{mathrm{Bio}}} + A_{{mathrm{RNA}}}^{mathrm{P}}left( {A_{{mathrm{Pho}}}B_{{mathrm{Chl}}} + Q_{mathrm{C}}^{{mathrm{Pro}}_{mathrm{Other}}}} right)Y_{{mathrm{RNA}}}^{{mathrm{C:P}}}\ c_{mathrm{C}} = left( {1 + A_{{mathrm{Pho}}} + A_{{mathrm{Pho}}}^{{mathrm{P:Chl}}}Y_{{mathrm{Plip}}}^{{mathrm{C:P}}}} right)B_{{mathrm{Chl}}} + Q_{mathrm{C}}^{{mathrm{Pro}}_{rm{Other}}}\ + Q_{{mathrm{P}},{mathrm{min}}}^{{mathrm{RNA}}}Y_{{mathrm{RNA}}}^{{mathrm{C:P}}} + Q_{mathrm{C}}^{{mathrm{DNA}}} + Q_{mathrm{C}}^{{mathrm{Other}}}end{array}$$Relating Fe quota and growth rateIn order to capture global scale biogeochemical dynamics including the iron-limited high-nitrogen, low chlorophyll regimes, CFM-Phyto21 is extended to resolve Fe quota and allocation. The model is guided by a laboratory proteomic study58 in which the major Fe allocations are to photosystems, storage and nitrogen-fixing enzymes (nitrogenase). As we do not resolve nitrogen-fixing organisms here, Fe allocation (mol Fe (mol C)−1) represents only the first two:$$Q_{{mathrm{Fe}}} = Q_{{mathrm{Fe}}}^{{mathrm{Pho}}} + Q_{{mathrm{Fe}}}^{{mathrm{Sto}}}$$
(19)
As for equation (7) and equation (14), we relate the allocation of Fe to photosystems to the investment in chlorophyll, (Q_{mathrm{C}}^{{mathrm{Chl}}}):$$Q_{{mathrm{Fe}}}^{{mathrm{Pho}}} = A_{{mathrm{Pho}}}^{{mathrm{Fe}}}Q_{mathrm{C}}^{{mathrm{Chl}}}$$
(20)
This is a strong simplification because the pigment to photosystem ratio is observed to vary59, but enables an explicit, mechanistically motivated representation of Fe limitation, which, a posteriori, results in global scale regimes of iron limitation that resemble those observed43 (Extended Data Fig. 4). With equations (10), (19) and (20), we obtain the following relationship between QFe and μ:$$Q_{{mathrm{Fe}}} = A_{{mathrm{Pho}}}^{{mathrm{Fe}}}A_{{mathrm{Chl}}}mu + A_{{mathrm{Pho}}}^{{mathrm{Fe}}}B_{{mathrm{Chl}}} + Q_{{mathrm{Fe}}}^{{mathrm{Sto}}}$$
(21)
Evaluating the growth rateWe assume that the cellular growth rate is constrained by the most limiting element within the cell (and its associated functional macromolecules). Thus, at each time step and location, and for each cell type, the evaluation of growth rate is based on the following two steps: (1) computation of the growth rate for each element without storage; that is, the case when all of the elemental quotas are allocated to functional macromolecules; and (2) selection of the lowest growth rate among these; Liebig’s Law of the Minimum. For the first step, we define (mu _i) (i = C, N, P, Fe) as the growth rate, assuming that nutrient i is limiting. Under this condition, (Q_i^{{mathrm{Sto}}}) should be small as element i is allocated to other essential molecules. We assume that (Q_{mathrm{N}}^{{mathrm{Sto}}}) is also small under C limitation because N storage molecules are rich in carbon. With these assumptions, the solution for (mu _i) is obtained by solving the standard quadratic relationships of equations (12), (15) and (18) for N, P and C, respectively, neglecting any (Q_i^{{mathrm{Sto}}}) terms:$$mu _i = frac{{ – b_i + sqrt {b_i^2 – 4a_ileft( {c_i – Q_i} right)} }}{{2a_i}}$$
(22)
where QC = 1. For μFe, equation (21) without (Q_{{mathrm{Fe}}}^{{mathrm{Sto}}}) leads to$$mu _{{mathrm{Fe}}} = frac{{Q_{{mathrm{Fe}}} – A_{{mathrm{Pho}}}^{{mathrm{Fe}}}B_{{mathrm{Chl}}}}}{{A_{{mathrm{Pho}}}^{{mathrm{Fe}}}A_{{mathrm{Chl}}}}}$$
(23)
Once the μi values are obtained, we determine the effective growth rate, μ, based on the lowest value, which identifies the limiting element based on current intracellular quotas:$$mu = {mathrm{min}}left( {mu _{mathrm{N}},mu _{mathrm{P}},mu _{mathrm{C}},mu _{{mathrm{Fe}}}} right)$$
(24)
Evaluating nutrient storageIn CFM-Phyto, non-limiting nutrients can be stored in an intracellular reserve21, reflecting commonly observed luxury uptake. Storage is evaluated as the difference between the total elemental quota (updated later) and the functionally allocated portion of that element:$$Q_i^{{mathrm{Sto}}} = Q_i – Q_i^{{mathrm{Non}}_{mathrm{Sto}}}$$
(25)
Here (Q_i^{{mathrm{Non}}_{mathrm{Sto}}}) represents the contribution to element i by functional, non-storage molecules. For N, P and C, (Q_i^{{mathrm{Non}}_{mathrm{Sto}}}) is represented by the non-(Q_i^{{mathrm{Sto}}}) terms on the right-hand side in equations (12), (15) and (18), respectively:$$Q_i^{{mathrm{Non}}_{mathrm{Sto}}} = a_imu ^2 + b_imu + c_i$$
(26)
Similarly, for Fe, from equation (21):$$Q_{{mathrm{Fe}}}^{{mathrm{Non}}_{mathrm{Sto}}} = A_{{mathrm{Pho}}}^{{mathrm{Fe}}}A_{{mathrm{Chl}}}mu + A_{{mathrm{Pho}}}^{{mathrm{Fe}}}B_{{mathrm{Chl}}}$$
(27)
When there is N storage, (Q_{mathrm{C}}^{{mathrm{Sto}}}) must be recomputed to consider the allocation of C to it:$$Q_{mathrm{C}}^{{mathrm{Sto}}} = Q_{mathrm{C}} – Q_{mathrm{C}}^{{mathrm{Non}}_{mathrm{Sto}}} – Q_{mathrm{N}}^{{mathrm{Sto}}}Y_{{mathrm{Nsto}}}^{{mathrm{C:N}}}$$
(28)
Evaluating the excess nutrientStorage capacity for any element is finite and we define excess nutrient as a nutrient (N, P, Fe) that is in beyond an empirically informed, imposed maximum phytoplankton storage capacity. Excess nutrient is assumed to be excreted (see below). Excess of element i ((Q_i^{{mathrm{Exc}}})) is computed:$$Q_i^{{mathrm{Exc}}} = {mathrm{max}}left( {Q_i – Q_i^{{mathrm{max}}},0} right)$$
(29)
where (Q_i^{{mathrm{max}}}) is maximum capacity for nutrient i. For N, CFM-Phyto computes (Q_i^{{mathrm{max}}}) as a sum of non-storage molecules and prescribed maximum nutrient storing capacity according to model–data comparison21:$$Q_i^{{mathrm{max}}} = Q_i^{{mathrm{Non}}_{mathrm{Sto}}} + Q_i^{{mathrm{Sto}}_{mathrm{max}}}$$
(30)
Laboratory studies suggest that when P is not limiting, the phosphorus quota maximizes to a value that is almost independent of growth rate21,39,44. Storage of each element is finite and the upper limit to storage is imposed by specifying the maximum cellular quotas ((Q_{mathrm{P}}^{{mathrm{max}}}) (ref. 21) and also (Q_{{mathrm{Fe}}}^{{mathrm{max}}})) with size and taxonomic dependencies (for example, refs. 27,41). Thus, the maximum storage is represented by the difference between the prescribed maximum quota and the actual quota21:$$Q_i^{{mathrm{Sto}}_{mathrm{max}}} = Q_i^{{mathrm{max}}} – Q_i$$
(31)
In the case where (Q_i^{{mathrm{Sto}}}) computed in the previous section exceeds (Q_i^{{mathrm{Sto}}_{mathrm{max}}}), the value of (Q_i^{{mathrm{Sto}}}) is replaced by (Q_i^{{mathrm{Sto}}_{mathrm{max}}}) and the difference is placed in the excess pool, (Q_i^{{mathrm{Exc}}}).Computing effective nutrient uptake rateOne factor that influences the cellular elemental quota is the effective nutrient uptake rate (mol i (mol C)−1 d−1) of N, P and Fe, which we define as follows:$$V_i^{{mathrm{Eff}}} = V_i – frac{{Q_i^{{mathrm{Exc}}}}}{{tau _i^{{mathrm{Exu}}}}}$$
(32)
where Vi (mol i (mol C)−1 d−1) is nutrient uptake rate and the second term represents the exudation of the excess nutrient based on the timescale (tau _i^{{mathrm{Exu}}}) (d−1). For Vi, we use Monod kinetics60,61:$$V_i = V_i^{{mathrm{max}}}frac{{[i]}}{{left[ i right] + K_i}}$$
(33)
where (V_i^{{mathrm{max}}}) is maximum nutrient uptake, [i] (mmol m−3) is the environmental concentration of nutrient i and Ki (mmol m−3) is the half-saturation constant of i. Previous models have resolved the relationship between nutrient uptake and allocation to transporters31,62. Here we do not explicitly resolve allocation to transporters, as proteomic studies indicate that it is a relatively minor component of the proteome compared with photosystems and biosynthesis in phytoplankton63. Transporter proteins could be represented in a model with a finer-scale resolution of the proteome64.Differentiating small and large phytoplanktonIn this model, ‘small’ phytoplankton broadly represent picocyanobacteria, which have high nutrient affinities and low maximum growth rates (for example, Prochlorococcus), whereas ‘large’ phytoplankton represent eukaryotes with higher maximum growth rates (for example, diatoms). The former are associated with a gleaner strategy adapted to oligotrophic regimes, while the latter are opportunistic, adapted to variable and nutrient-enriched regimes. To encapsulate this, the large phytoplankton have overall higher imposed (V_i^{{mathrm{max}}}) (~µmaxQi), Ki and (v_I^{mathrm{max}}) than for the small phytoplankton (Extended Data Table 1), consistent with the previous models (for example, ref. 10). In addition, the larger cells are assigned a higher (Q_{mathrm{P}}^{{mathrm{max}}}) following the observed trends (Fig. 1 and Extended Data Table 1).Computing the change in the elemental stoichiometryThe computation of the change in the elemental quotas is done based on the balance between the effective nutrient uptake rate and the loss of nutrient to the new cells:$$frac{{{mathrm{d}}Q_i}}{{{mathrm{d}}t}} = V_i^{{mathrm{Eff}}} – mu Q_i$$
(34)
This change in the elemental quotas based on the cellular processes and the passive transport of elements in phytoplankton by the flow field created by MITgcm governs the elemental stoichiometry of the next time step at a specific grid box, as in other versions of ecological models with MITgcm10.Calculation of CV valuesWe computed the CV values based on the following equation:$${mathrm{CV}} = frac{sigma }{{bar x}}$$
(35)
where σ is the standard deviation and (bar x) is the mean. The purpose of this computation is to quantify the latitudinal variation of the averaged elemental stoichiometry. Thus, we used the averaged values for each latitude (as plotted in Fig. 2) for the calculation of σ and (bar x).MITgcm-CFMThe biogeochemical and ecological component of the model resolves the cycling of C, P, N and Fe through inorganic, living, dissolved and particulate organic phases. The biogeochemical and biological tracers are transported and mixed by the MIT general circulation model (MITgcm)35,36, constrained to be consistent with altimetric and hydrographic observations (the ECCO-GODAE state estimates)65. This three-dimensional configuration has a coarse resolution (1° × 1° horizontally) and 23 depth levels ranging from 5 m at the surface to 5450 m at depth. The model was run for three years, and the results of the third year were analysed, by which time the modelled plankton distribution becomes quasi-stable. Equations for the biogeochemical processes are as described by equations and parameters in previous studies10,38. Here, however, we include only nitrate for inorganic nitrogen, and do not resolve the silica cycle. We simulated eukaryotic and prokaryotic analogues of phytoplankton (as ‘large’ and ‘small’ phytoplankton). The eukaryotic analogue has a higher maximum C fixation rate for the same macromolecular composition and higher maximum nutrient uptake rates, but also has overall higher half-saturation constants for nutrient uptake. We used light absorption spectra of picoeukaryotes, which sits in-between small prokaryotes and large eukaryotes10. In MITgcm, the mortality of phytoplankton is represented by the sum of a linear term (ml), representing sinking and maintenance losses, and quadratic terms representing the action of unresolved next-trophic levels66,67, implicitly assuming that the higher-trophic-level biomass scales with that of its prey. We assumed that the latter term is small to avoid introducing additional uncertainties. Similarly, we do not resolve the stoichiometric effects of prey selection due to the nutritional status of prey, or viral partitioning of nutrients in the environment50. Atmospheric iron deposition varies by orders of magnitude around the globe and has a large margin of uncertainty, including the bio-availability of the deposited iron, which in turn depends on the source and chemical history of the deposited material68. To realize a realistic global net primary production, we doubled the atmospheric iron input from ref. 10; this factor is well within the uncertainty of the iron supply estimates. Each of the two phytoplankton groups has variable C:N:P:Fe as determined by the component macromolecules at each time step. The pools of C, N, P and Fe are tracked within the modelled three-dimensional flow fields. More

HFRS distribution in China, 2004–2018From January 1, 2004 to December 31, 2018, 190 203 cases of HFRS were reported nationwide in China, with an average annual incidence rate of 0.950 per 100,000 people, with the highest incidence in 2004 (1.926 per 100,000) and the lowest in 2018 (0.86 per 100,000) (Fig. 1A), and the cases showed obvious seasonal fluctuations (Fig. 1B). HFRS cases existed every month and showed an obvious dual-season mode every year, with a spring peak from May to June and a winter peak from November to December. The highest number of cases were in May and November, with the composition ratios accounting of 9.51% and 17.06%, respectively (Fig. 1B).Figure 1The incidence and number of HFRS cases reported in China, 2004–2018. (A) Number of cases and incidence by year. Trend of the incidence rate of HFRS between 2004 and 2018 shown by the joinpoint regression (upper right corner). The red squares represent the observed crude incidence of HFRS and the lines represent the slope of the annual percentage change (APC). (B) The pink line represents the monthly incidence of HFRS. The bar chart shows the number of cases at peak and trough.Full size imageThe incidence of HFRS in northern regions was higher than that in the south, especially in Heilongjiang, Liaoning, Jining, Shaanxi, Shandong and Hebei provinces. Relatively few cases existed in south China, which were mainly concentrated in Jiangxi, Zhejiang, Hunan and Fujian (Figs. S1 and S2). Spatial autocorrelation analysis indicated that HFRS cases were positively correlated (Moran’s I = 0.09, p  More

Tallis, H. M. et al. Front. Ecol. Environ. 16, 563–570 (2018).Article 

Google Scholar 
Tu, C. et al. Nat. Sustain. https://doi.org/10.1038/s41893-022-01008-1 (2022).Article 

Google Scholar 
Lloyd, W. F. Two lectures on the checks to population (S. Collingwood, The University of Oxford, 1833).Gordon, H. S. J. Polit. Econ. 62, 124–142 (1954).Article 

Google Scholar 
Nash, J. F. Jr Proc. Natl Acad. Sci. USA 36, 48–49 (1950).Article 
CAS 

Google Scholar 
Hardin, G. Science 162, 1243–1248 (1968).Article 
CAS 

Google Scholar 
Ostrom, E. Governing the Commons: The Evolution of Institutions for Collective Action (Cambridge University Press, 1990).Sethi, R. & Somanathan, E. Am. Econ. Rev. 86, 766–788 (1996).
Google Scholar 
Schill, C. et al. Nat. Sustain. 2, 1075–1082 (2019).Article 

Google Scholar 
Levin, S. Philos. Trans R. Soc. B Biol. Sci. 365, 13–18 (2010).Article 

Google Scholar  More

Chironomids have the ability to survive and reproduce in polluted environments, and thus they are included in many ecological studies where approaches may be taxonomic or functional16. The diversity of most macroinvertebrates is controlled by the oxygen level of water, but chironomids may survive in hypoxic conditions where the oxygen concentration may be less than 3 mg l−117. The current study demonstrates that changing seasons, as well as anthropogenic activities, have a significant impact on the levels of DO in aquatic bodies. As observed from the result, DO highly influences the tube length of the chironomid larvae. Since KWC is a wastewater canal, the average oxygen level is lower (5.24 ± 1.14 mg l−1) than KFP (6.63 ± 1.28 mg l−1) which is a normal fish culturing pond. It has also been observed that the average tube length of the chironomid larvae of KWC (8.66 ± 0.88 mm) is higher than KFP (7.68 ± 0.62 mm), which indicates that a low concentration of DO promotes the building of longer tubes in natural conditions. Similar observations were also observed in laboratory conditions. When the oxygen level (7.03 ± 0.41 mg l−1) in the experiment was kept in the normal range, there was negligible variation in tube length (7.61 ± 0.31 mm). But when the concentration of oxygen is gradually reduced by dilution, the tube length starts to increase accordingly, which is explained graphically in Fig. 4. The regression model of both the experimental conditions also supports the hypothesis that the tube length has an inverse relationship with DO. The scatter plot and simple linear regression confirmed the inverse relationship between DO and tube length (Figs. 1 and 2).Chironomid larvae are able to grow in the polluted water of a wastewater pond as dominant macroinvertebrates18. It is observed that those larvae living in the sand tubes are more susceptible to chemical pollutants than the larvae living in silt tubes7. Sand particles are bigger than silt and are not suitable for the survival of larvae19. Chironomus riparius larvae make their tubes from different external particles and their own proteins20. Midge larvae are the inhabitants of sediments, and at the same time, sediment is the depository of different inorganic, organic, and heavy metals. In such cases, the tube of chironomid larvae may act as a defensive structure, which protects them from the adverse effects of undesirable pollutants and may increase their tolerance against such chemicals21,22,23.Larvae can thrive in benthic sediments with high decaying organic content and very low DO concentrations in water bodies24. In poor DO concentration, larvae can survive due to the presence of haemoglobin in their body tissue fluid, which plays an important physiological role in increasing respiratory efficiency, as was observed in Chironomus plumosus. Longer tube length may help larvae generate better respiratory currents so that they can cope with a low DO environment.Tube length is crucial for living in water because primarily tubes protect them from outer environmental factors like predators, and pollution. It was observed during this study that when the DO of water is low, larvae make elongated tubes to reach the upper layer of water, where the DO level is comparatively high. To get their required amount of oxygen, the larvae increase the tube length towards the water surface and increase the DO in tube water by undulating the body and other structures, creating a current inside the tube25,26. On contrary, when the DO level of the surrounding water of chironomid is sufficient, they can manage their normal physiological activities with the available oxygen. They need not to elongate their tube length. That’s why their tube length is inversely related to the DO of their surrounding medium.If tube length does not increase in size in hypoxic water, larvae will not be able to meet their oxygen demand. If the DO of water decreases, tube length will increase and vice versa. Behavioural and physiological adaptations of chironomids larvae make them successful to live in a hypoxic environment. Thus, in hypoxic conditions, larvae with longer tubes are able to gather more oxygen from the upper layer of water and get more space to create a current of water to increase the amount of O2 inside the tube by undulating the preanal papillae, anal gill, ventral gills. This would explain why the tube length of chironomids depends on the DO of water. Hence by measuring the tube length with a standard measuring scale, one may get an idea about the quality of water, especially DO, before doing any chemical analysis. The work seems to be unique and novel for its own kind. More