Philippa Kaur

Site description, chronology, and sampling
A detailed description of the site, coring methods, age-depth model reconstruction, and sampling strategy can be found in Alsos et al.41. Briefly, Lake Øvre Æråsvatnet is located on Andøya, Northern Norway (69.25579°N, 16.03517°E) (Fig. 1a, b). In 2013, two cores were collected from the deepest sediments, AND10 and AND11, which were stored at 4 °C prior to sampling. Macrofossil remains were dated, with those from AND10 all dating to within the LGM. For the longer core AND11, a Bayesian age-depth model was required to estimate the age of each layer41. In this study, we selected one sample of LGM sediments from each of the two cores. According to the Bayesian age-depth model, sample Andøya_LGM_B, from 1102 cm depth in AND11, was dated to a median age of 17,700 (range: 20,200–16,500) cal yr BP. The age of Andøya_LGM_A, from 938 cm depth in AND10, was estimated at 19,500 cal yr BP, based on the interpolated median date between two adjacent macrofossils (20 cm above: 19,940–18,980 cal yr BP, 30 cm below: 20,040–19,000 cal yr BP). As Andøya_LGM_A falls within the age range of Andøya_LGM_B, we consider the samples to be broadly contemporaneous.
Sampling, DNA extraction, library preparation, and sequencing
The two cores were subsampled at the selected layers under clean conditions, in a dedicated ancient DNA laboratory at The Arctic University Museum of Norway in Tromsø. We extracted DNA from 15 g of sediment following the Taberlet phosphate extraction protocol29 in the same laboratory. We shipped a 210 µL aliquot of each DNA extract to the ancient DNA dedicated laboratories at the Centre for GeoGenetics (University of Copenhagen, Denmark) for double-stranded DNA library construction. The extracts were concentrated to 80 µL using Amicon Ultra-15 30 kDa centrifugal filters (Merck Millipore, Darmstadt, Germany) and half of each extract (40 µL, totalling between 31.7 and 36.0 ng of DNA) was converted into Illumina-compatible libraries using established protocols10. Each library was dual-indexed via 12 cycles of PCR. The libraries were then purified using the AmpureBead protocol (Beckman Coulter, Indianapolis, IN, USA), adjusting the volume ratio to 1:1.8 library:AmpureBeads, and quantified using a BioAnalyzer (Agilent, Santa Clara, CA, USA). The indexed libraries were pooled equimolarly and sequenced on a lane of the Illumina HiSeq 2500 platform using 2 × 80 cycle paired-end chemistry.
Raw read filtering
For each sample, we merged and adapter-trimmed the paired-end reads with SeqPrep (https://github.com/jstjohn/SeqPrep/releases, v1.2) using default parameters. We only retained the resulting merged sequences, which were then filtered with the preprocess function of the SGA toolkit v0.10.15 (ref. 57) by the removal of those shorter than 35 bp or with a DUST complexity score > 1.
Metagenomic analysis of the sequence data
We first sought to obtain an overview of the taxonomic composition of the samples and therefore carried out a BLAST-based metagenomic analysis on the two filtered sequence datasets. To make the datasets more computationally manageable, we subsampled the first and last one-million sequences from the filtered dataset of each sample and analysed each separately. The data subsets were each identified against the NCBI nucleotide database (release 223) using the blastn function from the NCBI-BLAST+ suite v2.2.18+58 under default settings. For each sample, the results from the two subsets were checked for internal consistency, merged into one dataset, and loaded into MEGAN v6.12.3 (ref. 59). Analysis and visualization of the Last Common Ancestor (LCA) was carried out for the taxonomic profile using the following settings: min score = 35, max expected = 1.0E−5, min percent identity = 95%, top percent = 10%, min support percentage = 0.01, LCA = naive, min percent sequence to cover = 95%. We define sequences as the reads with BLAST hits assigned to taxa post-filtering, thus ignoring “unassigned” and “no hit” categories.
Alignment to reference genome panels
We mapped our filtered data against three different reference panels to help improve taxonomic identifications and provide insight into the sequence abundance of the identified taxa (Supplementary Data 2 and 3). The first reference panel consisted of 42 nuclear genomes that included genera expected from Northern Norway, exotic/implausible taxa for LGM Andøya, human, six Nannochloropsis species, and four strains of Mycobacterium. The inclusion of exotic taxa was to give an indication of the background spurious mapping rate, which can result from mappings to conserved parts of the genome and/or short and damaged ancient DNA molecules22,23. We included Nannochloropsis, Mycobacterium, and human genomes, due to their overrepresentation in the BLAST-based metagenomic analysis. The other two reference panels were based on either all 8486 mitochondrial or 2495 chloroplast genomes on NCBI GenBank (as of January 2018). The chloroplast dataset was augmented with 247 partial or complete chloroplast genomes generated by the PhyloNorway project60 for 2742 chloroplast genomes in total. The filtered data were mapped against each reference genome or organellar genome set individually using bowtie2 v2.3.4.1 (ref. 61) under default settings. The resulting bam files were processed with SAMtools v0.1.19 (ref. 62). We removed unmapped sequences with SAMtools view and collapsed PCR duplicate sequences with SAMtools rmdup.
For the nuclear reference panel, we reduced potential spurious or nonspecific sequence mappings by comparing the mapped sequences to both the aligned reference genome and the NCBI nucleotide database using NCBI-BLAST+, following the method used by Graham et al.9, as modified by Wang et al.11. The sequences were aligned using the following NCBI-BLAST+ settings: num_alignments = 100 and perc_identity = 90. Sequences were retained if they had better alignments, based on bit score, to reference genomes as compared to the NCBI nucleotide database. If a sequence had a better or equal match against the NCBI nucleotide database, it was removed, unless the LCA of the highest NCBI nucleotide bit score was from the same genus as the reference genome (based on the NCBI taxonID). To standardize the relative mapping frequencies to genomes of different size, we calculated the number of retained mapped sequences per Mb of genome sequence.
The sequences mapped against the chloroplast and mitochondrial reference panels were filtered and reported in a different manner than the nuclear genomes. First, to exclude any non-eukaryotic sequences, we used NCBI-BLAST+ to search sequence taxonomies and retained sequences if the LCA was, or was within, Eukaryota. Second, for the sequences that were retained, the LCA was calculated and reported in order to summarize the mapping results across the organelle datasets. LCAs were chosen as the reference sets are composed of multiple genera.
Within the Nannochloropsis nuclear reference alignments, the relative mapping frequency was highest for N. limnetica. In addition, the relative mapping frequency for other Nannochloropsis taxa was higher than those observed for the exotic taxa. This could represent the mapping of sequences that are conserved between Nannochloropsis genomes or suggest the presence of multiple Nannochloropsis taxa in a community sample. We therefore cross-compared mapped sequences to determine the number of uniquely mapped sequences per reference genome. First, we individually remapped the filtered data to six available Nannochloropsis nuclear genomes, the accession codes of which are provided in Supplementary Data 2. For each sample, we then calculated the number of sequences that uniquely mapped to, or overlapped, between each Nannochloropsis genome. We repeated the above analysis with six available chloroplast sequences (Supplementary Data 2) to get a comparable overlap estimation for the chloroplast genome.
Reconstruction of the Andøya Nannochloropsis community organellar palaeogenomes
To place the Andøya Nannochloropsis community taxon into a phylogenetic context, and provide suitable reference sequences for variant calling, we reconstructed environmental palaeogenomes for the Nannochloropsis mitochondria and chloroplast. First, the raw read data from both samples were combined into a single dataset and re-filtered with the SGA toolkit to remove sequences shorter than 35 bp, but retain low complexity sequences to assist in the reconstruction of low complexity regions in the organellar genomes. This re-filtered sequence dataset was used throughout the various steps for environmental palaeogenome reconstruction.
The re-filtered sequence data were mapped onto the N. limnetica reference chloroplast genome (NCBI GenBank accession: NC_022262.1) with bowtie2 using default settings. SAMtools was used to remove unmapped sequences and PCR duplicates, as above. We generated an initial consensus genome from the resulting bam file with BCFtools v1.9 (ref. 62), using the mpileup, call, filter, and consensus functions. For variable sites, we produced a majority-rule consensus using the –variants-only and –multiallelic-caller options, and for uncovered sites the reference genome base was called. The above steps were repeated until the consensus could no longer be improved. The re-filtered sequence data was then remapped onto the initial consensus genome sequence with bowtie2, using the above settings. The genomecov function from BEDtools v2.17.0 (ref. 63) was used to identify gaps and low coverage regions in the resulting alignment.
We attempted to fill the identified gaps, which likely consisted of diverged or difficult-to-assemble regions. For this, we assembled the re-filtered sequence dataset into de novo contigs with the MEGAHIT pipeline v1.1.4 (ref. 64), using a minimum k-mer length of 21, a maximum k-mer length of 63, and k-mer length increments of six. The MEGAHIT contigs were then mapped onto the initial consensus genome sequence with the blastn tool from the NCBI-BLAST+ toolkit. Contigs that covered the gaps identified by BEDtools were incorporated into the initial consensus genome sequence, unless a blast comparison against the NCBI nucleotide database suggested a closer match to non-Nannochloropsis taxa. We repeated the bowtie2 gap-filling steps iteratively, using the previous consensus sequence as reference, until a gap-free consensus was obtained. The re-filtered sequence data were again mapped, the resulting final assembly was visually inspected, and the consensus was corrected where necessary. This was to ensure the fidelity of the consensus sequence, which incorporated de novo-assembled contigs that could potentially be problematic, due to the fragmented nature and deaminated sites of ancient DNA impeding accurate assembly65.
Annotation of the chloroplast genome was carried out with GeSeq v1.77 (ref. 66), using the available annotated Nannochloropsis chloroplast genomes (accession codes provided in Supplementary Table 7). The resulting annotated chloroplast was visualized with OGDRAW v1.3.1 (ref. 67).
The same assembly and annotation methods outlined above were used to reconstruct the mitochondrial palaeogenome sequence, where the initial mapping assembly was based on the N. limnetica mitochondrial sequence (NCBI GenBank accession: NC_022256.1). The final annotation was carried out by comparison against all available annotated Nannochloropsis mitochondrial genomes (accession codes provided in Supplementary Table 7).
If the Nannochloropsis sequences derived from more than one taxon, then alignment to the N. limnetica chloroplast genome could introduce reference bias, which would underestimate the diversity of the Nannochloropsis sequences present. We therefore reconstructed Nannochloropsis chloroplast genomes, but using the six available Nannochloropsis chloroplast genome sequences, including N. limnetica, as seed genomes (accession codes for the reference genomes are provided in Supplementary Table 3). The assembly of the consensus sequences followed the same method outlined above, but with two modifications to account for the mapping rate being too low for complete genome reconstruction based on alignment to the non-N. limnetica reference sequences. First, consensus sequences were called with SAMtools, which does not incorporate reference bases into the consensus at uncovered sites. Second, neither additional gap filling nor manual curation was implemented.
Analysis of ancient DNA damage patterns
We checked for the presence of characteristic ancient DNA damage patterns for sequences aligned to three nuclear genomes: human, Nannochloropsis limnetica and Mycobacterium avium. We further analysed damage patterns for sequences aligned to both the reconstructed N. limnetica composite organellar genomes. Damage analysis was conducted with mapDamage v2.0.8 (ref. 68) using the following settings: –merge-reference-sequences and –length = 160.
Assembly of high- and low-frequency variant consensus sequences
The within-sample variants in each reconstructed organellar palaeogenome was explored by creating two consensus sequences, which included either high- or low-frequency variants at multiallelic sites. For each sample, the initial filtered sequence data were mapped onto the reconstructed Nannochloropsis chloroplast palaeogenome sequence with bowtie2 using default settings. Unmapped and duplicate sequences were removed with SAMtools, as above. We used the BCFtools mpileup, call, and normalize functions to identify the variant sites in the mapped dataset, using the –skip-indels, –variants-only, and –multiallelic-caller options. The resulting alleles were divided into two sets, based on either high- or low-frequency variants. High-frequency variants were defined as those present in the reconstructed reference genome sequence. Both sets were further filtered to only include sites with a quality score of 30 or higher and a coverage of at least half the average coverage of the mapping assembly (minimum coverage: Andøya_LGM_A = 22×, Andøya_LGM_B = 14×). We then generated the high- and low-frequency variant consensus sequences using the consensus function in BCFTools. The above method was repeated for the reconstructed Nannochloropsis mitochondrial genome sequence in order to generate comparable consensus sequences of high- and low-frequency variants (minimum coverage: Andøya_LGM_A = 16×, Andøya_LGM_B = 10×).
Phylogenetic analysis of the reconstructed organellar palaeogenomes
We determined the phylogenetic placement of our high- and low-frequency variant organellar palaeogenomes within Nannochloropsis, using either full mitochondrial and chloroplast genome sequences or three short loci (18S, ITS, rbcL). We reconstructed the 18S and ITS1-5.8S-ITS2 complex using DQ977726.1 (full length) and EU165325.1 (positions 147:1006, corresponding to the ITS complex) as seed sequences following the same approach that was used for the organellar palaeogenome reconstructions, except that the first and last 10 bp were trimmed to account for the lower coverage due to sequence tiling. We then called high and low variant consensus sequences as described above.
We created six alignments using available sequence data from NCBI GenBank (Supplementary Data 4) with the addition of: (1 + 2) the high- and low-frequency variant chloroplast or mitochondrial genome consensus sequences, (3) an ~1100 bp subset of the chloroplast genome for the rbcL alignment, (4 + 5) ~1800 and ~860 bp subsets of the nuclear multicopy complex for the 18S and ITS alignments, respectively, and (6) the reconstructed chloroplast genome consensus sequences derived from the alternative Nannochloropsis genome starting points. Full details on the coordinates of the subsets are provided in Supplementary Data 4. We generated alignments using MAFFT v7.427 (ref. 69) with the maxiterate = 1000 setting, which was used for the construction of a maximum likelihood tree in RAxML v8.1.12 (ref. 70) using the GTRGAMMA model and without outgroup specified. We assessed branch support using 1000 replicates of rapid bootstrapping.
Nannochloropsis variant proportions and haplogroup diversity estimation
To estimate major haplogroup diversity, we calculated the proportions of high and low variants in the sequences aligned to our reconstructed Nannochloropsis mitochondrial and chloroplast genomes. For each sample, we first mapped the initial filtered sequence data onto the high- and low-frequency variant consensus sequences with bowtie2. To avoid potential reference biases, and for each organellar genome, the sequence data were mapped separately against both frequency consensus sequences. The resulting bam files were then merged with SAMtools merge. We removed exact sequence duplicates, which may have been mapped to different coordinates, from the merged bam file by randomly retaining one copy. This step was replicated five times to examine its impact on the estimated variant proportions. After filtering, remaining duplicate sequences—those with identical mapping coordinates—were removed with SAMtools rmdup. We then called variable sites from the duplicate-removed bam files using BCFTools under the same settings as used in the assembly of the high- and low-frequency variant consensus sequences. We restricted our analyses to transversion-only variable positions to remove the impact of ancient DNA deamination artifacts. For each variable site, the proportion of reference and alternative alleles was calculated, based on comparison to the composite N. limnetica reconstructed organellar palaeogenomes. We removed rare alleles occurring at a proportion of More

Viruses and mammalian data
Viral genomic data
Complete sequences of coronaviruses were downloaded from Genbank44. Sequences labelled with the terms: ‘vaccine’, ‘construct’, ‘vector’, ‘recombinant’ were removed from the analyses. In addition, we removed those associated with experimental infections where possible. This resulted in a total of 3264 sequences for 411 coronavirus species or strains (i.e., viruses below species level on NCBI taxonomy tree). Of those, 88 were sequences of coronavirus species, and 307 sequences of strains (in 25 coronavirus species, with total number of species included = 92). Of our included species, six in total were unclassified Coronavirinae (unclassified coronaviruses).
Selection of potential mammalian hosts of coronaviruses
We processed meta-data accompanying all sequences (including partial sequences but excluding vaccination and experimental infections) of coronaviruses uploaded to GenBank to extract information on hosts (to species level) of these coronaviruses. We supplemented these data with species-level hosts of coronaviruses extracted from scientific publications via the ENHanCEd Infectious Diseases Database (EID2)45. This resulted in identification of 313 known terrestrial mammalian hosts of coronaviruses (regardless of whether a complete genome was available or not, n = 185 mammalian species for which an association with a coronavirus with complete genome was identified). We expanded this set of potential hosts by including terrestrial mammalian species in genera containing at least one known host of coronavirus, and which are known to host one or more other virus species (excluding coronaviruses, information of whether the host is associated with a virus were obtained from EID2). This results in total of 876 mammalian species which were selected.
Quantification of viral similarities
We computed three types of similarities between each two viral genomes as summarised below.
Biases and codon usage
We calculated proportion of each nucleotide of the total coding sequence length. We computed dinucleotide and codon biases46 and codon-pair bias, measured as the codon-pair score46,47 in each of the above sequences. This enabled us to produce for each genome sequence (n = 3264) the following feature vectors: nucleotide bias, dinucleotide bias, codon biased and codon-pair bias.
Secondary structure
Following alignment of sequences (using AlignSeqs function in R package Decipher48), we predicted the secondary structure for each sequence using PredictHEC function in the R package Decipher48. We obtained both states (final prediction), and probability of secondary structures for each sequence. We then computed for each 1% of the genome length both the coverage (number of times a structure was predicted) and mean probability of the structure (in the per cent of the genome considered). This enabled us to generate six vectors (length = 100) for each genome representing: mean probability and coverage for each of three possible structures—Helix (H), Beta-Sheet (E) or Coil (C).
Genome dissimilarity (distance)
We calculated pairwise dissimilarity (in effect a hamming distance) between each two sequences in our set using the function DistanceMatrix in the R package Decipher48. We set this function to penalise gap-to-gap and gap-to-letter mismatches.
Similarity quantification
We transformed the feature (traits) vectors described above into similarities matrices between coronaviruses (species or strains). This was achieved by computing cosine similarity between these vectors in each category (e.g., codon-pair usage, H coverage, E probability). Formally, for each genomic feature (n = 10) presented by vector as described above, this similarity was calculated as follows:

$${mathrm{sim}}_{{mathrm{genomic}}_{mathrm{l}}}left( {s_m,s_n} right) = {mathrm{sim}}_{{mathrm{genomic}}_{mathrm{l}}}left( {{mathbf{V}}_m^{f_l},{mathbf{V}}_n^{f_l}} right) = frac{{mathop {sum}nolimits_{i = 1}^d {left( {{mathbf{V}}_m^{f_l}[i] times {mathbf{V}}_n^{f_l}[i]} right)} }}{{sqrt {mathop {sum}nolimits_{i = 1}^d {{mathbf{V}}_m^{f_l}[i]^2} } times sqrt {mathop {sum}nolimits_{i = 1}^d {{mathbf{V}}_n^{f_l}[i]^2} } }}$$
(1)

where sm and sn are two sequences presented by two feature vectors ({mathbf{V}}_m^{f_l}) and ({mathbf{V}}_n^{f_l}) from the genomic feature space fl (e.g., codon-pair bias) of the dimension d (e.g., d = 100 for H coverage).
We then calculated similarity between each pair of virus strains or species (in each category) as the mean of similarities between genomic sequences of the two virus strains or species (e.g., mean nucleotide bias similarity between all sequences of SARS-CoV-2 and all sequences of MERS-CoV presented the final nucleotide bias similarity between SARS-CoV-2 and MERS-CoV). This enabled us to generate 11 genomic features similarity matrices (the above 10 features represented by vectors and genomic similarity matrix) between our input coronaviruses. Supplementary Fig. 1 illustrates the process.
Similarity network fusion (SNF)
We applied SNF49 to integrate the following similarities in order to reduce our viral genomic feature space: (1) nucleotide, dinucleotide, codon and codon-pair usage biases were combined into one similarity matrix—genome bias similarity. And (2) Helix (H), Beta-Sheet (E) or Coil (C) mean probability and coverage similarities (six in total) were combined into one similarity matrix—secondary structure similarity.
SNF applies an iterative nonlinear method that updates every similarity matrix according to the other matrices via nearest neighbour approach and is scalable and is robust to noise and data heterogeneity. The integrated matrix captures both shared and complementary information from multiple similarities.
Quantification of mammalian similarities
We calculated a comprehensive set of mammalian similarities. Table 3 summarises these similarities and provides justification for inclusion. Supplementary Note 1 provides full details.
Table 3 Mammalian phylogenetic, ecological and geospatial similarities.
Full size table

Quantification of network similarities
Network construction
We processed meta-data accompanying all sequences (including partial genome but excluding vaccination and experimental infections) of coronaviruses uploaded to Genbank44 (accessed 4 May 2020) to extract information on hosts (to species level) of these coronaviruses. We supplemented these data with virus–host associations extracted from publications via the EID2 Database45. This resulted in 1669 associations between 1108 coronaviruses and 545 hosts (including non-mammalian hosts). We transformed these associations into a bipartite network linking species and strains of coronaviruses with their hosts.
Quantification of topological features
The above constructed network summarises our knowledge to date of associations between coronaviruses and their hosts, and its topology expresses patterns of sharing these viruses between various hosts and host groups. Our analytical pipeline captures this topology, and relations between nodes in our network, by means of node embeddings. This approach encodes each node (here either a coronavirus or a host) with its own vector representation in a continuous vector space, which, in turn, enables us to calculate similarities between two nodes based on this representation.
We adopted DeepWalk23 to compute vectorised representations for our coronaviruses and hosts from the network connecting them. DeepWalk23 uses truncated random walks to get latent topological information of the network and obtains the vector representation of its nodes (in our case coronaviruses and their hosts) by maximising the probability of reaching a next node (i.e., probability of a virus–host association) given the previous nodes in these walks (Supplementary Note 2 lists further details).
Similarity calculations
Following the application of DeepWalk to compute the latent topological representation of our nodes, we calculated the similarity between two nodes in our network—n (vectorised as N) and m (vectorised as M), by using cosine similarity as follows24,25:

$${mathrm{sim}}_{{mathrm{network}}}left( {n,m} right) = {mathrm{sim}}_{{mathrm{network}}}left( {{mathbf{M}},{mathbf{N}}} right) = frac{{mathop {sum}nolimits_{i = 1}^d {left( {m_i times n_i} right)} }}{{sqrt {mathop {sum}nolimits_{i = 1}^d {m_i^2} } times sqrt {mathop {sum}nolimits_{i = 1}^d {n_i^2} } }}$$
(2)

where d is the dimension of the vectorised representation of our nodes: M, N; and mi and ni are the components of vectors M and N, respectively.
Similarity learning meta-ensemble—a multi-perspective approach
Our analytical pipeline stacks 12 similarity learners into testable meta-ensembles. The constituent learners can be categorised by the following three ‘points of view’ (see also Supplementary Fig. 4 for a visual description):
Coronaviruses—the virus point of view
We assembled three models derived from (a) genome similarity, (b) genome biases and (c) genome secondary structure. Each of these learners gave each coronavirus–mammalian association (( {v_i to m_j} )) a score, termed confidence, based on how similar the coronavirus vi is to known coronaviruses of mammalian species mj, compared to how similar vi is to all included coronaviruses. In other words, if vi is more similar (e.g., based on genome secondary structure) to coronaviruses observed in host mj than it is similar to all coronaviruses (both observed in mj and not), then the association (v_i to m_j) is given a higher confidence score. Conversely, if vi is somewhat similar to coronaviruses observed in mj, and also somewhat similar to viruses not known to infect this particular mammal, then the association (v_i to m_j) is given a medium confidence score. The association (v_i to m_j) is given a lower confidence score if vi is more similar to coronaviruses not known to infect mj than it is similar to coronaviruses observed in this host.
Formally, given an adjacency matrix A of dimensions (left| {mathbf{V}} right| times left| {mathbf{M}} right|) where (left| {mathbf{V}} right|) is number of coronaviruses included in this study (for which a complete genome could be found), and (left| {mathbf{M}} right|) is number of included mammals, such that for each (v_i in {mathbf{V}}) and (m_j in {mathbf{M}}), aij = 1 if an association is known to exist between the virus and the mammal, and 0 otherwise. Then for a similarity matrix simviral corresponding to each of the similarity matrices calculated above, a learner from the viral point of view is defined as follows24,25:

$${mathrm{confidence}}_{{mathrm{viral}}}( {v_i to m_j} ) = frac{{mathop {sum}nolimits_{l = 1,,l ne i}^{left| {mathbf{V}} right|} {( {{mathrm{sim}}_{{mathrm{viral}}}( {v_i,,v_l} ) times a_{lj}} )} }}{{mathop {sum}nolimits_{l = 1,,l ne i}^{left| {mathbf{V}} right|} {{mathrm{sim}}_{{mathrm{viral}}}left( {v_i,,v_l} right)} }}$$
(3)

Mammals—the host point of view
We constructed seven learners from the similarities summarised in Table 3. Each of these learners calculated for every coronavirus–mammalian association (( {v_i to m_j})) a confidence score based on how similar the mammalian species mj is to known hosts of the coronavirus vi, compared to how similar mj is to mammals not associated with vi. For instance, if mj is phylogenetically close to known hosts of vi, and also phylogenetically distant to mammalian species not known to be associated with this coronavirus, then the phylogenetic similarly learner will assign (v_i to m_j) a higher confidence score. However, if mj does not overlap geographically with known hosts of vi, then the geographical overlap learner will assign it a low (in effect 0) confidence score.
Formally, given the above-defined adjacency matrix A, and a similarity matrix simmammalian corresponding to each of the similarity matrices summarised in Table 3, a learner from the mammalian point of view is defined as follows24,25:

$${mathrm{confidence}}_{{mathrm{mammalian}}}( {v_i to m_j} ) = frac{{mathop {sum}nolimits_{l = 1,,l ne j}^{left| {mathbf{M}} right|} {( {{mathrm{sim}}_{{mathrm{mammalian}}}( {m_j,,m_l} ) times a_{il}} )} }}{{mathop {sum}nolimits_{l = 1,,l ne j}^{left| {mathbf{M}} right|} {{mathrm{sim}}_{{mathrm{mammalian}}}( {m_j,,m_l} )} }}$$
(4)

Network—the network point of view
We integrated two learners based on network similarities—one for mammals and one for coronaviruses. Formally, given the adjacency matrix A, our two learners from the network point of view as defined as follows24:

$${mathrm{confidence}}_{{mathrm{network}}_{mathbf{V}}}( {v_i to m_j} ) = frac{{mathop {sum}nolimits_{l = 1,,l ne i}^{left| {mathbf{V}} right|} {( {{mathrm{sim}}_{{mathrm{network}}}( {v_i,,v_l} ) times a_{lj}} )} }}{{mathop {sum}nolimits_{l = 1,,l ne i}^{left| {mathbf{V}} right|} {{mathrm{sim}}_{{mathrm{network}}}( {v_i,,v_l} )} }};;$$
(5)

$${mathrm{confidence}}_{{mathrm{network}}_{mathbf{M}}}( {v_i to m_j} ) = frac{{mathop {sum}nolimits_{l = 1,,l ne j}^{left| {mathbf{M}} right|} {( {{mathrm{sim}}_{{mathrm{network}}}( {m_j,,m_l} ) times a_{il}} )} }}{{mathop {sum}nolimits_{l = 1,,l ne j}^{left| {mathbf{M}} right|} {{mathrm{sim}}_{{mathrm{network}}}( {m_j,,m_l} )} }}$$
(6)

Ensemble construction
We combined the learners described above by stacking them into ensembles (meta-ensembles) using Stochastic Gradient Boosting (GBM). The purpose of this combination is to incorporate the three points of views, as well as varied aspects of the coronaviruses and their mammalian potential hosts, into a generalisable, robust model50. We selected GBM as our stacking algorithm following an assessment of seven machine-learning algorithms using held-out test sets (20% of known associations randomly selected, N = 5—Supplementary Fig. 14). In addition, GBM is known for its ability to handle non-linearity and high-order interactions between constituent learners51, and have been used to predict reservoirs of viruses46 and zoonotic hot-spots51.We performed the training and optimisation (tuning) of these ensembles using the caret R Package52.
Sampling
Our GBM ensembles comprised 100 replicate models. Each model was trained with balanced random samples using tenfold cross-validation (Supplementary Fig. 4). Final ensembles were generated by taking mean predictions (probability) of constituent models. Predictions were calculated form the mean probability at three cut-offs: >0.5 (standard), >0.75 and ≥0.9821. SD from mean probability was also generated and its values subtracted/added to predictions, to illustrate variation in the underlying replicate models.
Validation and performance estimation
We validated the performance of our analytical pipeline externally against 20 held-out test sets. Each test set was generated by splitting the set of observed associations between coronaviruses and their hosts into two random sets: a training set comprising 85% of all known associations and a test set comprising 15% of known associations. These test sets were held-out throughout the processes of generating similarity matrices; similarity learning, and assembling our learners, and were only used for the purposes of estimating performance metrics of our analytical pipeline. This resulted in 20 runs in which our ensemble learnt using only 85% of observed associations between our coronaviruses and their mammalian hosts. For each run, we calculated three performance metrics based on the mean probability across each set of 100 replicate models of the GBM meta-ensembles: AUC, true skill statistics (TSS) and F-score.
AUC is a threshold-independent measure of model predictive performance that is commonly used as a validation metric for host–pathogen predictive models21,46. Use of AUC has been criticised for its insensitivity to absolute predicted probability and its inclusion of a priori untenable prediction51,53, and so we also calculated the TSS (TSS = sensitivity + specificity − 1)54. F-score captures the harmonic mean of the precision and recall and is often used with uneven class distribution. Our approach is relaxed with respect to false positives (unobserved associations), hence the low F-score recorded overall.
We selected three probability cut-offs for our meta-ensemble: 0.50, 0.75 and 0.9821. One extreme of our cut-off range (0.5) maximises the ability of our ensemble to detect known associations (higher AUC, lower F-score). The other (0.9821) is calculated so that 90% of known positives are captured by our ensemble, while reducing the number of additional associations predicted (higher F-score, lower AUC).
Changes in network structure
We quantified the diversity of the mammalian hosts of each coronavirus in our input by computing mean phylogenetic distance between these hosts. Similarly, we captured the diversity of coronaviruses associated with each mammalian species by calculating mean (hamming) distance between the genomes of these coronaviruses. We termed these two metrics: mammalian diversity per virus and viral diversity per mammal, respectively. We aggregated both metrics at the network level by means of simple average. This enabled us to quantify changes in these diversity metrics, at the level of network, with addition of predicted links at three probability cut-offs: >0.5, >0.75 and ≥0.9821.
In addition, we captured changes in the structure of the bipartite network linking CoVs with their mammalian hosts, with the addition of predicted associations, by computing a comprehensive set of structural properties (Supplementary Note 3) at the probability cut-offs mentioned above, and comparing the results with our original network. Here we ignore properties that deterministically change with the addition of links (e.g., degree centrality, connectance; Supplementary Table 2 lists all computed metrics and changes in their values). Instead, we focus on non-trivial structural properties. Specifically, we capture changes in network stability, by measuring its nestedness55,56,57; and we quantify non-independence in interaction patterns by means of C-score58. Supplementary Note 3 provides full definition of these concepts as well as other metrics we computed for our networks.
Reporting summary
Further information on research design is available in the Nature Research Reporting Summary linked to this article. More

1.
Carré, P. & Pouzet, A. Rapeseed market, worldwide and in Europe. OCL 21(1), D102. https://doi.org/10.1051/ocl/201h3054 (2014).
Article  Google Scholar 
2.
Hammond, K. E. & Lewis, B. E. The timing and sequence of events leading to stem canker disease in populations of Brassica napus var. oleifera in the field. Plant Pathol. 35, 551–556. https://doi.org/10.1111/j.1365-3059.1986.tb02054.x (1986).
Article  Google Scholar 

3.
Deb, D., Khan, A. & Dey, N. Phoma diseases: Epidemiology and control. Plant. Pathol. 00, 1–15. https://doi.org/10.1111/ppa.13221 (2020).
CAS  Article  Google Scholar 

4.
Fitt, B. D. L., Brun, H., Barbetti, M. J. & Rimmer, S. R. World-wide importance of Phoma stem canker (Leptosphaeria maculans and L. biglobosa) on oilseed rape (Brassica napus). Eur. J. Plant Pathol. 114, 3–15. https://doi.org/10.1007/s10658-005-2233-5 (2006).
Article  Google Scholar 

5.
Winter, M. & Koopmann, B. Race spectra of Leptosphaeria maculans, the causal agent of blackleg disease of oilseed rape, in different geographic regions in northern Germany. Eur. J. Plant Pathol. 145, 629–641. https://doi.org/10.1007/s10658-016-0932-8 (2016).
Article  Google Scholar 

6.
Derbyshire, M. C. & Denton-Giles, M. The control of Sclerotinia stem rot on oilseed rape (Brassica napus): current practices and future opportunities. Plant. Pathol. 65, 859–877. https://doi.org/10.1111/ppa.12517 (2016).
CAS  Article  Google Scholar 

7.
Gladders, P., Symonds, B. V., Hardwick, N. V. & Sansford, C. E. Opportunities to control canker (Leptosphaeria maculans) in winter oilseed rape by improved spray timing. IOBC/WPRS Bull. 21, 111–120 (1998).
Google Scholar 

8.
Kuai, J. et al. The effect of nitrogen application and planting density on the radiation use efficiency and the stem lignin metabolism in rapeseed (Brassica napus L.). Field Crops Res. 199, 89–98. https://doi.org/10.1016/j.fcr.2016.09.025 (2016).
Article  Google Scholar 

9.
Card, S. D. et al. Beneficial endophytic microorganisms of Brassica —A review. Biol. Control 90, 102–112. https://doi.org/10.1016/j.biocontrol.2015.06.001 (2015).
Article  Google Scholar 

10.
Weyens, N., van der Lelie, D., Taghavi, S., Newman, L. & Vangronsveld, J. Exploiting plant–microbe partnerships to improve biomass production and remediation. Trends Biotechnol. 27, 591–598. https://doi.org/10.1016/j.tibtech.2009.07.006 (2009).
CAS  Article  PubMed  Google Scholar 

11.
Müller, H. & Berg, G. Impact of formulation procedures on the effect of the biocontrol agent Serratia plymuthica HRO-C48 on Verticillium wilt in oilseed rape. Biocontrol 53, 905–916. https://doi.org/10.1007/s10526-007-9111-3 (2008).
Article  Google Scholar 

12.
Granér, G., Persson, P., Meijer, J. & Alström, S. A study on microbial diversity in different cultivars of Brassica napus in relation to its wilt pathogen, Verticillium longisporum. FEMS Microbiol. Lett. 224, 269–276. https://doi.org/10.1016/S0378-1097(03)00449-X (2003).
CAS  Article  PubMed  Google Scholar 

13.
Croes, S. et al. Bacterial communities associated with Brassica napus L. grown on trace-element-contaminated and non-contaminated fields: a genotypic and phenotypic comparison. Microb. Biotechnol. 6, 371–384. https://doi.org/10.1111/1751-7915.12057 (2013).
CAS  Article  PubMed  PubMed Central  Google Scholar 

14.
Zhang, Q. et al. Diversity and biocontrol potential of endophytic fungi in Brassica napus. Biol. Control 72, 98–102. https://doi.org/10.1016/j.biocontrol.2014.02.018 (2014).
Article  Google Scholar 

15.
Berg, G. et al. The rhizosphere effect on bacteria antagonistic towards the pathogenic fungus Verticillium differs depending on plant species and site. FEMS Microbiol. Ecol. 56, 250–261. https://doi.org/10.1111/j.1574-6941.2005.00025.x (2006).
CAS  Article  PubMed  Google Scholar 

16.
Berg, G. et al. Impact of plant species and site on rhizosphere-associated fungi antagonistic to Verticillium dahliae Kleb. Appl. Environ. Microbiol. 71, 4203–4213. https://doi.org/10.1128/AEM.71.8.4203-4213.2005 (2005).
CAS  Article  PubMed  PubMed Central  Google Scholar 

17.
Robin, A. H. K. et al. Leptosphaeria maculans alters glucosinolate profiles in blackleg disease-resistant and -susceptible cabbage lines. Front. Plant Sci. 8, 1789. https://doi.org/10.3389/fpls.2017.01769 (2017).
Article  Google Scholar 

18.
Garrido-Sanz, D. et al. Genomic and genetic diversity within the Pseudomonas fluorescens complex. PLoS ONE 11(2), e0150183. https://doi.org/10.1371/journal.pone.0153733 (2016).
CAS  Article  PubMed  PubMed Central  Google Scholar 

19.
Taylor, A. Fungal diversity in ectotomycorrhizal communities: sampling effort and species distribution. Plant Soil 244, 19–28. https://doi.org/10.1023/A:1020279815472 (2002).
ADS  CAS  Article  Google Scholar 

20.
Schmidt, C. S. et al. Distinct communities of poplar endophytes on an unpolluted and a risk elements-polluted site and their plant growth promoting potential in vitro. Microb. Ecol. 75, 955–969. https://doi.org/10.1007/s00248-017-1103-y (2018).
CAS  Article  PubMed  Google Scholar 

21.
Jedryczka, M. Epidemiology and damage caused by stem canker of oilseed rape in Poland. Phytopathol. Pol. 45, 73–75 (2007).
Article  Google Scholar 

22.
Mazáková, J., Urban, J., Zouhar, M. & Ryšánek, P. Analysis of Leptosphaeria species complex causing Phoma leaf spot and stem canker of winter oilseed rape (Brassica napus) in the Czech Republic. Crop Pasture Sci. 68, 254–264. https://doi.org/10.1071/CP16308 (2017).
CAS  Article  Google Scholar 

23.
El Hadrami, A., Fernando, W. G. D. & Daayf, F. Variations in relative humidity modulate Leptosphaeria spp. pathogenicity and interfere with canola mechanisms of defence. Eur. J. Plant Pathol. 126, 187–202. https://doi.org/10.1007/s10658-009-9532-1 (2010).
Article  Google Scholar 

24.
Hilton, S., Bennett, A. J., Chandler, D., Mills, P. & Bending, G. D. Preceding crop and seasonal effects influence fungal, bacterial and nematode diversity in wheat and oilseed rape rhizosphere and soil. Appl. Soil Ecol. 126, 34–46. https://doi.org/10.1016/j.apsoil.2018.02.007 (2018).
Article  Google Scholar 

25.
Glynou, K. et al. The local environment determines the assembly of root endophytic fungi at a continental scale. Environ. Microbiol. 18, 2418–2434. https://doi.org/10.1111/1462-2920.13112 (2016).
CAS  Article  PubMed  Google Scholar 

26.
Croes, S., Weyens, N., Colpaet, J. & Vangronveld, J. Characterization of the cultivable bacterial populations associated with field grown Brassica napus L.: An evaluation of sampling and isolation protocols. Environ. Microbiol. 17, 2379–2392., https://doi.org/10.1111/1462-2920.12701 (2015).

27.
Alström, S. Characteristics of bacteria from oilseed rape in relation to their biocontrol activity against Verticillium dahliae. J. Phytopathol. 149, 57–64. https://doi.org/10.1046/j.1439-0434.2001.00585.x (2001).
Article  Google Scholar 

28.
Cope-Selby, N. et al. Endophytic bacteria in Miscanthus seed: Implications for germination, vertical inheritance of endophytes, plant evolution and breeding. GCB Bioenergy 9, 57–77. https://doi.org/10.1111/gcbb.12364 (2017).
CAS  Article  Google Scholar 

29.
Rathore, R. et al. Crop establishment practices are a driver of the plant microbiota in winter oilseed rape (Brassica napus). Front. Microbiol. 8, 1489. https://doi.org/10.3389/fmicb.2017.01489 (2017).
Article  PubMed  PubMed Central  Google Scholar 

30.
Lay, C. Y. et al. Canola-Root-Associated microbiomes in the Canadian prairies. Front. Microbiol. 9, 1189. https://doi.org/10.3389/fmicb.2018.01188 (2018).
Article  Google Scholar 

31.
Sundara-Rao, W. V. B. & Sinha, M. K. Phosphate dissolving microorganisms in the soil and rhizosphere. Indian J. Agric. Sci. 33, 272–278. https://doi.org/10.1007/BF01372637 (1963).
Article  Google Scholar 

32.
Bashan, Y., Kamnev, A. A. & de-Bashan, L. E. Tricalcium phosphate is inappropriate as a universal selection factor for isolating and testing phosphate-solubilizing bacteria that enhance plant growth: A proposal for an alternative procedure. Biol. Fertil. Soils 49, 465–479. https://doi.org/10.1007/s00374-012-0737-7 (2013).
CAS  Article  Google Scholar 

33.
Pii, Y. et al. Microbial interactions in the rhizosphere: beneficial influences of plant growth-promoting rhizobacteria on nutrient acquisition process. A review. Biol. Fertil. Soils 51, 403–415. https://doi.org/10.1007/s00374-015-0996-1 (2015).
CAS  Article  Google Scholar 

34.
Reddy, C. A. & Saravanan, R. S. Polymicrobial multi-functional approach for enhancement of crop productivity. in Advances in Applied Microbiology (eds. Gadd, G. M. & Sariaslani, S.) 53–113 (Oxford Academic, Oxford, 2013).

35.
Lally, R. D. et al. Application of endophytic Pseudomonas fluorescens and a bacterial consortium to Brassica napus can increase plant height and biomass under greenhouse and field conditions. Front. Plant Sci. 8, 2193. https://doi.org/10.3389/fpls.2017.02193 (2017).
Article  PubMed  PubMed Central  Google Scholar 

36.
Parikh, L., Eskelson, M. J. & Adesemoye, A. O. Relationship of in vitro and in planta screening: improving the selection process for biological control agents against Fusarium root rot in row crops. Arch. Phytopathol. Plant Protect. 51, 156–169. https://doi.org/10.1080/03235408.2018.1441098 (2018).
Article  Google Scholar 

37.
Bakker, P. A. H. M., Pieterse, C. M. J. & van Loon, L. C. Induced systemic resistance by fluorescent Pseudomonas sp. Phytopathology 97, 239–243. https://doi.org/10.1094/PHYTO-97-2-0239 (2007).
Article  PubMed  Google Scholar 

38.
Youssef, S. A., Tartoura, K. A. & Greash, A. G. Serratia proteamaculans mediated alteration of tomato defense system and growth parameters in response to early blight pathogen Alternaria solani infection. Physiol. Mol. Plant Pathol. 103, 16–22. https://doi.org/10.1016/j.pmpp.2018.04.004 (2018).
CAS  Article  Google Scholar 

39.
Li, H. et al. The use of Pseudomonas fluorescens P13 to control Sclerotinia stem rot (Sclerotinia sclerotiorum) of oilseed rape. J. Microbiol. 49, 884–889. https://doi.org/10.1007/s12275-011-1261-4 (2011).
Article  PubMed  Google Scholar 

40.
Smolińska, U. & Kowalska, B. Biological control of the soil-borne fungal pathogen Sclerotinia sclerotiorum—A review. J. Plant Pathol. 100, 1–12. https://doi.org/10.1007/s42161-018-0023-0 (2018).
Article  Google Scholar 

41.
Shaukat, M. F. Seed bio-priming with Serratia plymuthica HRO-C48 for the control of Verticillium longisporum and Phoma lingam in Brassica napus L. spp. oleifera. (PhD Dissertation, University of Uppsala, Sweden, 2013).

42.
Castellano-Hinojosa, A., Pérez-Tapia, V., Bedmar, E. J. & Santillana, N. Purple corn-associated rhizobacteria with potential for plant growth promotion. J. Appl. Microbiol. 124, 1254–1264. https://doi.org/10.1111/jam.13708 (2018).
CAS  Article  PubMed  Google Scholar 

43.
Li, L. et al. Synergistic plant–microbe interactions between endophytic bacterial communities and the medicinal plant Glycyrrhiza uralensis F. Antonie Van Leeuwenhoek 111, 1735–1748. https://doi.org/10.1007/s10482-018-1062-4 (2018).
Article  PubMed  Google Scholar 

44.
Barnawal, D., Bharti, N., Maji, D., Chanotiya, C. S. & Kalra, A. 1-Aminocyclopropane-1-carboxylic acid (ACC) deaminase-containing rhizobacteria protect Ocimum sanctum plants during waterlogging stress via reduced ethylene generation. Plant Physiol. Biochem. 58, 227–235. https://doi.org/10.1016/j.plaphy.2012.07.008 (2012).
CAS  Article  PubMed  Google Scholar 

45.
Egamberdieva, D., Wirth, S., Behrendt, U., Ahmad, P. & Berg, G. Antimicrobial activity of medicinal plants correlates with the proportion of antagonistic endophytes. Front. Microbiol. 8, 199. https://doi.org/10.3389/fmicb.2017.00199 (2017).
Article  PubMed  PubMed Central  Google Scholar 

46.
Joe, M. M. et al. Resistance responses of rice to rice blast fungus after seed treatment with the endophytic Achromobacter xylosoxidans AUM54 strains. Crop Protect. 42, 141–148. https://doi.org/10.1016/j.cropro.2012.07.006 (2012).
Article  Google Scholar 

47.
Bertrand, H. et al. Stimulation of the ionic transport system in Brassica napus by a plant growth-promoting rhizobacterium (Achromobacter sp.). Can. J. Microbiol. 46, 229–236 (2000).
CAS  Article  Google Scholar 

48.
Abuamsha, R., Salman, M. & Ehlers, R. U. Role of different additives on survival of Serratia plymuthica HRO-C48 on oilseed rape seeds and control of Phoma lingam. Br. Microbiol. Res. J. 4, 737–748 (2014).
Article  Google Scholar 

49.
Garrity, G. M., Winters, M. & Searles, D. B. Taxonomic outline of the prokaryotes. in Bergey’s Manual of Systematic Bacteriology, 2nd Edn, Release 1.0 (Springer, New York, 2001).

50.
Unterseher, M. & Schnittler, M. Dilution-to-extinction cultivation of leaf-inhabiting endophytic fungi in beech (Fagus sylvatica L.)—Different cultivation techniques influence fungal biodiversity assessment. Mycol. Res. 113, 645–654. https://doi.org/10.1016/j.mycres.2009.02.002 (2009).
Article  PubMed  Google Scholar 

51.
Zadok, J. C., Chang, T. T. & Konzak, A. A decimal code for the growth stages of cereals. Weed Res. 14, 415–421. https://doi.org/10.1111/j.1365-3180.1974.tb01084.x (1974).
Article  Google Scholar 

52.
Schmidt, C. S., Mrnka, L., Frantík, T., Lovecká, P. & Vosátka, M. Plant growth promotion of Miscanthus × giganteus by endophytic bacteria and fungi on non-polluted and polluted soils. World J. Microbiol. Biotechnol. 34, 48. https://doi.org/10.1007/s11274-018-2426-7 (2018).
CAS  Article  PubMed  Google Scholar 

53.
Koubek, J. et al. Whole-cell MALDI-TOF: Rapid screening method in environmental microbiology. Int. Biodeter. Biodegr. 69, 82–86. https://doi.org/10.1016/j.ibiod.2011.12.007 (2012).
CAS  Article  Google Scholar 

54.
Uhlik, O. et al. Matrix-assisted laser desorption ionization (MALDI)–time of flight mass spectrometry- and MALDI biotyper-based identification of cultured biphenyl-metabolizing bacteria from contaminated horseradish rhizosphere soil. Appl. Environ. Microb. 77, 6858–6866. https://doi.org/10.1128/AEM.05465-11 (2011).
CAS  Article  Google Scholar 

55.
Štorchová, H. et al. An improved method of DNA isolation from plants collected in the field and conserved in saturated NaCl/CTAB solution. Taxon 49, 79–84. https://doi.org/10.2307/1223934 (2000).
Article  Google Scholar 

56.
White, T. J., Bruns, T. D., Lee, S. & Taylor, J. Analysis of phylogenetic relationship by amplification and direct sequencing of ribosomal RNA genes. in PCR Protocols: A Guide to Methods and Applications (eds. Innis, M. A., Gelfand, D. H., Sninsky, J. J. & White, T. J.) 315–322 (Academic Press Inc., New York, 1990).

57.
Gardes, M. & Bruns, T. D. ITS primers with enhanced specificity for basidiomycetes—Application to the identification of mycorrhizae and rusts. Mol. Ecol. 2, 113–118. https://doi.org/10.1111/j.1365-294X.1993.tb00005.x (1993).
CAS  Article  PubMed  Google Scholar 

58.
McLaughlin, D. J., Hibbett, D. S., Lutzoni, F., Spatafora, J. W. & Vilgalys, R. The search for the fungal tree of life. Trends Microbiol. 11, 488–497. https://doi.org/10.1016/j.tim.2009.08.001 (2009).
CAS  Article  Google Scholar 

59.
Alexander, D. B. & Zuberer, D. A. Use of chrome azurol S reagents to evaluate siderophore production by rhizosphere bacteria. Biol. Fertil. Soils 12, 39–45. https://doi.org/10.1007/BF00369386 (1991).
CAS  Article  Google Scholar 

60.
Penrose, D. M. & Glick, B. R. Methods for isolating and characterizing ACC deaminase-containing plant growth-promoting rhizobacteria. Physiol. Plant. 118, 10–15. https://doi.org/10.1034/j.1399-3054.2003.00086.x (2003).
CAS  Article  PubMed  Google Scholar 

61.
Li, Z., Chang, S., Lin, L., Li, Y. & An, Q. A colorimetric assay of 1-aminocyclopropane-1-carboxylate (ACC) based on ninhydrin reaction for rapid screening of bacteria containing ACC deaminase. Lett. Appl. Microbiol. 53, 178–185. https://doi.org/10.1111/j.1472-765X.2011.03088.x (2011).
CAS  Article  PubMed  Google Scholar 

62.
Villano, D., Fernandez-Pachon, M. S., Moya, M. L., Troncoso, A. M. & Garcıa-Parrilla, M. C. Radical scavenging ability of polyphenolic compounds towards DPPH free radical. Talanta 71, 230–235. https://doi.org/10.1016/j.talanta.2006.03.050 (2007).
CAS  Article  PubMed  Google Scholar 

63.
Hajšlová, J., Fenclová, M. & Zachariašová, M. Methodology for the Rapid Screening of Isolates of Endophytic Microorganisms and Identification of Strains with Phytohormonal Activity (in Czech, ISBN 978-80-7080-869-6 ) (2013).

64.
Veprikova, Z. et al. Mycotoxins in plant-based dietary supplements: Hidden health risk for consumers. J. Agric. Food Chem. 63, 6633–6643. https://doi.org/10.1021/acs.jafc.5b02105 (2015).
CAS  Article  PubMed  Google Scholar 

65.
Zhou, Q. Untersuchungen zum Infektionsmodus, immunologischen Nachweis und zur biologischen Bekämpfung von Leptosphaeria maculans (Desm) Ces. & de Not., dem Erreger der Wurzelhals- und Stängelfäule an Winterraps (Brassica napus L.). (Ph.D Dissertation, University of Göttingen, Göttingen, 2001).

66.
Chèvre, A. M. et al. Stabilization of resistance to Leptosphaeria maculans in Brassica napus–B. juncea recombinant lines and its introgression into spring-type Brassica napus. Plant Dis. 92, 1208–1214. https://doi.org/10.1094/PDIS-92-8-1208 (2008).
Article  PubMed  Google Scholar 

67.
El-Tarabily, K. A. et al. Biological control of Sclerotinia minor using a chitinolytic bacterium and actinomycetes. Plant Pathol. 49, 573–583. https://doi.org/10.1046/j.1365-3059.2000.00494.x (2000).
Article  Google Scholar 

68.
Clarke, K. R. & Warwick, R. M. Change in Marine Communities: An Approach to Statistical Analysis and Interpretation 2nd edn. (Primer-E, Plymouth, 2001).
Google Scholar 

69.
Frisvad, J. C., Smedsgaard, J., Larsen, T. O. & Samson, R. A. Mycotoxins, drugs and other extrolites produced by species in Penicillium subgenus Penicillium. Stud. Mycol. 49, 201–241 (2004).
Google Scholar 

70.
Romero, F. M., Rossi, F. R., Gárriz, A., Carrasco, P. & Ruíz, O. A. A bacterial endophyte from apoplast fluids protects canola plants from different pathogens via antibiosis and induction of host resistance. Phytopathology 109, 375–383 (2019).
CAS  Article  Google Scholar 

71.
Kamal, M. M., Lindbeck, K. D., Savocchia, S. & Ash, G. J. Biological control of Sclerotinia stem rot of canola using antagonistic bacteria. Plant Pathol. 64, 1375–1384 (2015).
CAS  Article  Google Scholar 

72.
Fernando, W. G. D., Nakkeeran, S., Zhang, Y., Savchuk, S. Biological control of Sclerotinia sclerotiorum (Lib.) de Bary by Pseudomonas and Bacillus species on canola petals. Crop Protect. 26, 100–107. https://doi.org/10.1016/j.cropro.2006.04.007 (2007)

73.
Peng, G., McGregor, L., Lahlali, R., Gossen, B. D., Hwang, S. F., Adhikari, K. K., Strelkov, S. E., McDonald, M. R. Potential biological control of clubroot on canola and crucifer vegetable crops. Plant Pathol. 60, 566–574. https://doi.org/10.1111/j.1365-3059.2010.02400.x (2011)

74.
Wu, Y., Yuan, J., Raza, W., Shen, Q., Huang, Q. Biocontrol traits and antagonistic potential of Bacillus amyloliquefaciens strain NJZJSB3 against Sclerotinia sclerotiorum, a causal agent of canola stem rot. J. Microbiol. Biotechnol. 24, 1327–1336. https://doi.org/10.4014/jmb.1402.02061 (2014)

75.
Auer, S. & Ludwig-Müller, J. Biological control of clubroot (Plasmodiophora brassicae) by an endophytic fungus. Integrated control in oilseed crops. IOBC-WPRS Bull. 136, 155–156 (2018).
Google Scholar 

76.
Huang, H.-C. & Erickson, R. S. Biological control of Sclerotinia stem rot of canola using Ulocladium atrum. Plant Pathol. Bull. 16, 55–59 (2007).
CAS  Google Scholar 

77.
Marques, A. P. G. C., Pires, C., Moreira, H., Rangel, A. O. S. S., Castro, P.M.L. Assessment of the plant growth promoting abilities of six bacterial isolates using Zea mays as indicator plant. Soil Biol. Biochem. 42, 1229–1235. https://doi.org/10.1016/j.soilbio.2010.04.014 (2010) More

1.
Agren, G. I. Stoichiometry and nutrition of plant growth in natural communities. Annu. Rev. Ecol. Evol. Syst. 39, 153–170 (2008).
Article  Google Scholar 
2.
Strand, A. E., Pritchard, S. G., McCormack, M. L., Davis, M. A. & Oren, R. Irreconcilable differences: Fine-root life spans and soil carbon persistence. Science 319, 456–458 (2008).
ADS  CAS  PubMed  Article  Google Scholar 

3.
Norby, R. J. & Jackson, R. B. Root dynamics and global change: Seeking an ecosystem perspective. New Phytol. 147, 3–12 (2000).
CAS  Article  Google Scholar 

4.
Valliere, J. M. & Allen, E. B. Interactive effects of nitrogen deposition and drought-stress on plant-soil feedbacks of Artemisia californica seedlings. Plant Soil 403, 277–290 (2016).
CAS  Article  Google Scholar 

5.
Schulte-Uebbing, L. & de Vries, W. Global-scale impacts of nitrogen deposition on tree carbon sequestration in tropical, temperate, and boreal forests: A meta-analysis. Glob. Change Biol. 24, E416–E431 (2018).
Article  Google Scholar 

6.
Vanguelova, E. I. & Pitman, R. M. Nutrient and carbon cycling along nitrogen deposition gradients in broadleaf and conifer forest stands in the east of England. For. Ecol. Manage. 447, 180–194 (2019).
Article  Google Scholar 

7.
Wang, L. X., Mou, P. P. & Jones, R. H. Nutrient foraging via physiological and morphological plasticity in three plant species. Can. J. For. Res. 36, 164–173 (2006).
CAS  Article  Google Scholar 

8.
Yu, G. et al. Stabilization of atmospheric nitrogen deposition in China over the past decade. Nat. Geosci. 12, 424–429 (2019).
ADS  CAS  Article  Google Scholar 

9.
Wang, W. J., Mo, Q. F., Han, X. G., Hui, D. F. & Shen, W. J. Fine root dynamics responses to nitrogen addition depend on root order, soil layer, and experimental duration in a subtropical forest. Biol. Fertil. Soils 55, 723–736 (2019).
CAS  Article  Google Scholar 

10.
Ostonen, I. et al. Fine root foraging strategies in Norway spruce forests across a European climate gradient. Glob. Change Biol. 17, 3620–3632 (2011).
ADS  Article  Google Scholar 

11.
Makita, N. et al. Fine root morphological traits determine variation in root respiration of Quercus serrata. Tree Physiol. 29, 579–585 (2009).
CAS  PubMed  Article  Google Scholar 

12.
Craine, J. M., Froehle, J., Tilman, G. D., Wedin, D. A. & Chapin, F. S. The relationships among root and leaf traits of 76 grassland species and relative abundance along fertility and disturbance gradients. Oikos 93, 274–285 (2001).
Article  Google Scholar 

13.
Liu, R. Q. et al. Plasticity of fine-root functional traits in the litter layer in response to nitrogen addition in a subtropical forest plantation. Plant Soil 415, 317–330 (2017).
CAS  Article  Google Scholar 

14.
Nadelhoffer, K. J. The potential effects of nitrogen deposition on fine-root production in forest ecosystems. New Phytol. 147, 131–139 (2000).
CAS  Article  Google Scholar 

15.
Noguchi, K., Nagakura, J. & Kaneko, S. Biomass and morphology of fine roots of sugi (Cryptomeria japonica) after 3 years of nitrogen fertilization. Front. Plant Sci. 4, 7 (2013).
Article  Google Scholar 

16.
Lu, X. K., Mao, Q. G., Gilliam, F. S., Luo, Y. Q. & Mo, J. M. Nitrogen deposition contributes to soil acidification in tropical ecosystems. Glob. Change Biol. 20, 3790–3801 (2014).
ADS  Article  Google Scholar 

17.
Li, W. B. et al. The effects of simulated nitrogen deposition on plant root traits: A meta-analysis. Soil Biol. Biochem. 82, 112–118 (2015).
CAS  Article  Google Scholar 

18.
Comas, L. H. & Eissenstat, D. M. Patterns in root trait variation among 25 co-existing North American forest species. New Phytol. 182, 919–928 (2009).
CAS  PubMed  Article  Google Scholar 

19.
Zhang, X. et al. Effects of long-term nitrogen addition and decreased precipitation on the fine root morphology and anatomy of the main tree species in a temperate forest. For. Ecol. Manag. 455, 117664 (2020). https://doi.org/10.1016/j.foreco.2019.117664.
Article  Google Scholar 

20.
Burton, A. J., Pregitzer, K. S. & Hendrick, R. L. Relationships between fine root dynamics and nitrogen availability in Michigan northern hardwood forests. Oecologia 125, 389–399 (2000).
ADS  CAS  PubMed  Article  Google Scholar 

21.
Pregitzer, K. S. et al. Fine root architecture of nine North American trees. Ecol. Monogr. 72, 293–309 (2002).
Article  Google Scholar 

22.
Wurzburger, N. & Wright, S. J. Fine-root responses to fertilization reveal multiple nutrient limitation in a lowland tropical forest. Ecology 96, 2137–2146 (2015).
PubMed  Article  Google Scholar 

23.
Penuelas, J., Sardans, J., Rivas-Ubach, A. & Janssens, I. A. The human-induced imbalance between C, N and P in Earth’s life system. Glob. Change Biol. 18, 3–6 (2012).
ADS  Article  Google Scholar 

24.
Loewe, A., Einig, W., Shi, L., Dizengremel, P. & Hampp, R. Mycorrhiza formation and elevated CO2 both increase the capacity for sucrose synthesis in source leaves of spruce and aspen. New Phytol. 145, 565–574 (2000).
CAS  Article  Google Scholar 

25.
Grechi, I. et al. Effect of light and nitrogen supply on internal C:N balance and control of root-to-shoot biomass allocation in grapevine. Environ. Exp. Bot. 59, 139–149 (2007).
CAS  Article  Google Scholar 

26.
Jing, H. et al. Effect of nitrogen addition on the decomposition and release of compounds from fine roots with different diameters: The importance of initial substrate chemistry. Plant Soil 438, 281–296 (2019).
CAS  Article  Google Scholar 

27.
Mucha, J. et al. Fine root classification matters: Nutrient levels in different functional categories, orders and diameters of roots in boreal Pinus sylvestris across a latitudinal gradient. Plant Soil 447, 507–520 (2019). https://doi.org/10.1007/s11104-019-04395-1.
CAS  Article  Google Scholar 

28.
Aubrey, D. P. & Teskey, R. O. Stored root carbohydrates can maintain root respiration for extended periods. New Phytol. 218, 142–152 (2018).
CAS  PubMed  Article  Google Scholar 

29.
Kou, L. et al. Simulated nitrogen deposition affects stoichiometry of multiple elements in resource-acquiring plant organs in a seasonally dry subtropical forest. Sci. Total Environ. 624, 611–620 (2018). https://doi.org/10.1016/j.scitotenv.2017.12.080.
ADS  CAS  Article  PubMed  Google Scholar 

30.
Yan, X. L., Jia, L. M. & Dai, T. F. Fine root morphology and growth in response to nitrogen addition through drip fertigation in a Populus × euramericana “Guariento” plantation over multiple years. Ann. For. Sci. 76 (2019).

31.
Yan, G. et al. Spatial and temporal effects of nitrogen addition on root morphology and growth in a boreal forest. Geoderma 303, 178–187 (2017).
ADS  CAS  Article  Google Scholar 

32.
Lu, X. K. et al. Plant acclimation to long-term high nitrogen deposition in an N-rich tropical forest. Proc. Natl. Acad. Sci. USA 115, 5187–5192 (2018).
CAS  PubMed  Article  Google Scholar 

33.
Van der Sande, M. T. et al. Soil fertility and species traits, but not diversity, drive productivity and biomass stocks in a Guyanese tropical rainforest. Funct. Ecol. 32, 461–474 (2018).
Article  Google Scholar 

34.
Burton, A. J., Jarvey, J. C., Jarvi, M. P., Zak, D. R. & Pregitzer, K. S. Chronic N deposition alters root respiration-tissue N relationship in northern hardwood forests. Glob. Change Biol. 18, 258–266 (2012).
ADS  Article  Google Scholar 

35.
Eissenstat, D. M. & Yanai, R. D. The ecology of root lifespan. Adv. Ecol. Res. 27, 2–60 (1997).
Google Scholar 

36.
Chen, D. M., Lan, Z. C., Hu, S. J. & Bai, Y. F. Effects of nitrogen enrichment on belowground communities in grassland: Relative role of soil nitrogen availability vs soil acidification. Soil Biol. Biochem. 89, 99–108 (2015).
CAS  Article  Google Scholar 

37.
Vanguelova, E. I., Nortcliff, S., Moffat, A. J. & Kennedy, F. Morphology, biomass and nutrient status of fine roots of Scots pine (Pinussylvestris) as influenced by seasonal fluctuations in soil moisture and soil solution chemistry. Plant Soil 270, 233–247 (2005).
CAS  Article  Google Scholar 

38.
Zhang, H., Liu, Y., Zhou, Z. & Zhang, Y. Inorganic nitrogen addition affects soil respiration and belowground organic carbon fraction for a Pinus tabuliformis forest. Forests 10, 369 (2019).
Article  Google Scholar 

39.
Lalnunzira, C., Brearley, F. Q. & Tripathi, S. K. Root growth dynamics during recovery of tropical mountain forest in North-east India. J. Mt. Sci. 16, 2335–2347 (2019).
Article  Google Scholar 

40.
Kochsiek, A., Tan, S. & Russo, S. E. Fine root dynamics in relation to nutrients in oligotrophic Bornean rain forest soils. Plant Ecol. 214, 869–882 (2013).
Article  Google Scholar 

41.
Ostonen, I. et al. Specific root length as an indicator of environmental change. Plant Biosyst. 141, 426–442 (2007).
Article  Google Scholar 

42.
Fitter, A. H., Stickland, T. R., Harvey, M. L. & Wilson, G. W. Architectural analysis of plant root systems 1. Architectural correlates of exploitation efficiency. New Phytol. 118, 375–382 (1991).
Article  Google Scholar 

43.
Hodge, A. The plastic plant: Root responses to heterogeneous supplies of nutrients. New Phytol. 162, 9–24 (2004).
Article  Google Scholar 

44.
Zhou, Y. M., Tang, J. W., Melillo, J. M., Butler, S. & Mohan, J. E. Root standing crop and chemistry after six years of soil warming in a temperate forest. Tree Physiol. 31, 707–717 (2011).
PubMed  Article  CAS  Google Scholar 

45.
Fujita, Y., Robroek, B. J. M., de Ruiter, P. C., Heil, G. W. & Wassen, M. J. Increased N affects P uptake of eight grassland species: The role of root surface phosphatase activity. Oikos 119, 1665–1673 (2010).
CAS  Article  Google Scholar 

46.
Alvarez-Clare, S. & Mack, M. C. Do foliar, litter, and root nitrogen and phosphorus concentrations reflect nutrient limitation in a lowland tropical wet forest?. PLoS ONE 10, e0123796 (2015).
PubMed  PubMed Central  Article  CAS  Google Scholar 

47.
Koerselman, W. & Meuleman, A. F. M. The vegetation NP ratio a new tool to detect the nature. J. Appl. Ecol. 33, 1441 (1996).
Article  Google Scholar 

48.
Wang, Z. Q. et al. The scaling of fine root nitrogen versus phosphorus in terrestrial plants: A global synthesis. Funct. Ecol. 33, 2081–2094 (2019).
Article  Google Scholar 

49.
Desrochers, A., Landhausser, S. M. & Lieffers, V. J. Coarse and fine root respiration in aspen (Populus tremuloides). Tree Physiol. 22, 725–732 (2002).
PubMed  Article  Google Scholar 

50.
Son, Y. & Hwang, J. H. Fine root biomass, production and turnover in a fertilized Larix leptolepis plantation in central Korea. Ecol. Res. 18, 339–346 (2003).
Article  Google Scholar 

51.
Eissenstat, D. M. & Volder, A. The efficiency of nutrient acquisition over the life of a root, 185–220. In Nutrient Acquisition by Plants: An Ecological Perspective (ed. Barririrad, H.) (Springer, Berlin, Heidelberg, 2005).
Google Scholar 

52.
Chen, L., Deng, Q., Yuan, Z., Mu, X. & Kallenbach, R. L. Age-related C:N:P stoichiometry in two plantation forests in the Loess Plateau of China. Ecol. Eng. 120, 14–22 (2018).
Article  Google Scholar 

53.
Wapongnungsang, R. H. & Tripathi, S. K. Fine root growth and soil nutrient dynamics during shifting cultivation in tropical semi-evergreen forests of northeast India. J. Environ. Biol. 40, 45–52 (2019).
CAS  Article  Google Scholar 

54.
Razaq, M., Salahuddin, Shen, H., Sher, H. & Zhang, P. Influence of biochar and nitrogen on fine root morphology, physiology, and chemistry of Acer mono. Sci. Rep. 7, 5367 (2017). https://doi.org/10.1038/s41598-017-05721-2.

55.
Kobe, R. K., Iyer, M. & Walters, M. B. Optimal partitioning theory revisited: Nonstructural carbohydrates dominate root mass responses to nitrogen. Ecology 91, 166–179 (2010).
PubMed  Article  Google Scholar 

56.
Li, W. B. et al. Effects of nitrogen enrichment on tree carbon allocation: A global synthesis. Glob. Ecol. Biogeogr. 29, 573–589 (2020).
Article  Google Scholar 

57.
Wang, T. et al. Age structure of Picea schrenkiana forest along an altitudinal gradient in the central Tianshan Mountains, northwestern China. For. Ecol. Manage. 196, 267–274 (2004).
Article  Google Scholar 

58.
Li, J., Yutao, Z., Li, J., Li, X. & Lu, J. Effect of stimulated nitrogen deposition on the fine root decomposition and related nutrient release of Picea schrenkiana var. tianshanica. Acta Bot. Boreal. Occident. Sin. 35, 0182–0188 (2015) (in chinese).
Google Scholar 

59.
Bremner, J. & Mulvaney, R. Urease Activity in Soils (Academic Press, London, 1978).
Google Scholar 

60.
Liu, Y. et al. Nitrogen addition alleviates microbial nitrogen limitations and promotes soil respiration in a subalpine coniferous forest. Forests 10, 16 (2019).
Google Scholar 

61.
Bao, S. N. Soil Agrochemical Analysis 20–38 (China Agricultural Press, Beijing, 2000) (in Chinese).
Google Scholar 

62.
Buysse, J. & Merckx, R. An improved colorimetric method to quantify sugar content of plant tissue. J. Exp. Bot. 44, 1627–1629 (1993).
CAS  Article  Google Scholar 

63.
Su, L. et al. Soil and fine roots ecological stoichiometry in different vegetation restoration stages in a karst area, southwest China. J. Environ. Manag. 252, 109694 (2019).
CAS  Article  Google Scholar  More

We used a dataset from a grassland diversity experiment at Zürich-Reckenholz, Switzerland, in the Atlantic central climatic zone of Europe. The data contain measurements on many functions from 78 plots that comprised monocultures and mixtures sown at a wide range of species relative abundances, set up at three levels of N fertiliser application and maintained for 3 years following establishment, which is a typical time in grassland-crop rotations.
Monocultures and mixtures were sown following a simplex design64. Four perennial species, known to be key forage species in ruminant production, were selected based on the factorial combination of their functional traits related to temporal establishment (fast-establishing vs. temporally persistent), and N acquisition (non-fixing for grasses, N2-fixing for legumes). The species were Lolium perenne L. cultivar (cv.) Lacerta (fast-establishing grass), Dactylis glomerata L. cv. Accord (temporally persistent grass), Trifolium pratense L. cv. Merviot (fast-establishing legume), and Trifolium repens L. cv. Milo (temporally persistent legume). The type of stands were: monocultures (100% of one species), binary mixtures (50% of each of two species), an equi-proportional mixture (25% of each of the four species), dominant mixtures (70% of the dominant species, 10% of each of the other three), and co-dominant mixtures (40% of each of two species, 10% of each of the other two; see Supplementary Table S1). All types of stands were sown at two levels of overall sown density, with the high level being the recommended seed weight (100%) under conditions typical of Switzerland, and the low level being 60%.
The experiment was sown in August 2002 on plots of 3 m × 6 m and was maintained from 2003 (year 1) to 2005 (year 3). The plots were fertilised with N fertiliser (as NH4NO3) at rates following a geometric series: 50, 150, or 450 kg N ha−1 yr−1 (N50, N150, and N450, respectively), split into five equal applications. In early spring, all plots received phosphorus and potassium in amounts expected to be non-limiting for intensively managed grasslands on fertile soils in Switzerland. At the N150 treatment, all types of monocultures and mixtures were established, whereas the N50 and N450 treatments only included the monocultures, the equi-proportional mixture, and the dominant mixtures. The 78 plots were arranged in a fully randomised design. Consult Nyfeler et al.37 for full details of the experimental design, establishment, and maintenance.
Ten functions were measured representing (i) forage production: aboveground biomass yield, standard deviation of yield, temporal stability, weed biomass; (ii) N cycling: symbiotic N2 fixation, N efficiency, NO3 in soil solution; and (iii) forage quality: crude protein content, organic matter digestibility, metabolisable energy content (Table 1). To date, detailed analyses from the experiment have been published on two functions, namely biomass yield37 and symbiotic N2 fixation33.
Measurement of functions
Aboveground biomass yield and weed biomass
All plots were harvested five times annually at 5 cm above ground surface. Aboveground biomass yield at each harvest was determined by drying a representative subsample to constant weight (65° C for 48 h), and this data was summed to give total annual biomass yield. Biomass proportions of the four sown and pooled unsown species (weeds) were measured by manually separating samples from permanent sub-plots (0.8 m × 0.3 m), which was done at the first, third, and fifth harvest of each year. These data allowed for calculation of weed biomass per ha and year.
Standard deviation and stability of yield
Year-to-year standard deviation of yield (SDyield) was calculated from the annual yields of the three experimental years, and stability was defined as the ratio of averaged annual yields to year-to-year SDyield (following Lehman and Tilman65). To measure yield variation within each year, seasonal SDyield was calculated from the five annual harvests, and seasonal stability was defined as the ratio of total annual yield to seasonal SDyield. We purposely use both SDyield and stability as both measures are essential to evaluate yield variation66.
Symbiotic N 2fixation
Symbiotic N2 fixation (Nsym) was determined by the isotope dilution method67. Double-labelled 15N-enriched 15NH415NO3 was applied on a permanently defined, central part of each plot (1.4 m × 1.5 m). Plant samples were analysed for 15N and 14N abundance by gas isotope ratio mass spectrometry and by thermal conductometry. Nsym in the sward, as calculated here, comprises legume N derived from the atmosphere (Ndfa) plus N derived from apparent Ndfa transfer to the grass (Ntrans). See Supplementary Appendix S1 and Nyfeler et al.33 for full details of measurements and calculations.
N efficiency
N efficiency was defined as the ratio of total N yield to the amount of applied fertiliser N and therefore measures the total N output of the system in relation to the fertiliser N input. Total N yield was calculated by first multiplying N content from biomass samples with their total dry mass to give the N yield per harvest. Annual total N yield was then computed as the sum of all harvests.
NO 3in soil solution (NO3)
Porous cup tension lysimeters were installed to extract soil water from a depth of 60 cm below ground surface. In 2-week intervals from October 2004 to April 2006, a suction of 80 kPa was applied 1 day prior to sampling, and concentrations of nitrate–N (NO3-N) were determined by spectrophotometry. We note that NO3 data were only available for years 2 and 3. See Supplementary Appendix S1 for details of the measurements.
Crude protein content 
Crude protein content (CP) in stand biomass was calculated from the N content in biomass samples, multiplied by 6.25. The justification for the multiplicative factor is given by the fact that all biological proteins contain on average 16% N68.
Organic matter digestibility
Organic matter digestibility (OM digestibility) was determined from biomass samples of the second and fourth harvest following the two-stage in vitro fermentation process with rumen liquor and acidic pepsin solution according to Tilley and Terry69; see Supplementary Appendix S1 for details. Information on OM digestibility was only available for years 2 and 3 of the study.
Metabolisable energy content
Metabolisable energy content (ME) of stand biomass was calculated based on OM digestibility and CP following a reference manual of Agroscope70; see Supplementary Appendix S1 for calculation. Due to the connection with the measurement OM digestibility, ME data were only available for years 2 and 3.
Data for each function were computed at the plot level for each of three experimental years (the three exceptions as noted). For analyses across years, data was averaged across available years, except SDyield and stability (see above).
Data analyses
We applied the multivariate modelling framework13 to estimate simultaneously species identity and diversity effects of the ten functions along with effects of N fertilisation. To allow direct comparisons of the model terms, all functions’ data were standardised to a common scale by dividing them by their maximum value (at a single year) over the 3-year experiment and N fertilisation treatments. This scaling allowed for a direct comparison of results among years. Note that the multivariate approach is a generalisation of the univariate diversity interaction model61, and we refer to Supplementary Appendix S1 for a summary to the univariate regression.
In the following, we generally refer to the analysis of data averaged across experimental years, and all equations model the response at a single plot (plot subscripts are omitted). A preliminary regression equation was specified for the kth function (k = 1–10) with:

$${y}_{k}={alpha }_{k}mathrm{DENS}+sum_{f=1}^{3}sum_{i=1}^{4}{beta }_{ifk}{P}_{i}times {mathrm{N}_mathrm{Treat}}_{f}+sum_{begin{array}{c}i,j=1\ i More

arising from: P. B. Ngor et al.; Scientific Reports https://doi.org/10.1038/s41598-018-27340-1 (2018).
As one of the richest sources of fisheries-related data in the lower Mekong basin, the Tonle Sap dai fishery has received considerable attention in the literature in recent years as concerns grow over the impacts of hydropower dams on fisheries, which are important for livelihoods and food security1,2,3.
Ngor et al.4 reported a decline since 2000 in the catch of larger species which tend to occupy higher trophic levels; compensatory increases in the catch of smaller species; and declines in the mean body weight (and length) of common species in the Tonle Sap dai fishery, as evidence of the effects of indiscriminate fishing or “fishing-down” of the multi-species fish assemblage in the lower Mekong basin. We provide evidence below that suggest that these apparent recent changes are more likely to reflect changing hydrological conditions than fishing-down effects, possibly caused by climate change and recently also by hydropower development.
The dai fishery has been reliably monitored since 1997–98. Without explanation, Ngor et al. excluded the first three seasons (1997–98 to 1999–2000) of monitoring data which include one of the driest fishing seasons on record (1998–99). The authors thereby created a time series beginning with the three wettest seasons (largest floods) since monitoring began (2000–1 to 2002–3) that were followed by 12 seasons of variable, but decreasing flows caused by hydropower dam construction, low rainfalls possibly resulting from climate change, and abstractions for agriculture5,6 (Fig. 1).
Figure 1

Source: Mekong River Commission Secretariat.

The flood index (FI) or flood pulse14 in the Tonle Sap Great Lake System (1997/08–2014/15). The FI is a measure of the flood extent and duration, calculated as the sum of the flooded area days above the mean flooded area from April to March of the following year2. Whilst highly variable, a downward decline (p-value = 0.06) in the FI is observed between 2000/01 (Year 2001) and 2014/15 (Year 2015) shown by solid circles. Adding the most recent data for 2016–2018 (not shown here), confirmed that a downward linear trend in the FI since the 2000/01 season is statistically significant (p-value  45 cm) excluding those with zero catch in any year. These 28 species formed approximately 16% of the total catch during the study period. We also found negative regression coefficients for all 28 species, supporting the findings of Ngor et al. However, the combined annual catch of these 28 species did not decline significantly through time (R2 = 0.22; p-value = 0.07).
We did however find that the combined annual catch of these 28 larger species varied significantly with the annual flood index (FI)—a measure of flood extent and duration (R2 = 0.46; p-value  45 cm) species and the flood index (R2 = 0.46; p-value  More

1.
Paterson, S. et al. Antagonistic coevolution accelerates molecular evolution. Nature 464, 275–278 (2010).
CAS  PubMed  PubMed Central  Article  Google Scholar 
2.
Schulte, R. D., Makus, C., Hasert, B., Michiels, N. K. & Schulenburg, H. Multiple reciprocal adaptations and rapid genetic change upon experimental coevolution of an animal host and its microbial parasite. Proc. Natl Acad. Sci. USA 107, 7359–7364 (2010).
CAS  PubMed  Article  PubMed Central  Google Scholar 

3.
Koskella, B. & Lively, C. M. Evidence for negative frequency-dependent selection during experimental coevolution of a freshwater snail and a sterilizing trematode. Evolution 63, 2213–2221 (2009).
PubMed  Article  PubMed Central  Google Scholar 

4.
Decaestecker, E. et al. Host–parasite ‘Red Queen’ dynamics archived in pond sediment. Nature 450, 870–873 (2007).
CAS  PubMed  Article  PubMed Central  Google Scholar 

5.
Gómez, P. & Buckling, A. Bacteria–phage antagonistic coevolution in soil. Science 332, 106–109 (2011).
PubMed  Article  CAS  PubMed Central  Google Scholar 

6.
Refardt, D. & Ebert, D. Inference of parasite local adaptation using two different fitness components. J. Evol. Biol. 20, 921–929 (2007).
CAS  PubMed  Article  PubMed Central  Google Scholar 

7.
Duffy, M. A., Hall, S. R., Cáceres, C. E. & Ives, A. R. Rapid evolution, seasonality, and the termination of parasite epidemics. Ecology 90, 1441–1448 (2009).
PubMed  Article  Google Scholar 

8.
Springer, Y. P. Clinical resistance structure and pathogen local adaptation in a serpentine flax–flax rust interaction. Evolution 61, 1812–1822 (2007).
PubMed  Article  PubMed Central  Google Scholar 

9.
Tack, A. J. M., Laine, A.-L., Burdon, J. J., Bissett, A. & Thrall, P. H. Below-ground abiotic and biotic heterogeneity shapes above-ground infection outcomes and spatial divergence in a host–parasite interaction. New Phytol. 207, 1159–1169 (2015).
PubMed  PubMed Central  Article  Google Scholar 

10.
Wolinska, J. & King, K. C. Environment can alter selection in host–parasite interactions. Trends Parasitol. 25, 236–244 (2009).
PubMed  Article  PubMed Central  Google Scholar 

11.
Auld, S. K. J. R., Hall, S. R., Ochs, J. H., Sebastian, M. & Duffy, M. A. Predators and patterns of within-host growth can mediate both among-host competition and evolution of transmission potential of parasites. Am. Nat. 184, S77–S90 (2014).
PubMed  Article  PubMed Central  Google Scholar 

12.
Wright, R. C. T., Brockhurst, M. A. & Harrison, E. Ecological conditions determine extinction risk in co-evolving bacteria–phage populations. BMC Evol. Biol. 16, 227 (2016).
PubMed  PubMed Central  Article  CAS  Google Scholar 

13.
Duffy, M. A. et al. Ecological context influences epidemic size and parasite-driven evolution. Science 335, 1636–1638 (2012).
CAS  PubMed  Article  Google Scholar 

14.
Auld, S. K. J. R. & Brand, J. Environmental variation causes different (co) evolutionary routes to the same adaptive destination across parasite populations. Evol. Lett. 1, 245–254 (2017).
PubMed  PubMed Central  Article  Google Scholar 

15.
Su, M. & Boots, M. The impact of resource quality on the evolution of virulence in spatially heterogeneous environments. J. Theor. Biol. 416, 1–7 (2017).
PubMed  Article  PubMed Central  Google Scholar 

16.
Auld, S. K. J. R. & Tinsley, M. C. The evolutionary ecology of complex lifecycle parasites: linking phenomena with mechanisms. Heredity 114, 125–132 (2015).
CAS  PubMed  Article  PubMed Central  Google Scholar 

17.
Cardon, M., Loot, G., Grenouillet, G. & Blanchet, S. Host characteristics and environmental factors differentially drive the burden and pathogenicity of an ectoparasite: a multilevel causal analysis. J. Anim. Ecol. 80, 657–667 (2011).
PubMed  Article  PubMed Central  Google Scholar 

18.
Mahmud, M. A., Bradley, J. E. & MacColl, A. D. C. Abiotic environmental variation drives virulence evolution in a fish host–parasite geographic mosaic. Funct. Ecol. 31, 2138–2146 (2017).
Article  Google Scholar 

19.
Arruda, J. A., Marzolf, G. R. & Faulk, R. T. The role of suspended sediments in the nutrition of zooplankton in turbid reservoirs. Ecology 64, 1225–1235 (1983).
Article  Google Scholar 

20.
Mostowy, R. & Engelstädter, J. The impact of environmental change on host–parasite coevolutionary dynamics. Proc. R. Soc. B 278, 2283–2292 (2011).
PubMed  Article  Google Scholar 

21.
Thompson, J. N. The Geographic Mosaic of Coevolution (Univ. Chicago Press, 2005).

22.
Brett, M. T. Chaoborus and fish-mediated influences on Daphnia longispina population structure, dynamics and life history strategies. Oecologia 89, 69–77 (1992).
PubMed  Article  Google Scholar 

23.
Goss, L. B. & Bunting, D. L. Daphnia development and reproduction: responses to temperature. J. Therm. Biol. 8, 375–380 (1983).
Article  Google Scholar 

24.
Luijckx, P., Fienberg, H., Duneau, D. & Ebert, D. A matching-allele model explains host resistance to parasites. Curr. Biol. 23, 1085–1088 (2013).
CAS  PubMed  Article  Google Scholar 

25.
Bento, G. et al. The genetic basis of resistance and matching-allele interactions of a host–parasite system: the Daphnia magna–Pasteuria ramosa model. PLoS Genet. 13, e1006596 (2017).
PubMed  PubMed Central  Article  CAS  Google Scholar 

26.
Grosberg, R. K. Mate selection and the evolution of highly polymorphic self/nonself recognition genes. Science 289, 2111–2114 (2000).
CAS  PubMed  Article  Google Scholar 

27.
Hutchinson, G. E. The Ecological Theater and the Evolutionary Play (Yale Univ. Press, 1965).

28.
Stuart, Y. E. et al. Contrasting effects of environment and genetics generate a continuum of parallel evolution. Nat. Ecol. Evol. 1, 0158 (2017).
Article  Google Scholar 

29.
Klüttgen, B., Dülmer, U., Engels, M. & Ratte, H. ADaM, an artificial freshwater for the culture of zooplankton. Water Res. 28, 743–746 (1994).
Article  Google Scholar 

30.
Ebert, D., Zschokke-Rohringer, C. D. & Carius, H. J. Within- and between-population variation for resistance of Daphnia magna to the bacterial endoparasite Pasteuria ramosa. Proc. R. Soc. B 265, 2127–2134 (1998).
Article  Google Scholar 

31.
Auld, S. K. J. R. & Brand, J. Simulated climate change, epidemic size, and host evolution across host–parasite populations. Glob. Change Biol. 23, 5045–5053 (2017).
Article  Google Scholar 

32.
Holm, S. A simple sequentially rejective multiple test procedure. Scand. J. Stat. 6, 65–70 (1979).
Google Scholar 

33.
R Core Team R: A Language and Environment for Statistical Computing (R Foundation for Statistical Computing, 2019).

34.
Brereton, R. G. & Lloyd, G. R. Re-evaluating the role of the Mahalanobis distance measure. J. Chemom. 30, 134–143 (2016).
CAS  Article  Google Scholar 

35.
D’Orazio, M. StatMatch: Statistical Matching or Data Fusion. R package version 1.4.0 (2019).

36.
Goslee, S. C. & Urban, D. L. The ecodist package for dissimilarity-based analysis of ecological data. J. Stat. Softw. 22, 1–22 (2007).
Article  Google Scholar 

37.
Lefcheck, J. S. piecewiseSEM: piecewise structural equation modelling in R for ecology, evolution and systematics. Methods Ecol. Evol. 7, 573–579 (2016).
Article  Google Scholar 

38.
Auld, S. K. J. R., Wilson, P. J. & Little, T. J. Rapid change in parasite infection traits over the course of an epidemic in a wild host–parasite population. Oikos 123, 232–238 (2014).
Article  Google Scholar 

39.
Shocket, M. S. et al. Parasite rearing and infection temperatures jointly influence disease transmission and shape seasonality of epidemics. Ecology 99, 1975–1987 (2018).
PubMed  Article  PubMed Central  Google Scholar 

40.
Duncan, A. B., Mitchell, S. E. & Little, T. J. Parasite-mediated selection and the role of sex and diapause in Daphnia. J. Evol. Biol. 19, 1183–1189 (2006).
CAS  PubMed  Article  Google Scholar 

41.
Auld, S. K. J. R. et al. Variation in costs of parasite resistance among natural host populations. J. Evol. Biol. 26, 2479–2486 (2013).
CAS  PubMed  Article  Google Scholar 

42.
Laine, A.-L. Evolution of host resistance: looking for coevolutionary hotspots at small spatial scales. Proc. R. Soc. B 273, 267–273 (2006).
PubMed  Article  Google Scholar 

43.
Lohse, K., Gutierrez, A. & Kaltz, O. Experimental evolution of resistance in Paramecium caudatum against the bacterial parasite Holospora undulata. Evolution 60, 1177–1186 (2006).
Article  Google Scholar 

44.
Duffy, M. A. & Sivars-Becker, L. Rapid evolution and ecological host–parasite dynamics. Ecol. Lett. 10, 44–53 (2007).
PubMed  Article  Google Scholar 

45.
Brewer, M. J., Butler, A. & Cooksley, S. L. The relative performance of AIC, AICC and BIC in the presence of unobserved heterogeneity. Methods Ecol. Evol. 7, 679–692 (2016).
Article  Google Scholar 

46.
Shipley, B. A new inferential test for path models based on directed acyclic graphs. Struct. Equ. Model. 7, 206–218 (2000).
Article  Google Scholar  More

1.
Mayr E. Populations, species, and evolution: an abridgment of animal species and evolution. Cambridge: Belknap Press of Harvard University Press; 1970.
2.
Bickford D, Lohman DJ, Sodhi NS, Ng PKL, Meier R, Winker K, et al. Cryptic species as a window on diversity and conservation. Trends Ecol Evol. 2007;22:148–55.
PubMed  Article  Google Scholar 

3.
Fišer C, Robinson CT, Malard F. Cryptic species as a window into the paradigm shift of the species concept. Mol Ecol. 2018;27:613–35.
PubMed  Article  Google Scholar 

4.
Struck TH, Feder JL, Bendiksby M, Birkeland S, Cerca J, Gusarov VI, et al. Finding evolutionary processes hidden in cryptic species. Trends Ecol Evol. 2018;33:153–63.
PubMed  Article  Google Scholar 

5.
Sarno D, Kooistra WHCF, Medlin LK, Percopo I, Zingone A. Diversity in the genus Skeletonema (Bacillariophyceae). II. An assessment of the taxonomy of S. costatum-like species with the description of four new species. J Phycol. 2005;41:151–76.
Article  Google Scholar 

6.
Gaonkar CC, Kooistra WHCF, Lange CB, Montresor M, Sarno D. Two new species in the Chaetoceros socialis complex (Bacillariophyta): C. sporotruncatus and C. dichatoensis, and characterization of its relatives. J Phycol. 2017;53:889–907.
CAS  PubMed  Article  Google Scholar 

7.
Li Y, Boonprakob A, Gaonkar CC, Kooistra WHCF, Lange CB, Hernández-Becerril D, et al. Diversity in the globally distributed diatom genus Chaetoceros (Bacillariophyceae): three new species from warm-temperate waters. PLoS ONE. 2017;12:e0168887.
PubMed  PubMed Central  Article  CAS  Google Scholar 

8.
Finlay BJ, Clarke KJ. Ubiquitous dispersal of microbial species. Nature. 1999;400:828.
CAS  Article  Google Scholar 

9.
Finlay BJ, Fenchel T. Divergent perspectives on protist species richness. Protist. 1999;150:229–33.
CAS  PubMed  Article  Google Scholar 

10.
Fenchel T, Finlay BJ. The ubiquity of small species: patterns of local and global diversity. Bioscience. 2004;54:777.
Article  Google Scholar 

11.
Fenchel T. Cosmopolitan microbes and their ‘cryptic’ species. Aquat Microb Ecol. 2005;41:49–54.
Article  Google Scholar 

12.
Miglietta MP, Faucci A, Santini F. Speciation in the sea: overview of the symposium and discussion of future directions. Integr Comp Biol. 2011;51:449–55.
PubMed  Article  Google Scholar 

13.
Kooistra WHCF, Sarno D, Balzano S, Gu H, Andersen RA, Zingone A. Global diversity and biogeography of Skeletonema species (Bacillariophyta). Protist. 2008;159:177–93.
CAS  PubMed  Article  Google Scholar 

14.
Nanjappa D, Audic S, Romac S, Kooistra WHCF, Zingone A. Assessment of species diversity and distribution of an ancient diatom lineage using a DNA metabarcoding approach. PLoS ONE. 2014;9:e103810.
PubMed  PubMed Central  Article  CAS  Google Scholar 

15.
Kaczmarska I, Mather L, Luddington IA, Muise F, Ehrman JM. Cryptic diversity in a cosmopolitan diatom known as Asterionellopsis glacialis (Fragilariaceae): Implications for ecology, biogeography, and taxonomy. Am J Bot. 2014;101:267–86.
PubMed  Article  Google Scholar 

16.
Zhao Y, Yi Z, Gentekaki E, Zhan A, Al-Farraj SA, Song W. Utility of combining morphological characters, nuclear and mitochondrial genes: An attempt to resolve the conflicts of species identification for ciliated protists. Mol Phylogenet Evol. 2016;94:718–29.
PubMed  Article  Google Scholar 

17.
Weiner A, Aurahs R, Kurasawa A, Kitazato H, Kucera M. Vertical niche partitioning between cryptic sibling species of a cosmopolitan marine planktonic protist. Mol Ecol. 2012;21:4063–73.
PubMed  Article  Google Scholar 

18.
Lamari N, Ruggiero MV, d’Ippolito G, Kooistra WHCF, Fontana A, Montresor M. Specificity of lipoxygenase pathways supports species delineation in the marine diatom genus Pseudo-nitzschia. PLoS ONE. 2013;8:e73281.
CAS  PubMed  PubMed Central  Article  Google Scholar 

19.
Škaloud P, Friedl T, Hallmann C, Beck A, Dal Grande F. Taxonomic revision and species delimitation of coccoid green algae currently assigned to the genus Dictyochloropsis (Trebouxiophyceae, Chlorophyta). J Phycol. 2016;52:599–617.
PubMed  Article  CAS  Google Scholar 

20.
de Jesus PB, Costa AL, de Castro Nunes JM, Manghisi A, Genovese G, Morabito M, et al. Species delimitation methods reveal cryptic diversity in the Hypnea cornuta complex (Cystocloniaceae, Rhodophyta). Eur J Phycol. 2019;54:135–53.
Article  CAS  Google Scholar 

21.
Díaz-Tapia P, Ly M, Verbruggen H. Extensive cryptic diversity in the widely distributed Polysiphonia scopulorum (Rhodomelaceae, Rhodophyta): molecular species delimitation and morphometric analyses. Mol Phylogenet Evol. 2020;152:106909.
PubMed  Article  Google Scholar 

22.
Huson DH, Rupp R, Scornavacca C. Phylogenetic networks. Cambridge: Cambridge University Press; 2009.

23.
Huson DH, Bryant D. Application of phylogenetic networks in evolutionary studies. Mol Biol Evol. 2006;23:254–67.
CAS  PubMed  Article  Google Scholar 

24.
Solís-Lemus C, Yang M, Ané C. Inconsistency of species tree methods under gene flow. Syst Biol. 2016;65:843–51.
PubMed  Article  Google Scholar 

25.
Deiner K, Bik HM, Mächler E, Seymour M, Lacoursière-Roussel A, Altermatt F, et al. Environmental DNA metabarcoding: Transforming how we survey animal and plant communities. Mol Ecol. 2017;26:5872–95.
PubMed  Article  Google Scholar 

26.
Pawlowski J, Audic S, Adl S, Bass D, Belbahri L, Berney C, et al. CBOL protist working group: barcoding eukaryotic richness beyond the animal, plant, and fungal kingdoms. PLoS Biol. 2012;10:e1001419.
CAS  PubMed  PubMed Central  Article  Google Scholar 

27.
Trobajo R, Mann DG, Clavero E, Evans KM, Vanormelingen P, McGregor RC. The use of partial cox1, rbcL and LSU rDNA sequences for phylogenetics and species identification within the Nitzschia palea species complex (Bacillariophyceae). Eur J Phycol. 2010;45:413–25.
CAS  Article  Google Scholar 

28.
Decelle J, Suzuki N, Mahé F, De Vargas C, Not F. Molecular phylogeny and morphological evolution of the acantharia (Radiolaria). Protist. 2012;163:435–50.
PubMed  Article  Google Scholar 

29.
Stoeck T, Przybos E, Dunthorn M. The D1-D2 region of the large subunit ribosomal DNA as barcode for ciliates. Mol Ecol Resour. 2014;14:458–68.
CAS  PubMed  Article  Google Scholar 

30.
Moniz MBJ, Kaczmarska I. Barcoding of diatoms: nuclear encoded ITS revisited. Protist. 2010;161:7–34.
CAS  PubMed  Article  Google Scholar 

31.
Gile GH, Stern RF, James ER, Keeling PJ. DNA barcoding of chlorarachniophytes using nucleomorph ITS sequences. J Phycol. 2010;46:743–50.
CAS  Article  Google Scholar 

32.
Stern RF, Andersen RA, Jameson I, Küpper FC, Coffroth M-A, Vaulot D, et al. Evaluating the ribosomal internal transcribed spacer (ITS) as a candidate dinoflagellate barcode marker. PLoS ONE. 2012;7:e42780.
CAS  PubMed  PubMed Central  Article  Google Scholar 

33.
Saunders GW. Applying DNA barcoding to red macroalgae: a preliminary appraisal holds promise for future applications. Philos Trans R Soc B Biol Sci. 2005;360:1879–88.

34.
MacGillivary ML, Kaczmarska I. Survey of the efficacy of a short fragment of the rbcL gene as a supplemental DNA barcode for diatoms. J Eukaryot Microbiol. 2011;58:529–36.
CAS  PubMed  Article  Google Scholar 

35.
Zimmermann J, Jahn R, Gemeinholzer B. Barcoding diatoms: evaluation of the V4 subregion on the 18S rRNA gene, including new primers and protocols. Org Divers Evol. 2011;11:173–92.
Article  Google Scholar 

36.
Piredda R, Tomasino MP, D’Erchia AM, Manzari C, Pesole G, Montresor M, et al. Diversity and temporal patterns of planktonic protist assemblages at a Mediterranean Long Term Ecological Research site. FEMS Microbiol Ecol. 2016;93:fiw200.
PubMed  Article  CAS  Google Scholar 

37.
Pawlowski J, Lecroq B. Short rDNA barcodes for species identification in foraminifera. J Eukaryot Microbiol. 2010;57:197–205.
CAS  PubMed  Article  Google Scholar 

38.
Mordret S, Piredda R, Vaulot D, Montresor M. Kooistra WHCF, Sarno D. dinoref: a curated dinoflagellate (Dinophyceae) reference database for the 18S rRNA gene. Mol Ecol Resour. 2018;18:974–87.
CAS  Article  Google Scholar 

39.
Gaonkar CC, Piredda R, Minucci C, Mann DG, Montresor M, Sarno D, et al. Annotated 18S and 28S rDNA reference sequences of taxa in the planktonic diatom family Chaetocerotaceae. PLoS ONE. 2018;13:e0208929.
PubMed  PubMed Central  Article  Google Scholar 

40.
Balzano S, Percopo I, Siano R, Gourvil P, Chanoine M, Marie D, et al. Morphological and genetic diversity of Beaufort Sea diatoms with high contributions from the Chaetoceros neogracilis species complex. J Phycol. 2017;53:161–87.
CAS  PubMed  Article  Google Scholar 

41.
Kopf A, Bicak M, Kottmann R, Schnetzer J, Kostadinov I, Lehmann K, et al. The ocean sampling day consortium. Gigascience. 2015;4. https://doi.org/10.1186/s13742-015-0066-5.

42.
Pesant S, Not F, Picheral M, Kandels-Lewis S, Le Bescot N, Gorsky G, et al. Open science resources for the discovery and analysis of Tara Oceans data. Sci Data. 2015;2:150023.
CAS  PubMed  PubMed Central  Article  Google Scholar 

43.
Yau S, Lopes dos Santos A, Eikrem W, Gérikas Ribeiro C, Gourvil P, Balzano S, et al. Mantoniella beaufortii and Mantoniella baffinensis sp. nov. (Mamiellales, Mamiellophyceae), two new green algal species from the high arctic. J Phycol. 2020;56:37–51.
PubMed  Article  Google Scholar 

44.
Lopes Dos Santos A, Gourvil P, Tragin M, Noël M-H, Decelle J, Romac S, et al. Diversity and oceanic distribution of prasinophytes clade VII, the dominant group of green algae in oceanic waters. ISME J. 2017;11:512–28.
PubMed  Article  Google Scholar 

45.
Kuwata A, Yamada K, Ichinomiya M, Yoshikawa S, Tragin M, Vaulot D, et al. Bolidophyceae, a sister picoplanktonic group of diatoms—a review. Front Mar Sci. 2018;5:370.
Article  Google Scholar 

46.
Segawa T, Matsuzaki R, Takeuchi N, Akiyoshi A, Navarro F, Sugiyama S, et al. Bipolar dispersal of red-snow algae. Nat Commun. 2018;9:1–8.
CAS  Article  Google Scholar 

47.
Ichinomiya M, Dos Santos AL, Gourvil P, Yoshikawa S, Kamiya M, Ohki K, et al. Diversity and oceanic distribution of the Parmales (Bolidophyceae), a picoplanktonic group closely related to diatoms. ISME J. 2016;10:2419–34.
PubMed  PubMed Central  Article  Google Scholar 

48.
Tragin M, Vaulot D. Novel diversity within marine Mamiellophyceae (Chlorophyta) unveiled by metabarcoding. Sci Rep. 2019;9:1–14.
CAS  Article  Google Scholar 

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

50.
Pinseel E, Janssens SB, Verleyen E, Vanormelingen P, Kohler TJ, Biersma EM, et al. Global radiation in a rare biosphere soil diatom. Nat Commun. 2020;11:1–12.
Article  CAS  Google Scholar 

51.
Hasle GR, Syvertsen EE. Marine diatoms. In: Tomas CR, editor. Identifying marine phytoplankton. San Diego: Academic Press; 1997. pp 5–385.

52.
Kooistra WHCF, Sarno D, Hernández-Becerril DU, Assmy P, Di Prisco C, Montresor M. Comparative molecular and morphological phylogenetic analyses of taxa in the Chaetocerotaceae (Bacillariophyta). Phycologia. 2010;49:471–500.
Article  Google Scholar 

53.
De Luca D, Sarno D, Piredda R, Kooistra WHCF. A multigene phylogeny to infer the evolutionary history of Chaetocerotaceae (Bacillariophyta). Mol Phylogenet Evol. 2019;140:106575.
PubMed  Article  Google Scholar 

54.
Longhurst AR. Toward and ecological geography of the sea. In: Longhurst AR, editor. Ecological geography of the sea. 2nd ed. Cambridge: Academic Press; 2007. pp 1–17.

55.
Stoeck T, Bass D, Nebel M, Christen R, Jones MDM, Breiner H-W, et al. Multiple marker parallel tag environmental DNA sequencing reveals a highly complex eukaryotic community in marine anoxic water. Mol Ecol. 2010;19:21–31.
CAS  PubMed  Article  Google Scholar 

56.
Amaral-Zettler LA, McCliment EA, Ducklow HW, Huse SM. A method for studying protistan diversity using massively parallel sequencing of V9 hypervariable regions of small-subunit ribosomal RNA genes. PLoS ONE. 2009;4:e6372.
PubMed  PubMed Central  Article  CAS  Google Scholar 

57.
Schloss PD, Westcott SL, Ryabin T, Hall JR, Hartmann M, Hollister EB, et al. Introducing mothur: open-source, platform-independent, community-supported software for describing and comparing microbial communities. Appl Environ Microbiol. 2009;75:7537–41.
CAS  PubMed  PubMed Central  Article  Google Scholar 

58.
De Vargas C, Audic S, Tara Oceans Consortium C, Tara Oceans Expedition P. Total V9 rDNA information organized at the metabarcode level for the Tara Oceans Expedition (2009–12). 2017. PANGAEA. https://doi.org/10.1594/PANGAEA.873277.

59.
Ibarbalz FM, Henry N, Brandão MC, Martini S, Busseni G, Byrne H, et al. Global trends in marine plankton diversity across kingdoms of life. Cell. 2019;179:1084–97.
CAS  PubMed  PubMed Central  Article  Google Scholar 

60.
Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ. Basic local alignment search tool. J Mol Biol. 1990;215:403–10.
CAS  Article  Google Scholar 

61.
Katoh K, Rozewicki J, Yamada KD. MAFFT online service: multiple sequence alignment, interactive sequence choice and visualization. Brief Bioinform. 2019;20:1160–6.
CAS  PubMed  PubMed Central  Article  Google Scholar 

62.
Price MN, Dehal PS, Arkin AP. FastTree 2—approximately maximum-likelihood trees for large alignments. PLoS ONE. 2010;5:e9490.
PubMed  PubMed Central  Article  CAS  Google Scholar 

63.
Han MV, Zmasek CM. phyloXML: XML for evolutionary biology and comparative genomics. BMC Bioinform. 2009;10:356.
Article  CAS  Google Scholar 

64.
Templeton AR, Crandall KA, Sing CF. A cladistic analysis of phenotypic associations with haplotypes inferred from restriction endonuclease mapping and DNA sequence data. III. Cladogram estimation. Genetics. 1992;132:619–33.
CAS  PubMed  PubMed Central  Article  Google Scholar 

65.
Clement M, Posada D, Crandall KA. TCS: a computer program to estimate gene genealogies. Mol Ecol. 2000;9:1657–9.
CAS  PubMed  Article  Google Scholar 

66.
Leigh JW, Bryant D. popart: full‐feature software for haplotype network construction. Methods Ecol Evol. 2015;6:1110–6.
Article  Google Scholar 

67.
Nguyen L-T, Schmidt HA, von Haeseler A, Minh BQ. IQ-TREE: a fast and effective stochastic algorithm for estimating maximum-likelihood phylogenies. Mol Biol Evol. 2015;32:268–74.
CAS  Article  Google Scholar 

68.
Kalyaanamoorthy S, Minh BQ, Wong TKF, von Haeseler A, Jermiin LS. ModelFinder: fast model selection for accurate phylogenetic estimates. Nat Methods. 2017;14:587–9.
CAS  PubMed  PubMed Central  Article  Google Scholar 

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

70.
Jukes TH, Cantor CR. Evolution of protein molecules. Mamm Protein Metab. 1969;3:21–132.
CAS  Article  Google Scholar 

71.
Meyer CP, Paulay G. DNA barcoding: error rates based on comprehensive sampling. PLoS Biol. 2005;3:e422.
PubMed  PubMed Central  Article  CAS  Google Scholar 

72.
R Core Team. R: a language and environment for statistical computing. 2019. Vienna, Austria: R Foundation for Statistical Computing; 2019.

73.
McMurdie PJ, Holmes S. phyloseq: an R package for reproducible interactive analysis and graphics of microbiome census data. PLoS ONE. 2013;8:e61217.
CAS  PubMed  PubMed Central  Article  Google Scholar 

74.
Wickham H. ggplot2: elegant graphics for data analysis. New York: Springer-Verlag; 2016. https://ggplot2.tidyverse.org.

75.
Becker A, Wilks AR. Maps: draw geographical maps. 2018. https://CRAN.R-project.org/package=maps.

76.
Markmann M, Tautz D. Reverse taxonomy: an approach towards determining the diversity of meiobenthic organisms based on ribosomal RNA signature sequences. Philos Trans R Soc Lond B Biol Sci. 2005;360:1917–24.
CAS  PubMed  PubMed Central  Article  Google Scholar 

77.
López-Escardó D, Paps J, de Vargas C, Massana R, Ruiz-Trillo I, Del Campo J. Metabarcoding analysis on European coastal samples reveals new molecular metazoan diversity. Sci Rep. 2018;8:9106.
PubMed  PubMed Central  Article  CAS  Google Scholar 

78.
Álvarez I, Wendel JF. Ribosomal ITS sequences and plant phylogenetic inference. Mol Phylogenet Evol. 2003;29:417–34.

79.
Alverson AJ, Kolnick L. Intragenomic nucleotide polymorphism among small subunit (18S) rDNA paralogs in the diatom genus Skeletonema (Bacillariophyta). J Phycol. 2005;41:1248–57.
CAS  Article  Google Scholar 

80.
Gaonkar CC, Piredda R, Sarno D, Zingone A, Montresor M, Kooistra WHCF. Species detection and delineation in the marine planktonic diatoms Chaetoceros and Bacteriastrum through metabarcoding: making biological sense of haplotype diversity. Environ Microbiol. 2020;22:1917–29.
CAS  PubMed  Article  Google Scholar 

81.
Cleve PT. Pelagisk Diatomeer från Kattegat. In: Petersen CGJ, editor. Det Videnskabelige Udbytte af Kanonbaaden ‘Hauchs’ Togter i de Danske Have Indefor Skagen, I. Aarene 1883–86. Kjøbenhavn: Andr. Fred. Høst & Sons Forlag; 1889. pp 53–56.

82.
Gran HH. Den Norske Nordhaus-Expedition 1876-1878. Botanik, Protophyta: Diatomaceae, Silicoflagellata og Cilioflagellata. Christiania: Grøndal & Søns; 1897.

83.
De Luca D, Kooistra WHCF, Sarno D, Gaonkar CC, Piredda R. Global distribution and diversity of Chaetoceros (Bacillariophyta, Mediophyceae): integration of classical and novel strategies. PeerJ. 2019;7:e7410.
PubMed  PubMed Central  Article  Google Scholar 

84.
Wang J, Wu J. Occurrence and potential risks of harmful algal blooms in the East China Sea. Sci Total Environ. 2009;407:4012–21.
CAS  PubMed  Article  Google Scholar 

85.
Zhen Y, Mi T, Yu Z. Detection of several harmful algal species by sandwich hybridization integrated with a nuclease protection assay. Harmful Algae. 2009;8:651–7.
CAS  Article  Google Scholar 

86.
Richter DJ, Watteaux R, Vannier T, Leconte J, Frémont P, Reygondeau G, et al. Genomic evidence for global ocean plankton biogeography shaped by large-scale current systems. bioRxiv. 2019. https://doi.org/10.1101/867739.

87.
Sarno D, Kooistra WHCF, Balzano S, Hargraves PE, Zingone A. Diversity in the genus Skeletonema (Bacillariophyceae). III. Phylogenetic position and morphological variability of Skeletonema costatum and Skeletonema grevillei, with the description of Skeletonema ardens sp. nov. J Phycol. 2007;43:156–70.
CAS  Article  Google Scholar 

88.
Hasle GR. The biogeography of some marine planktonic diatoms. Deep Sea Res Oceanogr Abstr. 1976;23:319–338, IN1-IN6.

89.
Pargana A. Functional and molecular diversity of the diatom family Leptocylindraceae. 2017. PhD Thesis, The Open University, Milton Keynes, UK.

90.
Novis PM. Taxonomy of Klebsormidium (Klebsormidiales, Charophyceae) in New Zealand streams and the significance of low-pH habitats. Phycologia. 2006;45:293–301.
Article  Google Scholar 

91.
Rindi F, Guiry MD, López-Bautista JM. Distribution, morphology, and phylogeny of Klebsormidium (Klebsormidiales, Charophyceae) in urban environments in Europe. J Phycol. 2008;44:1529–40.
PubMed  Article  Google Scholar 

92.
Rindi F, Mikhailyuk TI, Sluiman HJ, Friedl T, López-Bautista JM. Phylogenetic relationships in Interfilum and Klebsormidium (Klebsormidiophyceae, Streptophyta). Mol Phylogenet Evol. 2011;58:218–31.
PubMed  Article  Google Scholar 

93.
Škaloud P, Rindi F. Ecological differentiation of cryptic species within an asexual protist morphospecies: a case study of filamentous green alga Klebsormidium (Streptophyta). J Eukaryot Microbiol. 2013;60:350–62.
PubMed  Article  CAS  Google Scholar 

94.
Baas Becking LGM. Geobiologie of Inleiding tot de Milieukunde. The Hague: Van Stockum & Zoon; 1934).

95.
Shapiro BJ, Leducq J-B, Mallet J. What is speciation? PLoS Genet. 2016;12:e1005860.
PubMed  PubMed Central  Article  CAS  Google Scholar 

96.
Godhe A, Rynearson T. The role of intraspecific variation in the ecological and evolutionary success of diatoms in changing environments. Philos Trans R Soc Lond B Biol Sci. 2017;372:20160399.
PubMed  PubMed Central  Article  Google Scholar 

97.
de Vargas C, Norris R, Zaninetti L, Gibb SW, Pawlowski J. Molecular evidence of cryptic speciation in planktonic foraminifers and their relation to oceanic provinces. Proc Natl Acad Sci USA. 1999;96:2864–8.
PubMed  Article  Google Scholar 

98.
Amato A, Kooistra WHCF, Levialdi Ghiron JH, Mann DG, Pröschold T, Montresor M. Reproductive isolation among sympatric cryptic species in marine diatoms. Protist. 2007;158:193–207.
CAS  PubMed  Article  Google Scholar 

99.
Weisse T. Distribution and diversity of aquatic protists: an evolutionary and ecological perspective. Biodivers Conserv. 2007;17:243–59.
Article  Google Scholar 

100.
Vanelslander B, Créach V, Vanormelingen P, Ernst A, Chepurnov VA, Sahan E, et al. Ecological differentiation between sympatric pseudocryptic species in the estuarine benthic diatom Navicula phyllepta (Bacillariophyceae). J Phycol. 2009;45:1278–89.
CAS  PubMed  Article  Google Scholar  More