Ecology

This study uses SincNet according to the instructions provided by the authors for its application in a different dataset32. This section provides an introduction to SincNet and NIPS4Bplus before detailing the experimental procedure.SincNetThe first convolutional layer of a standard CNN trained on the raw waveform learns filters from the data, where each filter has a number of parameters that matches the filter length (Eq. 1).$$yleft[ n right] = xleft[ n right] times fleft[ n right] = mathop sum limits_{i = 0}^{I – 1} xleft[ i right] cdot fleft[ {n – i} right],$$
(1)

where (xleft[ n right]) is the chunk of the sound, (fleft[ n right]) is the filter of length (I), and (yleft[ n right]) is the filtered output. All the elements of the filter ((i)) are learnable parameters. SincNet replaces (fleft[ n right]) with another function (g) that only depends on two parameters per filter: the lower and upper frequencies of a rectangular bandpass filter (Eq. 2).$$gleft[ {n,f_{l} ,f_{h} } right] = 2f_{h} sincleft( {2pi f_{h} n} right) – 2f_{l} sincleft( {2pi f_{l} n} right),$$
(2)

where (f_{l} text{ and } f_{h}) are the learnable parameters corresponding to the low and high frequencies of the filter and (sincleft( x right) = frac{sinleft( x right)}{x}). The function (g) is smoothed with a Hamming window and the learnable parameters are initialised with given cut-off frequencies in the interval (left[ {0,frac{{f_{s} }}{2}} right]), where (f_{s}) is the sampling frequency.This first layer of SincNet performs the sinc-based convolutions for a set number and length of filters, over chunks of the raw waveform of given window size and overlap. A conventional CNN architecture follows the first layer, that in this study maintains the architecture and uses both standard and enhanced settings. The standard settings used are those of the TIMIT speaker recognition experiment27,32. They include two convolutional layers after the first layer with 60 filters of length 5. All three convolutions use layer normalisation. Next, three fully-connected (leaky ReLU) layers with 2048 neurons each follow, normalised with batch normalisation. To obtain frame-level classification, a final softmax output layer, using LogSoftmax, provides a set of posterior probabilities over the target classes. The classification for each file derives from averaging the frame predictions and voting for the class that maximises the average posterior. Training uses the RMSprop optimiser with the learning rate set to 0.001 and minibatches of size 128. A sample of sinc-based filters generated during this study shows their response both in the time and the frequency domains (Fig. 4).Figure 4Examples of learned SincNet filters. The top row (a–c) shows the filters in the time domain, the bottom row (d–f) shows their respective frequency response.Full size imageThe SincNet repository32 provides an alternative set of settings used in the Librispeech speaker recognition experiment27. Tests of the alternative settings, which include changes in the hidden CNN layers, provided similar results to those of the TIMIT settings and are included as Supplementary Information 1.NIPS4BplusNIPS4Bplus includes two parts: sound files and rich labels. The sound files are the training files of the 2013 NIPS4B challenge for bird song classification23. They are a single channel with a 44.1 kHz sampling rate and 32 bit depth. They comprise field recordings collected from central and southern France and southern Spain15. There are 687 individual files with lengths from 1 to 5 s for a total length of 48 min. The tags in NIPS4Bplus are based on the labels released with the 2013 Bird Challenge but annotated in detail by an experienced bird watcher using dedicated software15. The rich labels include the name of the species, the class of sound, the starting time and the duration of each sound event for each file. The species include 51 birds, 1 amphibian and 9 insects. For birds there can be two types of vocalisations: call and song; and there is also the drumming of a woodpecker. Calls are generally short sounds with simple patterns, while songs are usually longer with greater complexity and can have modular structures or produced by one of the sexes8,13. In the dataset, only bird species have more than one type of sound, with a maximum of two types. The labels in NIPS4Bplus use the same 87 tags present in the 2013 Bird Challenge training dataset with the addition of two other tags: “human” and “unknown” (for human sounds and calls which could not be identified). Tagged sound events in the labels typically correspond to individual syllables although in some occasions the reviewer included multiple syllables into single larger events15. The tags cover only 569 files of the original training set of 687 files. Files without tags include 100 that, for the purpose of the challenge, had no bird sounds but only background noise. Other files were excluded for different reasons such as vocalisations hard to identify or containing no bird or only insect sounds15. The 2013 Bird Challenge also includes a testing dataset with no labels that we did not use15.The total number of individual animal sounds tagged in the NIPS4Bplus labels is 5478. These correspond to 61 species and 87 classes (Fig. 5). The mean length of each tagged sound ranges from ~ 30 ms for Sylcan_call (the call of Sylvia cantillans, subalpine warbler) to more than 4.5 s for Plasab_song (the song of Platycleis sabulosa, sand bush-cricket). The total recording length for a species ranges from 0.7 s for Turphi_call (the call of Turdus philomelos, song thrush) to 51.4 s for Plasab_song. The number of individual files for each call type varies greatly from 9 for Cicatr_song (the call of Cicadatra atra, black cicada) to 282 for Sylcan_call.Figure 5Distribution of sound types by number of calls (number of files) and total length in seconds. Sound types are sorted first by taxonomic group and then by alphabetical order.Full size imageProcessing NIPS4BplusThe recommended pre-processing of human speech files for speaker recognition using SincNet includes the elimination of silent leading and trailing sections and the normalisation of the amplitude27. This study attempts to replicate this by extracting each individual sound as a new file according to the tags provided in the NIPS4Bplus labels. A Python script42 uses the content of the labels to read each wavefile, apply normalisation, select the time of origin and length specified in each individual tag and save it as a new wavefile. The name of the new file includes the original file name and a sequential number suffix according to the order in which tags are listed in the label files (the start time of the sound) to match the corresponding call tags at the time of processing. Each wavefile in the new set fully contains a sound according to the NIPS4Bplus labels. A cropped file may contain sounds from more than one species15, with over 20% of the files in the new set overlapping, at least in part, with sound from another species. The machine learning task does not use files containing background noise or the other parts of the files that are not tagged in the NIPS4Bplus labels. A separate Python script42 generates the lists of files and tags that SincNet requires for processing. The script randomly generates a 75:25 split into lists of train and test files and a Python dictionary key that assigns each file to the corresponding tag according to the file name. The script selects only files confirmed as animal sounds (excluding the tags “unknown” and “human”) and generates three different combinations of tags, as follows: (1) “All classes”: includes all the 87 types of tags originally included in the 2013 Bird Challenge training dataset; (2) “Bird classes”: excludes tags for insects and one amphibian species for a total of 77 classes; and (3) “Bird species”: one class for each bird species independently of the sounds type (call, songs and drumming are merged for each species) for a total of 51 classes. The script also excludes three very short files (length shorter than 10 ms) which could not be processed without code modifications.To facilitate the repeatability of the results, this study attempts to maintain the default parameters of SincNet used in the TIMIT speaker identification task27,32. The number and length of filters in the first sinc-based convolutional layer was set to the same values as the TIMIT experiment (80 filters of length 251 samples) as was the architecture of the CNN. The filters were initialised following the Mel scale cut-off frequencies. We did change the following parameters: (1) reduced the window frame size (cw_len) from 200 to 10 ms to accommodate for the short duration of some of the sounds in the NIPS4Bplus tags (such as some bird vocalisations); (2) reduced the window shift (cw_shift) from 10 to 1 ms in proportion to the reduction in window size (a value a 0.5 could not be given without code modifications); (3) updated the sampling frequency setting (fs) from the TIMIT 16,000 to the 44,100 Hz of the present dataset; and (4) updated the number of output classes (class_lay) to match the number of classes in each training run.To evaluate performance, the training sequence was repeated with the same settings and different random train and test file splits. Five training runs took place for each of the selection of tags: “All classes”, “Bird classes” and “Bird species”.Enhancements and comparisonsChanges in the parameters of SincNet result in different levels of performance. To assess possible improvements and provide baselines to compare against other models we attempted to improve the performance by adjusting a series of parameters, but did not modify the number of layers or make functional changes to the code other than the two outlined below. The parameters tested include: the length of the window frame size, the number and length of the filters in the first layer, number of filters and lengths of the other convolutional and fully connected layers, the length and types of normalisation in the normalisation layers, alternative activation and classification functions, and the inclusion of dropouts (Supplementary Information 1). In addition the SincNet code includes a hard-coded random amplification of each sound sequence; we also tested changing the level and excluding this random amplification through changes in the code. In order to process window frames larger than some of the labelled calls in the NIPS4Bplus dataset, the procedure outlined earlier in which files are cut according to the labels was replaced by a purpose-built process. The original files were not cut, instead a custom python script42 generated train and test file lists that contain the start and length of each labelled call. A modification of the SincNet code42 uses these lists to read the original files and select the labelled call. When the call is shorter than the window frame the code randomly includes the surrounding part of the file to complete the length of the window frame. Grid searches for individual parameters or combinations of similar parameters, over a set number of epochs, selected the best performing values. We also tested the use of the Additive Margin Softmax (AM-softmax) as a cost function37. The best performing models reported in the results use combinations of the best parameter values (Supplementary Information 1). All enhancements and model comparisons use the same dataset selection, that is the same train and test dataset split, of the normalised files for each set of tagged classes.The comparison using waveform + CNN models trained directly on the raw waveform, replaces the initial sinc-based convolution of SincNet with a standard 1d convolutional layer27, thus retaining the same network architecture as SincNet. As with SincNet enhancements, a series of parameter searches provided the best parameter combinations to obtain the best performing models.The pre-trained models used for comparison are DenseNet121, ResNet50 and VGG16 with architectures and weights sourced from the Torchvision library of PyTorch33. We tested three types of spectrograms: Fast Fourier Transform (FFT), Mel spectrum (Mel) and Mel-frequency cepstral coefficient (MFCC) to fine-tune the pre-trained models. FFT calculations used a frame length of 1024 samples, 896 samples overlap and a Hamming window. Mel spectrogram calculations used 128 Mel bands. Once normalised and scaled to 255 pixel intensity three repeats of the same spectrogram represented each of the three input channels of the pre-trained models. The length of sound used to generate the spectrograms was 3 s, and similarly with routines above, for labelled calls shorter than 3 s the spectrogram would randomly include the surrounding sounds. That is, the extract would randomly start in the interval between the end of the labelled call minus 3 s and the start of the call plus 3 s. This wholly includes the labelled call but its position is random within the 3 s sample. A fully connected layer replaced the final classifying layer of the pre-trained models to output the number of labelled classes. In the fine-tuning process the number of trainable layers of the model was not limited to the final fully connected layer, but also included an adjustable number of final layers to improve the results. The learning rate set initially to 0.0001 was halved if the validation loss stopped decreasing for 10 epochs.MetricsMeasures of performance include accuracy, ROC AUC, precision, recall, F1 score, top 3 accuracy and top 5 accuracy. Accuracy, calculated as part of the testing routine, is the ratio between the number of correctly predicted files of the test set and the total number of test files. The calculation of the other metrics uses the Scikit-learn module43 relying on the predicted values provided by the model and performing weighted averages. The ROC AUC calculation uses the mean of the posterior probabilities provided by SincNet for each tagged call. In the pre-trained models the ROC AUC calculations used the probabilities obtained after normalising the output with a softmax function. More

Sample collection and filtrationWe used 40 water samples for eDNA metabarcoding from 27 sites in 9 rivers and 13 lakes in Japan from 2016 to 2018 (Fig. 4). Sampling ID and detailed information for each site are listed in Supplementary Table S1. In the river water sampling, 1-L water samples were collected from the surface of at the shore of each river using bleached plastic bottles. In the field, a 1-ml Benzalkonium chloride solution (BAC, Osvan S, Nihon Pharmaceutical, Tokyo, Japan)33 was added to each water sample to suppress eDNA degeneration before filtering the water samples. We did not include field negative control samples in the HTS library, considering the aim of the presents study. The lake samples were provided by Doi et al. (2020)34 as DNA extracted samples. In the lake samples, 1-L water samples were collected from the surface at shore sites at each lake. The samples were then transported to the laboratory in a cooler at 4 °C. Each of the 1-L water samples was filtered through GF/F glass fiber filter (normal pore size = 0.7 μm; diameter = 47 mm; GE Healthcare Japan Corporation, Tokyo, Japan) and divided into two parts (maximum 500-ml water per 1 GF/F filter). To prevent cross-contamination among the water samples, the filter funnels, and the measuring cups were bleached after filtration. All filtered samples were stored at -20 ℃ in the freezer until the DNA extraction step.Figure 4Sampling sites used in the present study. Blue circles and orange triangles show the locations of the river and lake samples, respectively. Detailed information on each site is listed in Supplementary Table S1. This map has been illustrated using QGIS ver.3.10 (http://www.qgis.org/en/site/) based on the Administrative Zones Data (http://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-N03-v2_3.html) which were obtained from free download service of the National Land Numerical Information (http://nlftp.mlit.go.jp/ksj/index.html, edited by RN). There was no need of obtaining permissions for editing and publishing of map data.Full size imageDNA extraction and library preparationThe total eDNA was extracted from each filtered sample using the DNeasy Blood and Tissue Kit (QIAGEN, Hilden, Germany). Extraction methods were according to Uchii et al.35, with a few modifications. A filtered sample was placed in the upper part of a Salivette tube and 440 μL of a solution containing 400 μL Buffer AL and 40 μL Proteinase K added. The tube with the filtered sample was incubated at 56 °C for 30 min. Afterward, the tube was centrifuged at 5000 × g for 3 min, and the solution at the bottom part of the tube was collected. To increase eDNA yield, 220-μL Tris–EDTA (TE) buffer was added to the filtered sample and the sample re-centrifuged at 5000 × g for 1 min. Subsequently, 400 μL of ethanol was added to the collected solution, and the mixture was transferred to a spin column. Afterward, the total eDNA was eluted in 100-μL buffer AE according to the manufacturer’s instructions. All eDNA samples were stored at -20 °C until the library preparation step.In the present study, we used a universal primer set “MiFish” for eDNA metabarcoding9. The amplicon library was prepared according to the following protocols. In the first PCR, the total reaction volume was 12 μL, containing 6.0μL 2 × KOD buffer, 2.4 μL dNTPs, 0.2 μL KOD FX Neo (TOYOBO, Osaka, Japan), 0.35 μL MiFish-U-F (5ʹ-ACACTCTTTCCCTACACGACGCTCTTCCGATCTNNNNNNGTCGGTAAAACTCGTGCCAGC-3ʹ), MiFish-U-R (5ʹ-GTGACTGGAGTTCAGACGTGTGCTCTTCCGATCTNNNNNNCATAGTGGGGTATCTAATCCCAGTTTG-3ʹ), MiFish-E-F (5ʹ-ACACTCTTTCCCTACACGACGCTCTTCCGATCTNNNNNNRGTTGGTAAATCTCGTGCCAGC-3ʹ) and MiFish-E-R (5ʹ-GTGACTGGAGTTCAGACGTGTGCTCTTCCGATCTNNNNNNGCATAGTGGGGTATCTAATCCTAGTTTG -3ʹ) primers with Illumina sequencing primer region and 6-mer Ns, and 2 μL template DNA. The thermocycling conditions were 94 ℃ for 2 min, 35 cycles of 98 ℃ for 10 s, 65 ℃ for 30 s, 68 ℃ for 30 s, and 68 ℃ for 5 min. The first PCR was repeated four times for each sample, and the replicated samples were pooled as a single first PCR product for use in the subsequent step. The pooled first PCR products were purified using the Solid Phase Reversible Immobilization select Kit (AMPure XP; BECKMAN COULTER Life Sciences, Indianapolis, IN, USA) according to the manufacturer’s instructions. The DNA concentrations of purified first PCR products were measured using a Qubit dsDNA HS assay kit and a Qubit 3.0 fluorometer (Thermo Fisher Scientific, Waltham, MA, USA). All purified first PCR products were diluted to 0.1 ng/μL with H2O, and the diluted samples were used as templates for the second PCR. In the first PCR step, the PCR negative controls (four replicates) were included in each experiment. A total of three PCR negative controls were included in the library (PCR Blank 1–3 samples in Supplementary Table S1, S2, S4, and S5).The second PCR was performed to add HTS adapter sequences with 8-bp dual indices. The total reaction volume was 12 μL, containing 6.0 μL 2 × KAPA HiFi HotStart ReadyMix, 1.4 μL forward and reverse primer (2.5 μM), 1 μL purified first PCR product, and 2.2 μL H2O. The thermocycling conditions were 95 ℃ for 3 min, 12 cycles of 98 ℃ for 20 s, 72 ℃ for 15 s, and 72 ℃ for 5 min.Each Indexed second PCR product was pooled in the equivalent volume, and 25 μL of the pooled libraries were loaded on a 2% E-Gel SizeSelect agarose gels (Thermo Fisher Scientific), and a target library size (ca. 370 bp) was collected. The quality of the amplicon library was checked using an Agilent 2100 Bioanalyzer and Agilent 2100 Expert (Agilent Technologies Inc., Santa Clara, CA, USA), and the DNA concentrations of the amplicon library were measured using Qubit dsDNA HS assay Kit using a Qubit 3.0 fluorometer.High-throughput sequencingAmplicon library was sequenced using iSeq and MiSeq platforms (Illumina, San Diego, CA, USA). To normalize the percentage of pass-filtered read numbers, the sequencing runs using the same libraries were performed using iSeq i1 Reagent and MiSeq Reagent Kit v2 Micro. Both sequencing was performed with 8 million pair-end reads and 2 × 150 bp read lengths. Each library was spiked with approximately 20% PhiX control (PhiX Control Kit v3, Illumina, San Diego, CA, USA) before sequencing runs according to the recommendation of Illumina. The wells of cartridges in the iSeq run were loaded with 20 μL of 50 pM library pool, and sequencing performed at Yamaguchi University, Yamaguchi, Japan. The wells of cartridges for MiSeq runs were loaded with 600 μL of 16 pM library pool, and sequencing performed at Illumina laboratories (Minato-ku, Tokyo, Japan). Subsequently, the sequencing dataset outputs from iSeq and MiSeq were subjected to pre-processing and taxonomic assignments. All sequence data are registered in the DNA Data Bank of Japan (DDBJ) Sequence Read Archive (DRA, Accession number: DRA10593).Pre-processing and taxonomic assignmentsWe used the USEARCH v11.066736 for all data pre-processing activities and taxonomic assignment of the HTS datasets obtained from the iSeq and MiSeq platforms16,37. First, pair-end reads (R1 and R2 reads) generated from iSeq and MiSeq platforms were assembled using the “fastq_mergepairs” command with a minimum overlap of 10 bp. In the process, the low-quality tail reads with a cut-off threshold at a Phred score of 2, and the paired reads with too many mismatches ( > 5 positions) in the aligned regions were discarded38. Secondly, the primer sequences were removed from the merged reads using the “fastx_truncate” command. Afterward, read quality filtering was performed using the “fastq_filter” command with thresholds of max expected error  > 1.0 and  > 50 bp read length. The pre-processed reads were dereplicated using the “fastx_uniques” command, and the chimeric reads and less than 10 reads were removed from all samples as the potential sequence errors. Finally, an error-correction of amplicon reads, which checks and discards the PCR errors and chimeric reads, was performed using the “unoise3” command in the unoise3 algorithm39. Before the taxonomic assignment, the processed reads from the above steps were subjected to sequence similarity search using the “usearch_global” command against reference databases of fish species that had been established previously (MiFish local database v34). The sequence similarity and cut off E-value were 99% and 10–5, respectively. If there was only one species with ≧ 99% similarity, the sequence was assigned to the top-hit species. Conversely, sequences assigned to two or more species in the ≧ 99% similarity were merged as species complex and listed in the synonym group. Generally, the species complexes were assigned to the genus level (e.g., Asian crucian carp Carassius spp.). Species that were unlikely to inhabit Japan were excluded from the candidate list of species complexes. For example, the sequence of one of bitterling Acheilignathus macropterus included other different two species, A. barbatus and A. chankaensis, as the species of the 2nd hit candidate; however, the two species are not currently found in Japan. Therefore, the sequence was assigned to A. macropterus in the present study. Because we used only freshwater fish species, we removed the operational taxonomic units (OTUs) assigned to marine and brackish fishes from each sample. Finally, sequence reads of each fish species were arranged into the matrix, with the rows and columns representing the number of sites and fish species (or genus), respectively.We evaluated sequence quality based on (1) the percentage of clustering passing filter (% PF) and (2) sequencing quality score ≧ % Q30 (Read1 and Read2) between iSeq and MiSeq platforms. The % PF value is an indicator of signal purity for each cluster40. The condition leads to poor template generation, which decreases the % PF value40. In the present study, a  > 80% PF value was set as the threshold of sequence quality in iSeq and MiSeq runs. Sequence quality scores (Q score) measure the probability that a base is called incorrectly. Higher Q scores indicate lower probability of sequencing error, and lower Q scores indicate probability of false-positive variant calls resulting in inaccurate conclusions41. In the present study, the % Q30 values (error rate = 0.001%) were used for the comparison of sequence quality between iSeq and MiSeq. The parameters were collected directly using Illumina BaseSpace Sequence Hub. We also evaluated changes in sequence reads in pre-processing steps between iSeq and MiSeq platforms. Sequence reads were assessed based (1) merge pairs, (2) quality filtering, and (3) denoising. In each step, the change in the number of reads before and after processing was calculated. The calculated numbers of sequence reads are listed in Supplementary Table S2 and S3 in series.Comparing sequence quality and fish fauna between iSeq and MiSeqTo test a relationship of remained sequence reads between iSeq and MiSeq in each pre-processing part, we performed spearman’s rank correlation test in each step. In the present study, however, the sequencing run by iSeq and MiSeq was performed only once each for the same sample. Therefore, we could not assess the variabilities of the sequence read in quality checks and taxonomic assignment in the same samples between iSeq and MiSeq.Before the comparison of fish fauna, rarefaction curves were illustrated for each sample in both iSeq and MiSeq to confirm that the sequencing depth adequately covered the species composition using the “rarecurve” function of the “vegan” package ver. 2.5–6 (https://github.com/vegandevs/vegan) in R ver. 3.6.242. In the present study, the differences in the numbers of sequence reads among samples were confirmed in the two sequencers, but rarefaction curves were saturated in all iSeq and MiSeq samples (Supplementary Fig. S6 and S7). We performed a rarefaction using the “rrarefy” function in “vegan” package to match up the iSeq sequence depths of each sample with that of MiSeq. However, the number of species in each sample on the iSeq have not changed before or after the rarefaction. Therefore, we have used the raw data set before the rarefaction for the subsequent analyses.We compared the species detection capacities of iSeq and MiSeq based on environmental DNA metabarcoding. Using fish faunal data obtained from iSeq and MiSeq, non-metric multidimensional scaling (NMDS) was performed in 1000 separate runs using the “metaMDS” functions in the “vegan” package ver. 2.5–6. For NMDS, the dissimilarity of the fish fauna was calculated based on the incidence-based Jaccard indices. To evaluate the differences in species composition and variance across sites between the two HTS, we performed a permutational multivariate analysis of variance (PERMANOVA) and the permutational analyses of multivariate dispersions (PERMDISP) with 10,000 permutations, respectively. For the PERMANOVA and PERMDISP, we used the “adonis”, and “betadisper” functions in the “vegan” package ver. 2.5-6.Comparison of fish species detectability between eDNA metabarcoding and conventional methodsWe evaluated species detectability between the two HTS by comparing the fish species lists of the two HTS with lists from conventional methods. Five sampling sites were selected from Kyushu and Chugoku districts (R23–27 in Fig. 4). The fish fauna data obtained by conventional methods were based on the results of a previous study43. The conventional surveys were conducted through hand-net sampling and visual observation by snorkeling (see a previous study43 for the detailed methods). The count data of each species were replaced with the incidence-based datasets (presence or absence) for comparing with the eDNA metabarcoding datasets. Fish sequence reads of each sampling site obtained by eDNA metabarcoding were also replaced with the incidence-based data.To test the detectability of species observed by conventional methods, the fish species compositions in five rivers were compared between the eDNA metabarcoding (iSeq and MiSeq) and the conventional methods. To visualize the differences in the species composition between HTSs and conventional methods, heat maps were illustrated for each sampling site. To assess differences in the number of species among methods at each river, the repeated measures analysis of variance (ANOVA) was performed among iSeq, MiSeq, and conventional methods. If a significant difference was found in repeated measures ANOVA, the Tukey–Kramer multiple comparison test was performed to analyze differences among methods.Using fish faunal data obtained from iSeq, MiSeq, and conventional methods, the NMDS was performed in 1000 separate runs with Jaccard indices. The PERMANOVA was performed with 1000 permutations to assess the differences in fish fauna among the methods and sites. Furthermore, to evaluate variance across sites among methods, the PERMDISP was also performed with 1000 permutations. To visualize the number of species in each method and the number of common species between methods, Venn diagrams were illustrated for each river using the “VennDiagram” package ver. 1.6.2 in R44. More

Site descriptionField experiments were conducted at the Changwu Experimental Station (35° 12′ N, 107° 40′ E, altitude 1200 m) located in Shaanxi Province, China. The experimental site was in the typical dryland farming area on the Loess Plateau. Annual precipitation in the area averaged 582 mm between 1957 and 2013, with a mean annual temperature of 9.7 °C over that period. Rainfall and temperature during the two study years are shown in Fig. S1. Soils were generally of the Calcaric Regosol group, according to the FAO/UNESCO soil classification system52, and were composed of 4% sand, 59% silt, and 37% clay53. The 0–20 cm soil properties were the following: pH, 8.4; organic matter content, 11.8 g kg−1; total N content, 0.87 g kg−1; and Olsen-P, 14.4 mg kg−1.Experimental design and field managementTwo-year experiment was arranged in a randomized complete block design with three replicate plots during 2012 and 2013 growing seasons25,54_ENREF_53. The study was conducted using the soybean cultivar (Glycine max L.) cv. Zhonghuang 24 and the maize cultivar (Zea mays L.) cv. Zhengdan 958 grown in cereal–legume agricultural systems. Zhonghuang 24 was bred from Jilin 21 and fendou 31 × Zhongdou 19 (deposition number 2008003); Zhengdan 958 was the offspring of inbred Zheng 58 and Chang 7-2 (deposition number 20000009), which are approved in China. The cropping system treatments were as follows:

1.

Sole-cropped soybean (S).

2.

Sole-cropped maize (M).

3.

Two rows of maize intercropped with two rows of soybean (M2S2).

4.

Two rows of maize intercropped with four rows of soybean (M2S4).

Each plot measured 6 m × 4 m, with row spacing of 50 cm for maize and soybean both in sole crops and intercrops. Individual plants were spaced at 22 cm and 19 cm for maize and soybean, respectively, with one plant per stand for maize and two plants per stand for soybean to attain densities of 90,000 and 210,000 plants ha−1, respectively. In 2012, seeds of maize and soybean were sown on 25 April and harvested on 28 September, and in 2013, seeds were sown on 20 April and harvested on 25 September. Before sowing, basal fertilizer was applied at a rate of 90 kg N ha−1 as urea (46% N) and 150 kg P2O5 ha−1 as superphosphate (12%, P2O5), and then additional fertilizers were uniformly spread in each plot, which were then ploughed into the 0–30 cm soil layer using a rotary tiller. All of the plots received 67.5 kg N ha−1 as urea at the bell and silking stages using a hole-seeding machine. No irrigation was applied, and weeds were removed by hand when sighted. The research on plants complied with relevant institutional, national, and international guidelines and legislation.Above- and below-ground measurementsThe Pn was measured with a LI-6400 portable photosynthesis system (LI-COR Inc., Lincoln, NE, USA) from 9:00 to 11:00 h at 120 days after sowing, which corresponds to the milk stage in maize and full seed stage in soybean7,13. We measured photosynthesis of ear leaves of maize, the first spreading leaves at the top of soybean in both the sole crops and intercrops. The Pn values were calculated as the sum of the mean readings for five leaves in each plot. The LAI values, DIFN were recorded using a Plant Canopy Analyzer (Li-2200, LiCor Inc., Lincoln, NE, USA) without direct sunlight at milk stage of maize. One above-canopy measurement and three below-canopy measurements at the soil surface were taken for four replicates in each plot. SPAD were collected using a hand-held dual wavelength meter (SPAD 502, Chlorophyll meter, Minolta Camera Co., Ltd., Japan) at milk stage of maize. Measurements were taken midway along the ear leaves of maize and the first spreading leaves at the top of soybean from five adjacent plants at the center of row in each plot.The SWS was measured gravimetrically using a soil auger at 10 cm intervals over a depth of 100 cm and at 20 cm intervals over a depth of 200 cm at milk stage of maize for three replicates in each plot. The SWS was calculated for each plot in the 0–200 cm soil profile for the soil moisture using the following formula: SWS = SWC × SD × SBD, where SWC represents soil water content, SD represents soil depth, and SBD represents soil bulk density. Apparent water use during crop growth season was expressed as evapotranspiration (ET), which was determined according to the following formula: ET = ΔSWS + P, where ΔSWS is the change in soil water storage in the top 200 cm and P is the rainfall (mm) between planting and at milk stage in maize. The six adjacent plant samples were collected at milk stage of maize in the middle two rows of each plots (Fig. S2). The sampling included shoots and roots of maize and soybean. At the cotyledonary node, above-ground parts were separated from below-ground parts. Soil core samples (9 cm diameter × 15 cm) at the intra-row of crop were collected to a depth of 100 cm using an auger and separated in 10-cm sections to determine the root growth in sole-cropping and intercropping systems. The samples were exposed to 105 °C for 30 min and then dried to a constant weight at 75 °C. The oven-dried samples were put in small plastic bags after grinding. The study of N and P uptake are the most common among mineral elements55,56. Concentrations of N and P in the plant dry matter were determined after digestion with H2SO4 and H2O2; N concentration was measured according to the Kjeldahl method20, whereas P concentration was measured by the molybdenum-antimony anti-spectrophotometric method16. Crop N and P uptake were calculated by the actual above-ground biomass multiplied by plant tissue N and P concentrations. Grain yield was estimated at harvest from 6 m2 for maize and soybean based on the average of three plot replicates.Data analysisThe LER for assessment of land use advantage. LER is sum of ratio of intercrop to sole crop for maize and soybean yield57:$$ LER = LER_{m} + LER_{s} ,;LER_{m} = frac{{Y_{im} }}{{Y_{sm} }}, ;LER_{s} = frac{{Y_{is} }}{{Y_{ss} }} $$where LERm and LERs are patial LER for maize and soybean, respectively. Yim and Yis are yields of maize and soybean under intercrops, respectively. Ysm and Yss are the yield of maize and soybean under sole crop, respectively.The water equivalent ratio (WER) was calculated to measure water use advantage of intercropping58:$$ WER = WER_{m} + WER_{s} ,;WER_{m} = frac{{Y_{im} /ET_{im} }}{{Y_{sm} /ET_{sm} }},;WER_{s} = frac{{Y_{is} /ET_{is} }}{{Y_{ss} /ET_{ss} }} $$where WERm and WERs are patial WER for maize and soybean, respectively. ETim and ETis are ET of maize and soybean under intercrops, respectively. ETsm and ETss are the ET of maize and soybean under sole crop, respectively.All analyses were conducted in SPSS Statistics 17.0 (SPSS Inc., Chicago, IL, USA). Treatment means showing significant differences among different cropping systems were separated using one-way ANOVA or least significant difference (LSD) at a threshold of 5% to compare the effect of yield, above- and below-ground related parameters (Pn, LAI, SPAD, DIFN, SWS, N and P uptake) in different maize–soybean intercropping. The variation in Pn, LAI, SPAD, DIFN, SWS, N, and P uptake of crop, and the effects of cropping system × year were made using Univariate General Linear Models. Pearson’s correlation test was used to analyze between LER and above-and below-ground biomass of maize and soybean. The effects of above- and below-ground factors on biological yield were quantified, by calculating the contribution value of some key factors to yield. The effects of between above- (LAI, SPAD, DIFN) and below-ground (SWS, N and P uptake) competition on the biological yield and contribution rate were conducted by the linear regression model59:$$ Y = beta_{0} LAI + beta_{1} SPAD + beta_{2} DIFN + beta_{3} SWS + beta_{4} {text{N}} + beta_{5} {text{P}} + beta_{6} X + beta_{7} $$
(1)
where Y represents biological yield, LAI represents leaf area index, SPAD represents chlorophyll, DIFN represents diffuse non interceptance, SWS represents soil water storage, N represents crop nitrogen uptake, P represents crop phosphorus uptake, X represents interaction for LAI, SPAD, DIFN, SWS, N, and P, and β0, β1, β2, β3, β4, β5, β6 and β7 represent the fitted parameters. The standard regression coefficients (Beta) of LAI, SPAD, DIFN, SWS, N, and P were determined on the basis of Eq. (1) to split their influence on the biological yield by the following equations:$$ beta_{0}^{prime } = beta_{0} times left( {LAI^{prime } /Y^{prime } } right) $$
(2)
$$ beta_{1}^{prime } = beta_{1} times left( {SPAD^{prime } /Y^{prime } } right) $$
(3)
$$ beta_{2}^{prime } = beta_{2} times left( {DIFN^{prime } /Y^{prime } } right) $$
(4)
$$ beta_{3}^{prime } = beta_{3} times left( {SWS^{prime } /Y^{prime } } right) $$
(5)
$$ beta_{4}^{prime } = beta_{4} times left( {{text{N}}^{prime } /Y^{prime } } right) $$
(6)
$$ beta_{5}^{prime } = beta_{5} times left( {{text{P}}^{prime } /Y^{prime } } right) $$
(7)
where β0′, β1′, β2′, β3′, β4′, and β5′ represent the standard regression coefficients for LAI, SPAD, DIFN, SWS, N, and P. LAI′, SPAD′, DIFN′, SWS′, N′, and P′ represent the standard deviations for LAI, SPAD, DIFN, SWS, N, and P. Y′ is the standard deviation for the modeled biological yield. More

1.Litchman, E., Ohman, M. D. & Kiørboe, T. Trait-based approaches to zooplankton communities. J. Plankton Res. 35, 473–484. https://doi.org/10.1093/plankt/fbt019 (2013).Article 

Google Scholar 
2.Kiørboe, T. & Hirst, A. G. Shifts in mass scaling of respiration, feeding, and growth rates across life-form transitions in marine pelagic organisms. Am. Nat. 183, E118–E130. https://doi.org/10.1086/675241 (2014).Article 
PubMed 

Google Scholar 
3.Andersen, K. H. et al. Characteristic sizes of life in the oceans, from bacteria to whales. Annu. Rev. Mar. Sci. 8, 1–25. https://doi.org/10.1146/annurev-marine-122414-034144 (2015).Article 

Google Scholar 
4.Bergmann, C. Über die Verhältnisse der Wärmeökonomie der Thiere zu ihrer Grösse. Göttinger Studien 3, 595–708 (1847).
Google Scholar 
5.Woodson, C., Schramski, J. R. & Joye, S. B. A unifying theory for top-heavy ecosystem structure in the ocean. Nat. Commun. 9, 1–8. https://doi.org/10.1038/s41467-017-02450-y (2018).CAS 
Article 

Google Scholar 
6.Brown, J. H., Gillooly, J. F., Allen, A. P., Savage, V. M. & West, G. B. Toward a metabolic theory of ecology. Ecology 85, 1771–1789. https://doi.org/10.1890/03-9000 (2004).Article 

Google Scholar 
7.Gardner, J. L., Peters, A., Kearney, M. R., Joseph, L. & Heinsohn, R. Declining body size: A third universal response to warming?. Trends Ecol. Evol. 26, 285–291. https://doi.org/10.1016/j.tree.2011.03.005 (2011).Article 
PubMed 

Google Scholar 
8.Angilletta, M. J., Steury, T. D. & Sears, M. W. Temperature, growth rate, and body size in ectotherms: Fitting pieces of a life-history puzzle. Integr. Comp. Biol. 44, 498–509. https://doi.org/10.1093/icb/44.6.498 (2004).Article 
PubMed 

Google Scholar 
9.Atkinson, D. Temperature and organism size: A biological law for ectotherms?. Adv. Ecol. Res. 25, 1–58. https://doi.org/10.1016/S0065-2504(08)60212-3 (1994).Article 

Google Scholar 
10.Atkinson, D. & Sibly, R. M. Why are organisms usually bigger in colder environments? Making sense of a life history puzzle. Trends Ecol. Evol. 12, 235–239. https://doi.org/10.1016/S0169-5347(97)01058-6 (1997).CAS 
Article 
PubMed 

Google Scholar 
11.Sunagawa, S. et al. Structure and function of the global ocean microbiome. Science 348, 1261359. https://doi.org/10.1126/science.1261359 (2015).CAS 
Article 
PubMed 

Google Scholar 
12.Audzijonyte, A. et al. Is oxygen limitation in warming waters a valid mechanism to explain decreased body sizes in aquatic ectotherms?. Glob. Ecol. Biogeogr. 28, 64–77. https://doi.org/10.1111/geb.12847 (2018).Article 

Google Scholar 
13.Begon, M., Townsend, C. R. & Harper, J. L. Ecology: From Individuals to Ecosystems 4th edn. (Blackwell Publishing, New York, 2006).
Google Scholar 
14.Hirata, T., Aiken, J., Hardman-Mountford, N., Smyth, T. J. & Barlow, R. G. An absorption model to determine phytoplankton size classes from satellite ocean colour. Remote Sens. Environ. 112, 3153–3159. https://doi.org/10.1016/j.rse.2008.03.011 (2008).ADS 
Article 

Google Scholar 
15.Kostadinov, T., Siegel, D. & Maritorena, S. Global variability of phytoplankton functional types from space: Assessment via the particle size distribution. Biogeosciences 7, 3239–3257. https://doi.org/10.5194/bg-7-3239-2010 (2010).ADS 
Article 

Google Scholar 
16.Brun, P., Payne, M. R. & Kiørboe, T. Trait biogeography of marine copepods: An analysis across scales. Ecol. Lett. 19, 1403–1413. https://doi.org/10.1111/ele.12688 (2016).Article 
PubMed 

Google Scholar 
17.Horne, C. R., Hirst, A. G., Atkinson, D., Neves, A. & Kiørboe, T. A global synthesis of seasonal temperature–size responses in copepods. Glob. Ecol. Biogeogr. 25, 988–999. https://doi.org/10.1111/geb.12460 (2016).Article 

Google Scholar 
18.Garzke, J., Hansen, T., Ismar, S. M. H. & Sommer, U. Combined effects of ocean warming and acidification on copepod abundance, body size and fatty acid content. PLoS ONE 11, e0155952. https://doi.org/10.1371/journal.pone.0155952 (2006).CAS 
Article 

Google Scholar 
19.Stelzer, C. P. Phenotypic plasticity of body size at different temperatures in a planktonic rotifer: Mechanisms and adaptive significance. Funct. Ecol. 16, 835–841. https://doi.org/10.1046/j.1365-2435.2002.00693.x (2002).Article 

Google Scholar 
20.Riemer, K., Anderson-Teixeira, K. J., Smith, F. A., Harris, D. J. & Ernest, S. K. M. Body size shifts influence effects of increasing temperatures on ectotherm metabolism. Global Ecol. Biogeogr. 27, 958–967. https://doi.org/10.1111/geb.12757 (2018).Article 

Google Scholar 
21.Gorsky, G. et al. Digital zooplankton image analysis using the ZooScan integrated system. J. Plankton Res. 32, 285–303. https://doi.org/10.1093/plankt/fbp124 (2010).Article 

Google Scholar 
22.Ibarbalz, F. M. et al. Global trends in marine plankton diversity across kingdoms of life. Cell 179, 1084–1097. https://doi.org/10.1016/j.cell.2019.10.008 (2019).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
23.Hoefnagel, K. N. & Verberk, W. C. Is the temperature-size rule mediated by oxygen in aquatic ectotherms?. J. Therm. Biol. 54, 56–65. https://doi.org/10.1016/j.jtherbio.2014.12.003 (2015).Article 
PubMed 

Google Scholar 
24.Wojewodzic, M. W., Kyle, M., Elser, J. J., Hessen, D. O. & Andersen, T. Joint effect of phosphorus limitation and temperature on alkaline phosphatase activity and somatic growth in Daphnia magna. Oecologia 165, 837–846. https://doi.org/10.1007/s00442-010-1863-2 (2011).ADS 
Article 
PubMed 

Google Scholar 
25.Gillooly, J. F., Brown, J. H., West, G. B., Savage, V. M. & Charnov, E. L. Effects of size and temperature on metabolic rate. Science 293, 2248–2251. https://doi.org/10.1126/science.1061967 (2001).ADS 
CAS 
Article 
PubMed 

Google Scholar 
26.Czarnoleski, M., Ejsmont-Karabin, J., Angilletta, M. K. & Kozlowski, J. Colder rotifers grow larger but only in oxygenated waters. Ecosphere 6, 1–5. https://doi.org/10.1890/ES15-00024.1 (2015).Article 

Google Scholar 
27.Kiørboe, T. How zooplankton feed: Mechanisms, traits and trade-offs. Biol. Rev. 86, 311–339. https://doi.org/10.1111/j.1469-185X.2010.00148.x (2011).Article 
PubMed 

Google Scholar 
28.Benedetti, F., Gasparini, S. & Ayata, S.-D. Identifying copepod functional groups from species functional traits. J. Plankton Res. 38, 159–166. https://doi.org/10.1093/plankt/fbv096 (2016).Article 
PubMed 

Google Scholar 
29.Brun, P., Payne, M. R. & Kiørboe, T. A trait database for marine copepods. Earth Syst. Sci. Data 9, 99–113. https://doi.org/10.5194/essd-9-99-2017 (2017).ADS 
Article 

Google Scholar 
30.Anderson, T. R. Plankton functional type modelling: Running before we can walk?. J. Plankton Res. 27, 1073–1081. https://doi.org/10.1093/plankt/fbi076 (2005).ADS 
Article 

Google Scholar 
31.Biard, T. et al. In situ observations unveil an unexpectedly large biomass of Radiolaria and Phaeodaria (Rhizaria) in the oceans. Nature 532, 504–507. https://doi.org/10.1038/nature17652 (2016).ADS 
CAS 
Article 
PubMed 

Google Scholar 
32.Takagi, H. et al. Characterizing photosymbiosis in modern planktonic foraminifera. Biogeosciences 16, 3377–3396. https://doi.org/10.5194/bg-16-3377-2019 (2019).ADS 
CAS 
Article 

Google Scholar 
33.Rink, S., Kühl, M., Bijma, J. & Spero, H. J. Microsensor studies of photosynthesis and respiration in the symbiotic foraminifer Orbulina universa. Mar. Biol. 131, 583–595. https://doi.org/10.1007/s002270050350 (1998).Article 

Google Scholar 
34.Lombard, F., Erez, J., Michel, E. & Labeyrie, L. Temperature effect on respiration and photosynthesis of the symbiont-bearing planktonic foraminifera Globigerinoides ruber, Orbulina universa, and Globigerinella siphonifera. Limnol. Oceanogr. 54, 210–218. https://doi.org/10.4319/lo.2009.54.1.0210 (2009).ADS 
CAS 
Article 

Google Scholar 
35.Lesser, M. P. Coral Bleaching: Causes and Mechanisms. In Coral Reefs: An Ecosystem in Transition (eds Dubinsky, Z. & Stambler, N.) (Springer, 2011).
Google Scholar 
36.Villar, E. et al. Symbiont chloroplasts remain active during bleaching-like response induced by thermal stress in Collozoum pelagicum (Collodaria, Retaria). Front. Mar. Sci. 5, 387. https://doi.org/10.3389/fmars.2018.00387 (2018).Article 

Google Scholar 
37.Hemleben, C., Spindler, M. & Anderson, O. R. Modern Planktonic Foraminifera (Springer-Verlag, 1989).Book 

Google Scholar 
38.Suzuki, N. & Not, F. Biology and ecology of radiolaria. In Marine Protists: Diversity and Dynamics (eds Ohtsuka, S. et al.) (Springer, 2015).
Google Scholar 
39.de Puelles, F. et al. Zooplankton abundance and diversity in the tropical and subtropical ocean. Diversity 11, 203. https://doi.org/10.3390/d11110203 (2019).CAS 
Article 

Google Scholar 
40.Beaugrand, G., Edwards, M. & Legendre, L. Marine biodiversity, ecosystem functioning, and carbon cycles. PNAS 107, 10120–10124. https://doi.org/10.1073/pnas.0913855107 (2010).ADS 
Article 
PubMed 
PubMed Central 

Google Scholar 
41.Brun, P. et al. Climate change has altered zooplankton-fuelled carbon export in the North Atlantic. Nat. Ecol. Evol. 3, 416–423. https://doi.org/10.1038/s41559-018-0780-3 (2019).Article 
PubMed 

Google Scholar 
42.Buitenhuis, E. T., Le Quéré, C., Bednaršek, N. & Schiebel, R. Large contribution of pteropods to shallow CaCO3 export. Glob. Biogeochem. Cyc. 33, 458–468. https://doi.org/10.1029/2018GB006110 (2019).ADS 
CAS 
Article 

Google Scholar 
43.Follows, M. J., Dutkiewicz, J., Grant, S. & Chisholm, S. W. Emergent biogeography of microbial communities in a model ocean. Science 315, 1843–1846. https://doi.org/10.1126/science.1138544 (2007).ADS 
CAS 
Article 
PubMed 

Google Scholar 
44.Ward, B. A., Dutkiewicz, S., Jahn, O. & Follows, J. F. A size-structured food-web model for the global ocean. Limnol. Oceanogr. 57, 1877–1891. https://doi.org/10.4319/lo.2012.57.6.1877 (2012).ADS 
Article 

Google Scholar 
45.Sailley, S. F. et al. Comparing food web structures and dynamics across a suite of global marine ecosystem models. Ecol. Model. 261, 43–57. https://doi.org/10.1016/j.ecolmodel.2013.04.006 (2013).Article 

Google Scholar 
46.Le Quéré, C. et al. Role of zooplankton dynamics for Southern Ocean phytoplankton biomass and global biogeochemical cycles. Biogeosciences 13, 4111–4133. https://doi.org/10.5194/bg-13-4111-2016 (2016).ADS 
CAS 
Article 

Google Scholar 
47.Kwiatkowski, L. et al. Emergent constraints on projections of declining primary production in the tropical oceans. Nat. Clim. Change 7, 355–358. https://doi.org/10.1038/nclimate3265 (2017).ADS 
CAS 
Article 

Google Scholar 
48.Sunagawa, S. et al. Tara Oceans: Towards global ocean ecosystems biology. Nat. Rev. Microbiol. 18, 428–445. https://doi.org/10.1038/s41579-020-0364-5 (2020).CAS 
Article 
PubMed 

Google Scholar 
49.Pesant, S. et al. Open science resources for the discovery and analysis of Tara Oceans data. Sci. Data 2, 150023. https://doi.org/10.1038/sdata.2015.23 (2015).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
50.Picheral, M. et al. Vertical profiles of environmental parameters measured on discrete water samples collected with Niskin bottles at station TARA_147 during the Tara Oceans expedition 2009–2013. PANGAEA https://doi.org/10.1594/PANGAEA.839235 (2014).51.Guidi, L. et al. Environmental context of all samples from the Tara Oceans Expedition (2009–2013), about sensor data in the targeted environmental feature. PANGAEA https://doi.org/10.1594/PANGAEA.875576 (2017).52.Guidi, L. et al. Environmental context of all samples from the Tara Oceans Expedition (2009–2013), about pigment concentrations (HPLC) in the targeted environmental feature. PANGAEA https://doi.org/10.1594/PANGAEA.875569 (2017).53.Guidi, L. et al. Environmental context of all samples from the Tara Oceans Expedition (2009–2013), about nutrients in the targeted environmental feature. PANGAEA https://doi.org/10.1594/PANGAEA.875575 (2017).54.Speich, S. et al. Environmental context of all samples from the Tara Oceans Expedition (2009–2013), about the water column features at the sampling location. PANGAEA https://doi.org/10.1594/PANGAEA.875579 (2017).55.de Boyer-Montegut, C., Madec, G., Fischer, A. S., Lazar, A. & Iudicone, D. Mixed layer depth over the global ocean: An examination of profile data and a profile-based climatology. J. Geophys. Res. 109, C12003. https://doi.org/10.1029/2004JC002378 (2004).ADS 
Article 

Google Scholar 
56.Aminot, A., Kérouel, R. & Coverly, S. C. Nutrients in seawater using segmented flow analysis. In Practical Guidelines for the Analysis of Seawater (ed. Wurl, O.) (CRC Press, 2009).
Google Scholar 
57.Uitz, J., Claustre, H., Morel, A. & Hooker, S. B. Vertical distribution of phytoplankton communities in Open Ocean: An assessment based on surface chlorophyll. J. Geophys. Res. 111, C08005. https://doi.org/10.1029/2005JC003207 (2006).ADS 
Article 

Google Scholar 
58.Pante, E. & Simon-Bouhet, B. marmap: A Package for importing, plotting and analyzing bathymetric and topographic data in R. PLoS ONE 8(9), e73051. https://doi.org/10.1371/journal.pone.0073051 (2013).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
59.Picheral, M., Colin, S. & Irisson J.-O. EcoTaxa, A Tool for the Taxonomic Classification of Images. http://ecotaxa.obs-vlfr.fr (2017).60.Wood, S. N. Generalized Additive Models: An Introduction with R 2nd edn. (Chapman and Hall/CRC, 2017).Book 

Google Scholar 
61.Dormann, C. F. et al. Collinearity: A review of methods to deal with it and a simulation study evaluating their performance. Ecography 36, 27–46. https://doi.org/10.1111/j.1600-0587.2012.07348.x (2013).Article 

Google Scholar 
62.Giorgino, T. Computing and visualizing dynamic time warping alignments in R: The dtw package. J. Stat. Softw. 31, 1–24. https://doi.org/10.18637/jss.v031.i07 (2009).Article 

Google Scholar 
63.Park, H.-S. & Jun, C.-H. A simple and fast algorithm for K-medoids clustering. Expert Syst. Appl. 36, 3336–3341. https://doi.org/10.1016/j.eswa.2008.01.039 (2009).Article 

Google Scholar 
64.R Core Team. R: A Language and Environment for Statistical Computing. (R Foundation for Statistical Computing, 2018). https://www.R-project.org/.65.Wickham, H. et al. Welcome to the Tidyverse. J. Open Source Softw. 4, 1686. https://doi.org/10.21105/joss.01686 (2019).ADS 
Article 

Google Scholar 
66.Heiberger, R. M. HH: Statistical Analysis and Data Display: Heiberger and Holland. R package version 3.1–40, https://CRAN.R-project.org/package=HH (2020).67.Lê, S., Josse, J. & Husson, F. FactoMineR: An R package for multivariate analysis. J. Stat. Softw. 25, 1–18. https://doi.org/10.18637/jss.v025.i01 (2008).Article 

Google Scholar 
68.Sarda-Espinosa, A. dtwclust: Time Series Clustering Along with Optimizations for the Dynamic Time Warping Distance. R package version 5.5.6, https://CRAN.R-project.org/package=dtwclust (2019). More

1.Luyssaert, S. et al. Old-growth forests as global carbon sinks. Nature 455, 213–215 (2008).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
2.Pan, Y. et al. A large and persistent carbon sink in the world’s forests. Science 333, 988–993 (2011).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
3.Rinne-Garmston, K. T. et al. Carbon flux from decomposing wood and its dependency on temperature, wood N2 fixation rate, moisture and fungal composition in a Norway spruce forest. Glob. Chang. Biol. 25, 1852–1867 (2019).ADS 
PubMed 
PubMed Central 
Article 

Google Scholar 
4.Šamonil, P. et al. Convergence, divergence or chaos? Consequences of tree trunk decay for pedogenesis and the soil microbiome in a temperate natural forest. Geoderma 376, 114499 (2020).ADS 
Article 

Google Scholar 
5.Tláskal, V. et al. Complementary roles of wood-inhabiting fungi and bacteria facilitate deadwood decomposition. mSystems 6, e01078–20 (2021).PubMed 
PubMed Central 
Article 

Google Scholar 
6.Odriozola, I. et al. Fungal communities are important determinants of bacterial community composition in deadwood. mSystems 6, e01017–20 (2021).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
7.Valášková, V., de Boer, W., Gunnewiek, P. J. A. K., Pospíšek, M. & Baldrian, P. Phylogenetic composition and properties of bacteria coexisting with the fungus Hypholoma fasciculare in decaying wood. ISME J. 3, 1218–1221 (2009).PubMed 
Article 
CAS 
PubMed Central 

Google Scholar 
8.Brunner, A. & Kimmins, J. P. Nitrogen fixation in coarse woody debris of Thuja plicata and Tsuga heterophylla forests on northern Vancouver Island. Can. J. For. Res. 33, 1670–1682 (2003).CAS 
Article 

Google Scholar 
9.Rinne, K. T. et al. Accumulation rates and sources of external nitrogen in decaying wood in a Norway spruce dominated forest. Funct. Ecol. 31, 530–541 (2016).Article 

Google Scholar 
10.Põlme, S. et al. FungalTraits: a user-friendly traits database of fungi and fungus-like stramenopiles. Fungal Divers. 105, 1–16 (2020).Article 

Google Scholar 
11.Tláskal, V. & Baldrian, P. Deadwood-inhabiting bacteria show adaptations to changing carbon and nitrogen availability during decomposition. Front. Microbiol. 12, 685303 (2021).PubMed 
PubMed Central 
Article 

Google Scholar 
12.Lemos, L. N., Mendes, L. W., Baldrian, P. & Pylro, V. S. Genome-resolved metagenomics is essential for unlocking the microbial black box of the soil. Trends Microbiol. 29, 279–282 (2021).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
13.Větrovský, T. et al. GlobalFungi, a global database of fungal occurrences from high-throughput-sequencing metabarcoding studies. Sci. Data 7, 228 (2020).PubMed 
PubMed Central 
Article 

Google Scholar 
14.Thompson, L. R. et al. A communal catalogue reveals Earth’s multiscale microbial diversity. Nature 551, 457–463 (2017).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
15.Anderson-Teixeira, K. J., Davies, S. J., Bennett, A. C., Muller-landau, H. C. & Wright, S. J. CTFS-ForestGEO: a worldwide network monitoring forests in an era of global change. Glob. Chang. Biol. 21, 528–549 (2015).Article 

Google Scholar 
16.Baldrian, P. et al. Fungi associated with decomposing deadwood in a natural beech-dominated forest. Fungal Ecol. 23, 109–122 (2016).Article 

Google Scholar 
17.Smyth, C. E. et al. Patterns of carbon, nitrogen and phosphorus dynamics in decomposing wood blocks in Canadian forests. Plant Soil 9, 46–62 (2016).
Google Scholar 
18.Král, K. et al. Local variability of stand structural features in beech dominated natural forests of Central Europe: Implications for sampling. For. Ecol. Manage. 260, 2196–2203 (2010).Article 

Google Scholar 
19.Caporaso, J. G. et al. Ultra-high-throughput microbial community analysis on the Illumina HiSeq and MiSeq platforms. ISME J. 6, 1621–1624 (2012).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
20.Lanzén, A. et al. CREST – Classification resources for environmental sequence tags. PLoS One 7, e49334 (2012).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
21.Quast, C. et al. The SILVA ribosomal RNA gene database project: improved data processing and web-based tools. Nucleic Acids Res. 41, D590–D596 (2013).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
22.Žifčáková, L., Větrovský, T., Howe, A. & Baldrian, P. Microbial activity in forest soil reflects the changes in ecosystem properties between summer and winter. Environ. Microbiol. 18, 288–301 (2016).PubMed 
Article 
CAS 
PubMed Central 

Google Scholar 
23.Bolger, A. M., Lohse, M. & Usadel, B. Trimmomatic: A flexible trimmer for Illumina sequence data. Bioinformatics 30, 2114–2120 (2014).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
24.Li, D., Liu, C. M., Luo, R., Sadakane, K. & Lam, T. W. MEGAHIT: An ultra-fast single-node solution for large and complex metagenomics assembly via succinct de Bruijn graph. Bioinformatics 31, 1674–1676 (2015).CAS 
Article 

Google Scholar 
25.Kang, D. D., Froula, J., Egan, R. & Wang, Z. MetaBAT, an efficient tool for accurately reconstructing single genomes from complex microbial communities. PeerJ 3, e1165 (2015).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
26.Parks, D. H., Imelfort, M., Skennerton, C. T., Hugenholtz, P. & Tyson, G. W. CheckM: assessing the quality of microbial genomes recovered from isolates, single cells, and metagenomes. Genome Res. 25, 1043–1055 (2015).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
27.Parks, D. H. et al. Recovery of nearly 8,000 metagenome-assembled genomes substantially expands the tree of life. Nat. Microbiol. 2, 1533–1542 (2017).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
28.Parks, D. H. et al. A standardized bacterial taxonomy based on genome phylogeny substantially revises the tree of life. Nat. Biotechnol. 36, 996–1004 (2018).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
29.Lee, M. D. GToTree: A user-friendly workflow for phylogenomics. Bioinformatics 35, 4162–4164 (2019).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
30.Hyatt, D. et al. Prodigal: prokaryotic gene recognition and translation initiation site identification. BMC Bioinformatics 11, 119 (2010).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
31.Eddy, S. R. Accelerated profile HMM searches. PLoS Comput. Biol. 7, e1002195 (2011).ADS 
MathSciNet 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
32.Edgar, R. C. MUSCLE: Multiple sequence alignment with high accuracy and high throughput. Nucleic Acids Res. 32, 1792–1797 (2004).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
33.Capella-Gutiérrez, S., Silla-Martínez, J. M. & Gabaldón, T. trimAl: A tool for automated alignment trimming in large-scale phylogenetic analyses. Bioinformatics 25, 1972–1973 (2009).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
34.Price, M. N., Dehal, P. S. & Arkin, A. P. FastTree 2 – approximately maximum-likelihood trees for large alignments. PLoS One 5, e9490 (2010).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
35.Ihrmark, K. et al. New primers to amplify the fungal ITS2 region – evaluation by 454-sequencing of artificial and natural communities. FEMS Microbiol. Ecol. 82, 666–677 (2012).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
36.Větrovský, T., Baldrian, P. & Morais, D. SEED 2: A user-friendly platform for amplicon high-throughput sequencing data analyses. Bioinformatics 34, 2292–2294 (2018).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
37.Aronesty, E. Comparison of sequencing utility programs. Open Bioinforma. J. 7, 1–8 (2013).MathSciNet 
Article 

Google Scholar 
38.Nilsson, R. H. et al. An open source software package for automated extraction of ITS1 and ITS2 from fungal ITS sequences for use in high-throughput community assays and molecular ecology. Fungal Ecol. 3, 284–287 (2010).Article 

Google Scholar 
39.Edgar, R. C. Search and clustering orders of magnitude faster than BLAST. Bioinformatics 26, 2460–2461 (2010).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
40.Edgar, R. C. UPARSE: highly accurate OTU sequences from microbial amplicon reads. Nat. Methods 10, 996–998 (2013).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
41.Nilsson, R. H. et al. The UNITE database for molecular identification of fungi: handling dark taxa and parallel taxonomic classification. Nucleic Acids Res. 47, D259–D264 (2018).PubMed Central 
Article 
CAS 

Google Scholar 
42.Wright, E. S. Using DECIPHER v2.0 to analyze big biological sequence data in R. R J. 8, 352–359 (2016).Article 

Google Scholar 
43.Murali, A., Bhargava, A. & Wright, E. S. IDTAXA: A novel approach for accurate taxonomic classification of microbiome sequences. Microbiome 6, 140 (2018).PubMed 
PubMed Central 
Article 

Google Scholar 
44.NCBI BioProject https://identifiers.org/ncbi/bioproject:PRJNA603240 (2020).45. NCBI Sequence Read Archive, https://identifiers.org/ncbi/bioproject:PRJNA672674 (2020).46.Sutela, S., Poimala, A. & Vainio, E. J. Viruses of fungi and oomycetes in the soil environment. FEMS Microbiol. Ecol. 95, fiz119 (2019).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
47.Woodcroft, B. J. et al. Genome-centric view of carbon processing in thawing permafrost. Nature 560, 49–54 (2018).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
48.Mackelprang, R. et al. Microbial community structure and functional potential in cultivated and native tallgrass prairie soils of the Midwestern United States. Front. Microbiol. 9, 1775 (2018).PubMed 
PubMed Central 
Article 

Google Scholar 
49.Hervé, V. et al. Phylogenomic analysis of 589 metagenome-assembled genomes encompassing all major prokaryotic lineages from the gut of higher termites. PeerJ 8, e8614 (2020).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
50.Clissmann, F. et al. First insight into dead wood protistan diversity: a molecular sampling of bright-spored Myxomycetes (Amoebozoa, slime-moulds) in decaying beech logs. FEMS Microbiol. Ecol. 91, fiv050 (2015).PubMed 
Article 
CAS 
PubMed Central 

Google Scholar 
51.Urich, T. et al. Simultaneous assessment of soil microbial community structure and function through analysis of the meta-transcriptome. PLoS One 3, e2527 (2008).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
52.Geisen, S. et al. Metatranscriptomic census of active protists in soils. ISME J. 9, 2178–2190 (2015).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
53.Tláskal, V., Zrůstová, P., Vrška, T. & Baldrian, P. Bacteria associated with decomposing dead wood in a natural temperate forest. FEMS Microbiol. Ecol. 93, fix157 (2017).Article 
CAS 

Google Scholar 
54.Moll, J. et al. Bacteria inhabiting deadwood of 13 tree species reveal great heterogeneous distribution between sapwood and heartwood. Environ. Microbiol. 20, 3744–3756 (2018).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
55.Christofides, S. R., Hiscox, J., Savoury, M., Boddy, L. & Weightman, A. J. Fungal control of early-stage bacterial community development in decomposing wood. Fungal Ecol. 42, 100868 (2019).Article 

Google Scholar 
56.Nayfach, S. et al. A genomic catalog of Earth’s microbiomes. Nat. Biotechnol. 39, 499–509 (2021).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
57.Seibold, S. et al. Experimental studies of dead-wood biodiversity — A review identifying global gaps in knowledge. Biol. Conserv. 191, 139–149 (2015).Article 

Google Scholar  More

1.Bartelink, E. J., Mackinnon, A. T., Prince-Buitenhuys, J. R., Tipple, B. J. & Chesson, L. A. Stable isotope forensics as an investigative tool in missing persons investigations BT. In Handbook of Missing Persons (eds Morewitz, S. J. & Sturdy Colls, C.) 443–462 (Springer, 2016).Chapter 

Google Scholar 
2.Ehleringer, J. R. et al. Hydrogen and oxygen isotope ratios in human hair are related to geography. Proc. Natl. Acad. Sci. USA 105, 2788–2793 (2008).CAS 
PubMed 
PubMed Central 
Article 
ADS 

Google Scholar 
3.Laffoon, J. E., Rojas, R. V. & Hofman, C. L. Oxygen and carbon isotope analysis of human dental enamel from the caribbean: Implications for investigating individual origins. Archaeometry 55, 742–765 (2013).CAS 
Article 

Google Scholar 
4.Glick, P. L. Patterns of enamel maturation. J. Dent. Res. 58, 883–895 (1979).CAS 
PubMed 
Article 

Google Scholar 
5.Lacruz, R. S., Habelitz, S., Wright, J. T. & Paine, M. L. Dental enamel formation and implications for oral health and disease. Physiol. Rev. 97, 939–993 (2017).PubMed 
PubMed Central 
Article 

Google Scholar 
6.Longinelli, A. Oxygen isotopes in mammal bone phosphate: A new tool for paleohydrological and paleoclimatological research?. Geochim. Cosmochim. Acta 48, 385–390 (1984).CAS 
Article 
ADS 

Google Scholar 
7.Luz, B., Kolodny, Y. & Horowitz, M. Fractionation of oxygen isotopes between mammalian bone-phosphate and environmental drinking water. Geochim. Cosmochim. Acta 48, 1689–1693 (1984).CAS 
Article 
ADS 

Google Scholar 
8.Levinson, A. A., Luz, B. & Kolodny, Y. Variations in oxygen isotopic compositions of human teeth and urinary stones. Appl. Geochemistry 2, 367–371 (1987).CAS 
Article 

Google Scholar 
9.Dansgaard, W. Stable isotopes in precipitation. Tellus A 16, 436–468 (1964).Article 
ADS 

Google Scholar 
10.Rozanski, K., Araguás-Araguás, L. & Gonfiantini, R. Isotopic patterns in modern global precipitation in Climate Change in Continental Isotopic Records 1–36 (American Geophysical Union, 1993). 11.Terzer, S., Wassenaar, L. I., Araguás-Araguás, L. J. & Aggarwal, P. K. Global isoscapes for δ18O and δ2H in precipitation: Improved prediction using regionalized climatic regression models. Hydrol. Earth Syst. Sci. Discuss. https://doi.org/10.5194/hessd-10-7351-2013 (2013).12.Bowen, G. J. Interpolating the isotopic composition of modern meteoric precipitation. Water Resour. Res. 39, 1–13 (2003).Article 
CAS 

Google Scholar 
13.West, A. G., February, E. C. & Bowen, G. J. Spatial analysis of hydrogen and oxygen stable isotopes (“isoscapes”) in ground water and tap water across South Africa. J. Geochem. Explor. 145, 213–222 (2014).CAS 
Article 

Google Scholar 
14.Ehleringer, J. R. et al. A Framework for the incorporation of isotopes and isoscapes in geospatial forensic investigations in Isoscapes SE-17 (eds. West, J. B., Bowen, G. J., Dawson, T. E. & Tu, K. P.) 357–387 (Springer, 2010).15.Laffoon, J. E. et al. Investigating human geographic origins using dual-isotope (87Sr/86Sr, δ18O) assignment approaches. PLoS ONE 12, e0172562 (2017).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
16.Daux, V. et al. Oxygen isotope fractionation between human phosphate and water revisited. J. Hum. Evol. 55, 1138–1147 (2008).PubMed 
Article 

Google Scholar 
17.Dotsika, E. Correlation between δ18Ow and δ18Οen for estimating human mobility and paleomobility patterns. Sci. Rep. 10, 15439 (2020).CAS 
PubMed 
PubMed Central 
Article 
ADS 

Google Scholar 
18.Chenery, C., Müldner, G., Evans, J., Eckardt, H. & Lewis, M. Strontium and stable isotope evidence for diet and mobility in Roman Gloucester, UK. J. Archaeol. Sci. 37, 150–163 (2010).Article 

Google Scholar 
19.Herrmann, N., Li, Z.-H., Warner, M., Weinand, D. & Soto, M. Isotopic and elemental analysis of the William Bass donated skeletal collection and other modern donated collections. Dep. Justice Rep. 1, 248669 (2015).
Google Scholar 
20.Keller, A. T., Regan, L. A., Lundstrom, C. C. & Bower, N. W. Evaluation of the efficacy of spatiotemporal Pb isoscapes for provenancing of human remains. Forensic Sci. Int. 261, 83–92 (2016).CAS 
PubMed 
Article 

Google Scholar 
21.Pollard, A. M., Pellegrini, M. & Lee-Thorp, J. A. Technical note: some observations on the conversion of dental enamel δ18O(p)values to δ18O(w)to determine human mobility. Am. J. Phys. Anthropol. 145, 499–504 (2011).CAS 
PubMed 
Article 

Google Scholar 
22.Posey, R. Development and Validation of a Spatial Prediction Model for Forensic Geographical Provenancing of Human Remains (University of East Anglia, 2011).23.Warner, M. M., Plemons, A. M., Herrmann, N. P. & Regan, L. A. Refining stable oxygen and hydrogen isoscapes for the identification of human remains in Mississippi. J. Forensic Sci. https://doi.org/10.1111/1556-4029.13575 (2018).Article 
PubMed 

Google Scholar 
24.Chenery, C. A., Pashley, V., Lamb, A. L., Sloane, H. J. & Evans, J. A. The oxygen isotope relationship between the phosphate and structural carbonate fractions of human bioapatite. Rapid Commun. Mass Spectrom. 26, 309–319 (2012).CAS 
PubMed 
Article 
ADS 

Google Scholar 
25.Bell, L. S., Lee-Thorp, J. A. & Elkerton, A. Sailing against the wind. J. Archaeol. Sci. 37, 683–686 (2010).Article 

Google Scholar 
26.Pellegrini, M., Pouncett, J., Jay, M., Pearson, M. P. & Richards, M. P. Tooth enamel oxygen “isoscapes’’ show a high degree of human mobility in prehistoric Britain. Sci. Rep. 6, 1–10 (2016).Article 
CAS 

Google Scholar 
27.Podlesak, D. W., Bowen, G. J., O’Grady, S., Cerling, T. E. & Ehleringer, J. R. δ2H and δ18O of human body water: A GIS model to distinguish residents from non-residents in the contiguous USA. Isotopes Environ. Health Stud. https://doi.org/10.1080/10256016.2012.644283 (2012).Article 
PubMed 

Google Scholar 
28.Bowen, G. J. & Wilkinson, B. Spatial distribution of δ18O in meteoric precipitation. Geology 30, 315–318 (2002).Article 
ADS 

Google Scholar 
29.Lee-Thorp, J. Two decades of progress towards understanding fossilization processes and isotopic signals in calcified tissue minerals. Archaeometry 44, 435–446 (2002).CAS 
Article 

Google Scholar 
30.Bryant, D., Koch, P. L., Froelich, P. N., Showers, W. J. & Genna, B. J. Oxygen isotope partitioning between phosphate and carbonate in mammalian apatite. Geochim. Cosmochim. Acta 60, 5145–5148 (1996).CAS 
Article 
ADS 

Google Scholar 
31.Iacumin, P., Bocherens, H., Mariotti, A. & Longinelli, A. Oxygen isotope analyses of co-existing carbonate and phosphate in biogenic apatite: A way to monitor diagenetic alteration of bone phosphate?. Earth Planet. Sci. Lett. 142, 1–6 (1996).CAS 
Article 
ADS 

Google Scholar 
32.Wright, L. E., Valdés, J. A., Burton, J. H., Douglas Price, T. & Schwarcz, H. P. The children of Kaminaljuyu: Isotopic insight into diet and long distance interaction in Mesoamerica. J. Anthropol. Archaeol. 29, 155–178 (2010).Article 

Google Scholar 
33.Chesson, L. A., Kenyhercz, M. W., Regan, L. A. & Berg, G. E. Addressing data comparability in the creation of combined data sets of bioapatite carbon and oxygen isotopic compositions. Archaeometry 61, 1193–1206 (2019).CAS 
Article 

Google Scholar 
34.International Atomic Energy Agency/World Meteorological Organization. Global Network of Isotopes in Precipitation http://www.iaea.org/water (2020).35.Bowen, G. J. The Online Isotopes in Precipitation Calculator, Version 3.1  https://wateriso.utah.edu/waterisotopes/pages/data_access/oipc.html (2018).36.Tipple, B. J. et al. Stable hydrogen and oxygen isotopes of tap water reveal structure of the San Francisco Bay Area’s water system and adjustments during a major drought. Water Res. 119, 212–224 (2017).CAS 
PubMed 
Article 

Google Scholar 
37.Wang, S. et al. Water source signatures in the spatial and seasonal isotope variation of chinese tap waters. Water Resour. Res. 54, 9131–9143 (2018).Article 
ADS 

Google Scholar 
38.Ghassemi, F. & White, I. Inter-Basin Water Transfer: Case Studies from Australia, United States, Canada, China and India (Cambridge University Press, 2007).Book 

Google Scholar 
39.World Health Organization. Progress on Sanitation and Drinking Water: 2010 Update http://apps.who.int/iris/bitstream/10665/81245/1/9789241505390_eng.pdf (2013).40.Good, S. P. et al. Patterns of local and nonlocal water resource use across the western U.S. determined via stable isotope intercomparisons. Water Resour. Res. https://doi.org/10.1002/2012WR013085 (2014).Article 

Google Scholar 
41.Bowen, G. J., Ehleringer, J. R., Chesson, L. A., Stange, E. & Cerling, T. E. Stable isotope ratios of tap water in the contiguous United States. Water Resour. Res. 43, 1–12 (2007).Article 
CAS 

Google Scholar 
42.Chesson, L. A. et al. Strontium isotopes in tap water from the coterminous USA. Ecosphere 3, 67 (2012).Article 

Google Scholar 
43.Ammer, S. T. M., Bartelink, E. J., Vollner, J. M., Anderson, B. E. & Cunha, E. M. Spatial distributions of oxygen stable isotope ratios in tap water from mexico for region of origin predictions of unidentified border crossers. J. Forensic Sci. 65, 1049–1055 (2020).CAS 
PubMed 
Article 

Google Scholar 
44.Daux, V. et al. Oxygen and hydrogen isotopic composition of tap waters in France. Geol. Soc. Lond. Spec. Publ. 507, 207 (2021).
Google Scholar 
45.Gautam, M. K., Song, B. Y., Shin, W. J., Bong, Y. S. & Lee, K. S. Spatial variations in oxygen and hydrogen isotopes in waters and human hair across South Korea. Sci. Total Environ. 726, 138365 (2020).CAS 
PubMed 
Article 
ADS 

Google Scholar 
46.Zhao, S. et al. Divergence of stable isotopes in tap water across China. Sci. Rep. 7, 43653 (2017).PubMed 
PubMed Central 
Article 
ADS 

Google Scholar 
47.Ueda, M. & Bell, L. S. A city-wide investigation of the isotopic distribution and source of tap waters for forensic human geolocation ground-truthing. J. Forensic Sci. 62, 655–667 (2017).CAS 
PubMed 
Article 

Google Scholar 
48.Jameel, Y. et al. Tap water isotope ratios reflect urban water system structure and dynamics across a semiarid metropolitan area. Water Resour. Res. 52, 5891–5910 (2016).CAS 
Article 
ADS 

Google Scholar 
49.Lee-Thorp, J., Manning, L. & Sponheimer, M. Problems and prospects for carbon isotope analysis of very small samples of fossil tooth enamel. Bull. Geol. Soc. Fr. 168, 767–773 (1997).CAS 

Google Scholar 
50.Connan, M. et al. Multidimensional stable isotope analysis illuminates resource partitioning in a sub-Antarctic island bird community. Ecography 42, 1948–1959 (2019).Article 

Google Scholar 
51.Coplen, T. B. Guidelines and recommended terms for expression of stable-isotope-ratio and gas-ratio measurement results. Rapid Commun. Mass Spectrom. 25, 2538–2560 (2011).CAS 
PubMed 
Article 
ADS 

Google Scholar 
52.Metro Vancouver. Water Treatment & Supply  http://www.metrovancouver.org/services/water/Pages/default.aspx (2020).53.International Atomic Energy Agency. International Atomic Energy Agency: RCWIP (Regionalized Cluster-Based Water Isotope Prediction) Model: Gridded precipitation δ18O|δ2H| δ18O and δ2H isoscape data http://www.iaea.org/water (2014).54.Gujarati, D. N. Linear regression: A mathematical introduction. Rapid Commun. Mass Spectrom. https://doi.org/10.4135/9781071802571 (2019).Article 

Google Scholar 
55.Coplen, T. B. Normalization of oxygen and hydrogen isotope data. Chem. Geol. Isot. Geosci. Sect. 72, 293–297 (1988).CAS 
Article 

Google Scholar 
56.Kruskal, W. H. & Wallis, W. A. Use of ranks in one-criterion variance analysis. J. Am. Stat. Assoc. 47, 583–621 (1952).MATH 
Article 

Google Scholar 
57.Shapiro, S. S. & Wilk, M. B. An analysis of variance test for normality (complete samples). Biometrika 52, 591–611 (1965).MathSciNet 
MATH 
Article 

Google Scholar 
58.R Core Team. R: A Language and Environment for Statistical Computing (2020).59.Kobayashi, K. & Salam, M. U. Comparing simulated and measured values using mean squared deviation and its components. Agron. J. 92, 345–352 (2000).Article 

Google Scholar 
60.Gauch, H. G., Hwang, J. T. G. & Fick, G. W. Model evaluation by comparison of model-based predictions and measured values. Agron. J. 95, 1442–1446 (2003).Article 

Google Scholar 
61.Metro Vancouver. Metro Vancouver Growth Projections: A Backgrounder http://www.metrovancouver.org/services/regional-planning/PlanningPublications/OverviewofMetroVancouversMethodsinProjectingRegionalGrowth.pdf (2018).62.City of Kelowna. City of Kelowna Water http://www.kelowna.ca/cm/page393.aspx (2009).63.City of Prince George. Water and Watersheds  https://www.princegeorge.ca/CityServices/Pages/Environment/WaterAndWatersheds.aspx (2017).64.City of Winnipeg. Water and Waste Department http://www.winnipeg.ca/waterandwaste/water/default.stm (2017).65.City of Toronto. Drinking Water http://www1.toronto.ca/wps/portal/contentonly?vgnextoid=e73cfe4eda8ae310VgnVCM10000071d60f89RCRD (2017).66.Jacklin, J. Assessment of Vanderhoof South Drinking Water Supply: Source Water Characteristics https://www2.gov.bc.ca/assets/gov/environment/air-land-water/water/waterquality/water-quality-reference-documents/dwa-vanderhoof_south.pdf (2005).67.Ministère du Développement durable de l’Environnement et de la Lutte contre les changements climatiques. Summary Profile of the Rivière des Outaouais Watershed http://www.environnement.gouv.qc.ca/eau/bassinversant/bassins/outaouais/portrait-sommaire-en.pdf (2015).68.Municipalité de L’Isle-aux-Allumettes. Municipality of L’Isle-Aux-Allumettes http://www.isle-aux-allumettes.com/index-en.php (2012).69.Mant, M., Nagel, A. & Prowse, T. investigating residential history using stable hydrogen and oxygen isotopes of human hair and drinking water. J. Forensic Sci. 61, 884–891 (2016).CAS 
PubMed 
Article 

Google Scholar 
70.Koehler, G. & Hobson, K. A. Tracking cats revisited: Placing terrestrial mammalian carnivores on δ2H and δ18O isoscapes. PLoS ONE 14, e0221876 (2019).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
71.Wassenaar, L. I., Athanasopoulos, P. & Hendry, M. J. Isotope hydrology of precipitation, surface and ground waters in the Okanagan Valley, British Columbia, Canada. J. Hydrol. 411, 37–48 (2011).CAS 
Article 
ADS 

Google Scholar 
72.Cameron, E. M., Hall, G. E. M., Veizer, J. & Krouse, H. R. Isotopic and elemental hydrogeochemistry of a major river system: Fraser River, British Columbia, Canada. Chem. Geol. 122, 149–169 (1995).CAS 
Article 
ADS 

Google Scholar 
73.Waterisotope Database. http://waterisotopesDB.org (2019).74.Mostapa, R., Ishak, A. K., Mohamad, K. & Demana, R. Identification of bottled zam zam water in Malaysian market using hydrogen and oxygen stable isotope ratios(d2H and d18O). J. Nucl. Relat. Technol. 11, 48–53 (2014).
Google Scholar 
75.Thompson, A. H. et al. Stable isotope analysis of modern human hair collected from Asia (China, India, Mongolia, and Pakistan). Am. J. Phys. Anthropol. https://doi.org/10.1002/ajpa.21162 (2010).Article 
PubMed 

Google Scholar 
76.Katebe, R., Musukwa, G., Mweetwa, B., Shaba, P. & Njovu, E. Assessment of Uranium in Drinking Water in Kitwe, Chambeshi and Chingola in the Copperbelt Region of Zambia (IAEA, 2015).77.Zurbrügg, R., Wamulume, J., Kamanga, R., Wehrli, B. & Senn, D. B. River-floodplain exchange and its effects on the fluvial oxygen regime in a large tropical river system (Kafue Flats, Zambia). J. Geophys. Res. Biogeosci. 117, 1–10 (2012).Article 
CAS 

Google Scholar 
78.The Government of the Hong Kong Special Administrative Region. Water Supplies Department https://www.wsd.gov.hk/en/publications-and-statistics/pr-publications/the-facts/index.html (2019).79.Demény, A., Gugora, A. D., Kesjár, D., Lécuyer, C. & Fourel, F. Stable isotope analyses of the carbonate component of bones and teeth: The need for method standardization. J. Archaeol. Sci. 109, 104979 (2019).Article 
CAS 

Google Scholar 
80.Wunder, M. B. Using Isoscapes to Model probability surfaces for determining geographic origins BT: Isoscapes: Understanding movement, pattern, and process on Earth through isotope mapping in (eds. West, J. B., Bowen, G. J., Dawson, T. E. & Tu, K. P.) 251–270 (Springer, 2010).81.Graul, C. leafletR: Interactive Web-Maps Based on the Leaflet JavaScript Library (2016).82.Pye, K. Isotope and trace element analysis of human teeth and bones for forensic purposes. Geol. Soc. Lond. Spec. Publ. 232, 215–236 (2004).CAS 
Article 
ADS 

Google Scholar  More

1.Bo, M. et al. Characteristics of a black coral meadow in the twilight zone of the Central Mediterranean Sea. Mar. Ecol. Prog. Ser. 397, 53–61 (2009).ADS 
Article 

Google Scholar 
2.Fabri, M. et al. Megafauna of vulnerable marine ecosystems in French mediterranean submarine canyons: Spatial distribution and anthropogenic impacts. Deep Sea Res. Part II Top. Stud. Oceanogr. 104, 184–207 (2014).ADS 
Article 

Google Scholar 
3.Opresko, D. M. & Baron-Szabo, R. Re-descriptions of the antipatharian corals described by E. J. C. ESPER with selected English translations of the original German text (Cnidaria, Anthozoa, Antipatharia). Senckenb. Biol. 81, 1–21 (2019).
Google Scholar 
4.Molodstova, T. N. Deep-sea fauna of European seas: An annotated species check-list of benthic invertebrates living deeper than 2000 m in the seas bordering Europe. Antipatharia. Invertebr. Zool. 11, 3–7 (2014).Article 

Google Scholar 
5.IUCN. No Title. (2020). Available at: http://www.iucn.it/scheda.php?id=-512293580. Accessed: 27th August 2020.6.Chimienti, G., Bo, M., Taviani, M. & Mastrototaro, F. 19 Occurrence and Biogeography of Mediterranean Cold-Water Corals. (2019). https://doi.org/10.1007/978-3-319-91608-8_197.Bo, M. et al. Persistence of pristine deep-sea coral gardens in the Mediterranean Sea (SW Sardinia). PLoS ONE 10, 1–21 (2015).Article 
CAS 

Google Scholar 
8.Cau, A. et al. Leiopathes glaberrima millennial forest from SW Sardinia as nursery ground for the small spotted catshark Scyliorhinus canicula. Aquat. Conserv. Mar. Freshw. Ecosyst. 27, 731–735 (2017).Article 

Google Scholar 
9.D’Onghia, G. Cold-water corals as shelter, feeding and life-history critical habitats for fish species: Ecological interactions and fishing impact. In Mediterranean cold-water corals: Past, present and future (eds. Covadonga, O. & Jiménez, C.) 335–356 (Springer, 2019).10.Etnoyer, P. J. et al. Models of habitat suitability, size, and age-class structure for the deep-sea black coral Leiopathes glaberrima in the Gulf of Mexico. Deep. Res. Part II Top. Stud. Oceanogr. 150, 218–228 (2018).ADS 
Article 

Google Scholar 
11.Vitale, S. et al. Black coral age and growth validation using 14C dating in the central Mediterranean Sea. In 2nd international radiocarbon in the environment conference (2017).12.Carlier, A. et al. Trophic relationships in a deep mediterranean cold-water coral bank (Santa Maria di Leuca, Ionian sea). Mar. Ecol. Prog. Ser. 397, 125–137 (2009).ADS 
CAS 
Article 

Google Scholar 
13.Wagner, D., Luck, D. G. & Toonen, R. J. The biology and ecology of black corals (Cnidaria: Anthozoa: Hexacorallia: Antipatharia. Adv. Mar. Biol. 63, 67–132 (2012).PubMed 
Article 
PubMed Central 

Google Scholar 
14.Bo, M. et al. Deep coral oases in the South Tyrrhenian Sea. PLoS ONE 7, e49870 (2012).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
15.Mytilineou, C., Smith, C., Anastasopoulou, A., Papadopoulou, K. & Christidis, G. New cold-water coral occurrences in the Eastern Ionian Sea: Results from experimental long line fishing. Deep. Res. Part I Oceanogr. Res. Pap. 99, 146–157 (2014).ADS 

Google Scholar 
16.Angeletti, L. et al. First report of live deep water cnidarian assemblages from the Malta Escarpment. Ital. J. Zool. 82, 291–297 (2015).
Google Scholar 
17.Costantini, F., Fauvelot, C. & Abbiati, M. Fine-scale genetic structuring in Corallium rubrum: Evidence of inbreeding and limited effective larval dispersal. Mar. Ecol. Prog. Ser. 340, 109–119 (2007).ADS 
CAS 
Article 

Google Scholar 
18.Miller, K. Short-distance dispersal of black coral larvae:inference from spatial analysis of colony genotypes. Mar. Ecol. Prog. Ser. 171, 225–233 (1998).ADS 
Article 

Google Scholar 
19.Bo, M. et al. The coral assemblages of an off-shore deep Mediterranean rocky bank (NW Sicily, Italy). Mar. Ecol. 35, 332–342 (2014).ADS 
Article 

Google Scholar 
20.Bo, M. et al. Fishing impact on deep Mediterranean rocky habitats as revealed by ROV investigation. Biol. Conserv. 171, 167–176 (2014).Article 

Google Scholar 
21.Deidun, A., Tsounis, G., Balzan, F. & Micallef, A. Records of black coral (Antipatharia) and red coral (Corallium rubrum) fishing activities in the Maltese Islands. Mar. Bioldivers. Rec. 3, 1–6 (2010).Article 

Google Scholar 
22.Tsounis, G. et al. The exploitation and conservation of precious corals. Oceanogr. Mar. Biol. An Annu. Rev. 48, 161–211 (2010).Article 

Google Scholar 
23.Deidun, A. et al. First characterisation of a Leiopathes glaberrima (Cnidaria: Anthozoa: Antipatharia) forest in Maltese exploited fishing grounds. Ital. J. Zool. 82, 271–280 (2015).
Google Scholar 
24.Thompson, A., Sanders, J., Tandstad, M., Carocci, F. & Fuller, M. Vulnerable Marine Ecosystems: Processes and Practices in the High Seas (2016).25.Fabri, M. C. et al. Megafauna of vulnerable marine ecosystems in French mediterranean submarine canyons: Spatial distribution and anthropogenic impacts. Deep. Res. Part II Top. Stud. Oceanogr. 104, 184–207 (2014).ADS 
Article 

Google Scholar 
26.Pearson, R. G. Species distribution modeling for conservation educators and practitioners synthesis. In (2008).27.Sundahl, H., Buhl-Mortensen, P. & Buhl-Mortensen, L. Distribution and Suitable Habitat of the Cold-Water Corals Lophelia pertusa, Paragorgia arborea, and Primnoa resedaeformis on the Norwegian continental shelf. Front. Mar. Sci. 7, 1–22 (2020).Article 

Google Scholar 
28.Tittensor, D. P. et al. Predicting global habitat suitability for stony corals on seamounts. J. Biogeogr. 36, 1111–1128 (2009).Article 

Google Scholar 
29.Davies, A. J. & Guinotte, J. M. Global habitat suitability for framework-forming cold-water corals. PLoS ONE 6, e18483 (2011).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
30.Bargain, A. et al. Predictive habitat modeling in two Mediterranean canyons including hydrodynamic variables. Prog. Oceanogr. 169, 151–168 (2018).ADS 
Article 

Google Scholar 
31.Lauria, V. et al. Species distribution models of two critically endangered deep-sea octocorals reveal fishing impacts on vulnerable marine ecosystems in central Mediterranean Sea. Sci. Rep. 7, 1–14 (2017).CAS 
Article 

Google Scholar 
32.Massi, D. et al. Spatial distribution of the black coral Leiopathes glaberrima (Esper, 1788) (Antipatharia: Leiopathidae) in the Mediterranean: A prerequisite for protection of Vulnerable Marine Ecosystems (VMEs). Eur. Zool. J. 85, 170–179 (2018).Article 

Google Scholar 
33.Hayes, D. R. et al. Review of the circulation and characteristics of intermediate water masses of the mediterranean: Implications for cold-water coral habitats. In Mediterranean cold-water corals: Past, present and future (eds. Orejas, C. & Jiménez, C.) 195–212 (Springer, 2019).34.Pinardi, N. & Masetti, E. Variabilityofthelarge-scalegeneralcirculationofthe, MediterraneanSeafromobservationsandmodelling:areview. Palaeogeogr. Palaeoclimatol. Palaeoecol. 158, 153–173 (2000).Article 

Google Scholar 
35.Pinardi, N. et al. Mediterranean Sea large-scale low-frequency ocean variability and water mass formation rates from 1987 to 2007: A retrospective analysis. Prog. Oceanogr. 132, 318–332 (2015).ADS 
Article 

Google Scholar 
36.Astraldi, M. et al. Water mass properties and chemical signatures in the central Mediterranean region. J. Mar. Syst. 133–134, 155–177 (2002).Article 

Google Scholar 
37.Taviani, M. et al. The “Sardinian cold-water coral province” in the context of the Mediterranean coral ecosystems. Deep. Res. Part II Top. Stud. Oceanogr. 145, 61–78 (2017).ADS 
Article 

Google Scholar 
38.EMODnet Bathymetry Consortium. EMODnet Digital Bathymetry (DTM 2016). EMODnet Bathymetry Consortium (2016). https://doi.org/10.12770/c7b53704-999d-4721-b1a3-04ec60c87238.39.McArthur, M. et al. On the use of abiotic surrogates to describe marine benthic biodiversity. Estuar. Coast. Shelf Sci. 88, 21–32 (2010).ADS 
Article 

Google Scholar 
40.Lauria, V. et al. Species distribution models of two critically endangered deep-sea octocorals reveal fishing impacts on vulnerable marine ecosystems in central Mediterranean Sea. Sci. Rep. 7, 8049 (2017).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
41.Wilson, M. F. J., O’Connell, B., Brown, C., Guinan, J. C. & Grehan, A. J. Multiscale terrain analysis of multibeam bathymetry data for habitat mapping on the continental slope. Mar. Geod. 30, 3–35 (2007).Article 

Google Scholar 
42.Tong, R., Purser, A., Unnithan, V. & Guinan, J. Multivariate statistical analysis of distribution of deep-water gorgonian corals in relation to seabed topography on the norwegian margin. PLoS ONE 7, 1–13 (2012).
Google Scholar 
43.Savini, A., Vertino, A., Marchese, F., Beuck, L. & Freiwald, A. Mapping cold-water coral habitats at different scales within the Northern Ionian Sea (central Mediterranean): An assessment of coral coverage and associated vulnerability. PLoS ONE 9, e102405 (2014).Article 
CAS 

Google Scholar 
44.Sbrocco, E. & Barber, P. MARSPEC: Ocean climate layers for marine spatial ecology. Ecology 94, 979 (2013).Article 

Google Scholar 
45.White, M., Mohn, C., Stigter, H. & Mottram, G. Deep-water coral development as a function of hydrodynamics and surface productivity around the submarine banks of the Rockall Trough, NE Atlantic. In eds. Freiwald, A. & Robert, JM. 503−514 (2005).46.Davies, A. J. et al. Downwelling and deep-water bottom currents as food supply mechanisms to the cold-water coral Lophelia pertusa (Scleractinia) at the Mingulay Reef complex. Limnol. Oceanogr. 54, 620–629 (2009).ADS 
Article 

Google Scholar 
47.Clementi, E. et al. Mediterranean Sea Analysis and Forecast (CMEMS MED-Currents, EAS5 system). (2019).48.Guinotte, J. Climate change and deep-sea corals. J. Mar. Educ. 21, 48–49 (2005).
Google Scholar 
49.Zuur, A. F., Ieno, E. N. & Smith, G. M. Analysing Ecological Data. (Springer2007).50.Phillips, S. J., Anderson, R. P. & Schapire, R. E. Maximum entropy modelling of species distributions. Ecol. Model. 190, 231–259 (2006).Article 

Google Scholar 
51.Elith, J. et al. Novel methods improve prediction of species’ distributions from occurrence data. Ecography (Cop.) 29, 129–151 (2006).Article 

Google Scholar 
52.Tong, Y., Chen, X. & Chen, Y. Evaluating alternative management strategies for bigeye tuna, Thunnus obesus in the Indian Ocean. Sci. Mar. 77, 449–460 (2013).Article 

Google Scholar 
53.Pearson, R. G., Raxworthy, C. J., Nakamura, M. & Peterson, A. T. Predicting species distributions from small numbers of occurrence records: a test case using cryptic geckos in Madagascar. J. Biogeogr. 34, 102–117 (2007).Article 

Google Scholar 
54.Wisz, M. S. et al. Effects of sample size on the performance of species distribution models. Divers. Distrib. 14, 763–773 (2008).Article 

Google Scholar 
55.Bargain, A., Marchese, F., Savini, A., Taviani, M. & Fabri, M. C. Santa Maria di Leuca province (Mediterranean Sea): Identification of suitable mounds for cold-water coral settlement using geomorphometric proxies and maxent methods. Front. Mar. Sci. 4, 1–17 (2017).Article 

Google Scholar 
56.Fabri, M. C., Bargain, A., Pairaud, I., Pedel, L. & Taupier-Letage, I. Cold-water coral ecosystems in Cassidaigne Canyon: An assessment of their environmental living conditions. Deep. Res. Part II Top. Stud. Oceanogr. 137, 436–453 (2017).ADS 
Article 

Google Scholar 
57.Phillips, S. J. & Dudı, M. Modeling of species distributions with Maxent: New extensions and a comprehensive evaluation. 161–175 (2008). https://doi.org/10.1111/j.2007.0906-7590.05203.x58.Elith, J. et al. A statistical explanation of MaxEnt for ecologists. Divers. Distrib. 17, 43–57 (2011).Article 

Google Scholar 
59.Phillips, S. J. & Dudik, M. Modeling of species distributions withMaxent: New extensions and a comprehensive evaluation. Ecography (Cop.) 31, 161–175 (2008).Article 

Google Scholar 
60.Young, N., Lane, C. & Evangelista, P. A MaxEnt Model Tutorial (ArcGISv10) (2011).61.Phillips, S. J., Anderson, R. P., Dudík, M., Schapire, R. E. & Blair, M. Opening the black box: an open-source release of Maxent. Ecography (Cop.) 40, 887–893 (2017).Article 

Google Scholar 
62.Mytilineou, C. et al. New cold-water coral occurrences in the Eastern Ionian Sea: Results from experimental long line fishing. Deep. Res. Part II Top. Stud. Oceanogr. 99, 146–157 (2014).ADS 
Article 

Google Scholar 
63.Chimienti, G., Bo, M., Taviani, M. & Mastrototaro, F. Coral reefs of the world. In Mediterranean cold-water corals: Past; present and future (eds. Orejas, C. & Jiménez C.) 213–243 (Springer, 2019).64.Mascle, J., Migeon, S., Coste, M., Hassoun, V. & Rouillard, P. Rocky vs sedimentary canyons around the mediterranean sea and the black sea. In Submarine canyon dynamics in the mediterranean and tributary seas—An integrated geological, oceanographic and biological perspective CIESM (ed. Briand, F.) (2015).65.Canals, M. et al. Flushing submarine canyons. Nature 444, 354–357 (2006).ADS 
CAS 
PubMed 
Article 

Google Scholar 
66.Huvenne, V. A. I. et al. picture on the wall: innovative mapping reveals coldwater coral refuge in submarine canyon. PLoS ONE 6, e28755 (2011).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
67.Bo, M. et al. Coral assemblages off the Calabrian Coast (South Italy) with new observations on living colonies of Antipathes dichotoma. Ital. J. Zool. 78, 231–242 (2010).Article 

Google Scholar 
68.Freiwald, A. & Roberts, M. J. Cold-water corals and ecosystems: preface. In Cold-water corals and ecosystems (ed. Freiwald, A. R.) 1243 (Springer, 2005).69.Khripounoff, A. et al. Deep cold-water coral ecosystems in the Brittany submarine canyons (Northeast Atlantic): Hydrodynamics, particle supply, respiration, and carbon cycling. Limnol. Oceanogr. 59, 87–98 (2014).ADS 
CAS 
Article 

Google Scholar 
70.Davies, A. J., Wisshak, M., Orr, J. C. & Murray Roberts, J. Predicting suitable habitat for the cold-water coral Lophelia pertusa (Scleractinia). Deep. Res. Part I Oceanogr. Res. Pap. 55, 1048–1062 (2008).ADS 
Article 

Google Scholar 
71.Greathead, C. F., Donnan, D. W., Mair, J. M. & Saunders, G. R. The sea pens Virgularia mirabilis, Pennatula phosphorea and Funiculina quadrangularis: Distribution and conservation issues in Scottish waters. J. Mar. Biol. Assoc. UK 87, 1095–1103 (2007).Article 

Google Scholar 
72.Barry, S. & Elith, J. Error and uncertainty in habitat models. J. Appl. Ecol. 43, 413–423 (2006).Article 

Google Scholar 
73.Wisz, M. S. et al. Effects of sample size on the performance of species distribution models. Divers. Distrib. 763–773 (2008).74.Parolo, G., Rossi, G. & Ferrarini, A. oward improved species niche mod-elling: arnica montana in the Alps as a case study. J. Appl. Ecol. 45, 1410–1418 (2008).Article 

Google Scholar 
75.Gogol-Prokurat, M. Predicting habitat suitability for rare plants at localspatial scales using a species distribution model. Ecol. Appl. 21, 33–47 (2011).PubMed 
Article 

Google Scholar 
76.Bean, W. T., Stafford, R. & Brashares, J. S. he effects of small sample sizeand sample bias on threshold selection and accuracy assessment of species dis-tribution models. Ecography (Cop.) 35, 250–258 (2012).Article 

Google Scholar 
77.Morales, N. S., Fernández, I. C. & Baca-González, V. MaxEnt’s parameter configuration and small samples: are we paying attention to recommendations? A systematic review. PeerJ (2017).78.Opresko, D. M. & Försterra, G. No Title. in El Mar Mediterraneo: fauna, flora, ecologia. (ed. Hofrichter, R.) 506–509 (2004).79.Aguilar, R. & Marín, P. Mediterranean deep-sea corals : reasons protection under the Barcelona Convention. 1–18 (2013).80.FAO. International Guidelines for the Management of Deep-sea Fisheries in the High Seas. 73 (2009).81.Marin, P. & Aguilar, R. Mediterranean submarine canyons 2012: pending protection. in Mediterranean Submarine Canyons: Ecology and Governance 191–206 (IUCN, 2012).82.Thurstan, R. H., Brockington, S. & Roberts, C. M. The effects of 118 years of industrial fishing on UK bottom trawl fisheries. Nat. Commun. 1, 1–6 (2010).CAS 
Article 

Google Scholar 
83.Worm, B. & Tittensor, D. P. Range contraction in large pelagic predators. Proc. Natl. Acad. Sci. USA 108, 11942–11947 (2011).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
84.Gross, M. Deep sea in deep trouble?. Curr. Biol. 25, 1019-R1021 (2015).Article 
CAS 

Google Scholar 
85.Morato, T., Watson, R., Pitcher, T. & Pauly, D. Fishing down the deep. Fish Fish. 7, 24–34 (2006).Article 

Google Scholar 
86.Roberts, C. M. Deep impact: The rising toll of fishing in the deep sea. Trends Ecol. Evol. 17, 242–245 (2002).MathSciNet 
Article 

Google Scholar  More

Sea urchinsSea urchins (~ 2 cm of test diameter) were transported from Dalian Haibao Fishery Company (121° 22′ E, 38° 77′ N) to the Key Laboratory of Mariculture and Stock Enhancement in North China’s Sea, Ministry of Agriculture and Rural Affairs (121° 56′ E, 38° 87′ N) in November 2020 and maintained in the fiberglass tank of 139 L (750 × 430 × 430 mm) of the recirculating system (Huixin Co., China) at ~ 14 °C. During the experiment period, the salinity of the recirculated seawater was 30.07–30.25‰. Sea urchins were fed the kelp Saccharina japonica ad libitum with aeration before the experiments. The seawater was changed every three days to remove feces and algal debris.Conspecific alarm cuesConspecific alarm cues were made from crushing one M. nudus in 50 mL of fresh seawater and then filtered using a fine silk net (260 μm of mesh size)6. Fresh conspecific alarm cues were prepared before each behavioral experiment. Five mL of conspecific alarm cues was used to simulate the natural conditions, in which a fraction of body fluids from a sea urchin instantly released into the surrounding water6.Experiment 1: whether M. nudus are attracted by the kelp exposed to conspecific alarm cuesA piece of wild fresh S. japonica (~ 5 g) was placed on the side of the device with the raceways (70 × 6 × 5 cm) and 5 mL of conspecific alarm cues were subsequently added above it in the experimental group (Fig. 4A). Five mL of conspecific alarm cues were added in the control group, while the kelp was not involved (Fig. 4B). The sea urchin was individually placed 15 cm away from the kelp for each group (Fig. 4A,B). Foraging behavior of sea urchins was recorded for 20 min using a digital video recorder (FDR-AXP55, SONY, Japan). The danger area (b area) refers to the position 15 cm away from the kelp and the escape area (a area) refers to the position 15 cm away from the sea urchin (Fig. 4A,B). The time of sea urchins spent in the danger (a area) and escape (b area) areas was calculated individually for each group. The experiment was repeated 20 times using different sea urchins for each group (N = 20). The seawater was changed and the experimental device was washed for each trial to avoid potential non-experimental effects.Figure 4Foraging devices for the experiments with Mesocentrotus nudus. Experiment 1: the kelp was attracted by sea urchins (A,B); the positions of a, c and f refer to the escape areas, while positions of b, d and e refer to the danger areas. Experiment 2: five mL of conspecific alarm cues were put in the middle of the device (C). Experiment 3 (D–F): foraging behavior of M. nudus at 70 cm in the first trial (D); five mL of conspecific alarm cues were put in the middle of the device in the second trial (E); the position of g refers to the danger area; different amounts of conspecific alarm cues (5 mL and 0.5 mL) were put in the left and right of the devices in the third trial (F); the positions of h and i refer to the safety areas close to the areas with more and less conspecific alarm cues, respectively. This figure was created using the Adobe Photoshop (version CS5) software.Full size imageExperiment 2: whether foraging behavior of fasted M. nudus is affected when they encountered conspecific alarm cuesThis experiment investigated the effects of conspecific alarm cues on the foraging behavior of fasted sea urchins. The experiment group was set to simulate that conspecific alarm cues appear on the way to kelp beds (Fig. 4C). Sea urchins were fasted for 7, 14 and 21 days in the experiment groups.An acrylic device with the raceways (70 × 6 × 5 cm) was designed according to the method of our previous study17 with some revisions. Five mL of conspecific alarm cues were added in the middle of the device, the wild fresh S. japonica (~ 10 g) was placed 15 cm away from the left side of the device, while the sea urchin was placed on 15 cm away from the right of the device (Fig. 4C). Foraging behavior of sea urchins was recorded for 20 min using a digital video recorder (FDR-AXP55, SONY, Japan). Arriving at the kelp within 20 min was defined as successfully foraging18. The danger area (e area) refers to the position between conspecific alarm cues and the sea urchin, while the escape area (f area) refers to the 15 cm position away from the right side of the sea urchin (Fig. 4C). The time of sea urchins spent in the e and f areas was calculated individually for all the groups. The individual experiment was repeated 20 times with different sea urchins for each group (N = 20).Experiment 3: whether foraging ability of M. nudus links to their responses to conspecific alarm cuesAn acrylic device with the raceways (70 × 6 × 5 cm) was designed for the measurement of foraging behavior of sea urchins. Experiment 3 included three trials.In the first trial, we measured the foraging ability of sea urchins. Ten grams of the kelp was placed on one side of the tank and one sea urchin was placed on the other side of the tank (Fig. 4D). Foraging behavior was recorded for 20 min using a digital video recorder (FDR-AXP55, SONY, Japan). Movement and velocity of sea urchins were calculated using the software ImageJ (version 1.51 n). The experiment was individually repeated 31 times using different sea urchins (N = 31). At the end of this foraging experiment, 31 sea urchins were recorded in individual cylindrical plastic cages in the tank. All the sea urchins were subsequently measured for the second and third trials (N = 31). To avoid potential fatigue influences, each behavioral experiment was carried out after another 24 h.The second trial was carried out to test the hypothesis that sea urchins with strong foraging ability are affected when conspecific alarm cues appear on the way to the kelp. A piece of wild fresh kelp (~ 10 g) was placed on the left side of the device (70 × 6 × 5 cm), while 5 mL of conspecific alarm cues were put in the middle of the device. The individual was placed on the right of the device (Fig. 4E). Foraging behavior of sea urchins was recorded for 20 min using a digital video recorder (FDR-AXP55, SONY, Japan). The danger area (g area) refers to the 15 cm position away from conspecific alarm cues. The time of the 31 sea urchins spent in the g area was individually calculated.We investigated whether different amounts of conspecific alarm cues around the kelp affected the sea urchins with strong foraging ability in the third trial. More conspecific alarm cues (5 mL) and a piece of wild fresh kelp (~ 10 g) were placed on the left side of the device (70 × 6 × 5 cm), while less conspecific alarm cues (0.5 mL) and wild fresh kelp (~ 10 g) were placed on the right side of the device (Fig. 4F). Sea urchins were individually placed in the middle of the device. Foraging behavior of sea urchins was recorded for 20 min using a digital video recorder (FDR-AXP55, SONY, Japan). The number of sea urchins that successfully foraged to areas with more and less conspecific alarm cues was recorded. The area h refers to the safety area that was 15 cm position away from the sea urchin, which was close to the side with more conspecific alarm cues. The area of i refers to the safety area that was 15 cm position away from the sea urchin, which was close to the side with less conspecifics alarm cues. The time spent in the h and i areas was individually calculated in 31 sea urchins.Statistical analysisNormal distribution and homogeneity of variance of the data were analyzed using the Kolmogorov–Smirnov test and Levene test, respectively. Independent-sample t-test was carried out to compare the differences of duration in danger (b and d) and escape (a and c) areas in experiment 1. One-way ANOVA was used to analyze the duration in danger (e) and escape (f) areas of the sea urchins fasted for 7, 14 and 21 days. Pairwise multiple comparisons were carried out using LSD test when significant differences were found in the ANOVAs. Correlations between foraging velocity and duration in the g, h and i areas were analyzed using Pearson correlation analysis. A probability level of P < 0.05 was considered significant. All data analyses were performed using SPSS 21.0 statistical software. Graphs 2−4 were performed using Origin 9.0 software (OriginLab, USA). More