Ecology

LookingGlass design and optimizationDataset generationThe taxonomic organization of representative Bacterial and Archaeal genomes was determined from the Genome Taxonomy Database, GTDB51 (release 89.0). The complete genome sequences were downloaded via the NCBI Genbank ftp52. This resulted in 24,706 genomes, comprising 23,458 Bacterial and 1248 Archaeal genomes.Each genome was split into read-length chunks. To determine the distribution of realistic read lengths produced by next-generation short-read sequencing machines, we obtained the BioSample IDs52 for each genome, where they existed, and downloaded their sequencing metadata from the MetaSeek53 database using the MetaSeek API. We excluded samples with average read lengths less than 60 or greater than 300 base pairs. This procedure resulted in 7909 BioSample IDs. The average read lengths for these sequencing samples produced the read-length distribution (Supplementary Fig. 1) with a mean read length of 136 bp. Each genome was split into read-length chunks (with zero overlap in order to maximize information density and reduce data redundancy in the dataset): a sequence length was randomly selected with replacement from the read-length distribution and a sequence fragment of that length was subset from the genome, with a 50% chance that the reverse complement was used. The next sequence fragment was chosen from the genome starting at the end point of the previous read-length chunk, using a new randomly selected read length, and so on. These data were partitioned into a training set used for optimization of the model; a validation set used to evaluate model performance during parameter tuning and as a benchmark to avoid overfitting during training; and a test set used for final evaluation of model performance. To ensure that genomes in the training, validation, and test sets had low sequence similarity, the sets were split along taxonomic branches such that genomes from the Actinomycetales, Rhodobacterales, Thermoplasmata, and Bathyarchaeia were partitioned into the validation set; genomes from the Bacteroidales, Rhizobiales, Methanosarcinales, and Nitrososphaerales were partitioned into the test set; and the remaining genomes remained in the training set. This resulted in 529,578,444 sequences in the training set, 57,977,217 sequences in the validation set, and 66,185,518 sequences in the test set. We term this set of reads the GTDB representative set (Table 1).Table 1 Summary table of datasets used.Full size tableThe amount of data needed for training was also evaluated (Supplementary Fig. 2). Progressively larger amounts of data were tested by selecting at random 1, 10, 100, or 500 read-length chunks from each of the GTDB representative genomes in the GTDB representative training set. Additionally, the performance of smaller but more carefully selected datasets, representing the diversity of the microbial tree of life, were tested by selecting for training one genome at random from each taxonomic class or order in the GTDB taxonomy tree. In general, better accuracy was achieved in fewer epochs with a greater amount of sequencing data (Supplementary Fig. 2); however, a much smaller amount of data performed better if a representative genome was selected from each GTDB taxonomy class.The final LookingGlass model was trained on this class-level partition of the microbial tree of life. We term this dataset the GTDB class set (Table 1). The training, validation, and test sets were split such that no classes overlapped across sets: the validation set included 8 genomes from each of the classes Actinobacteria, Alphaproteobacteria, Thermoplasmata, and Bathyarchaeia (32 total genomes); the test set included 8 genomes from each of the classes Bacteroidia, Clostridia, Methanosarcinia, and Nitrososphaeria (32 total genomes); and the training set included 1 genome from each of the remaining classes (32 archaeal genomes and 298 bacterial genomes for a total of 330 genomes). This resulted in a total of 6,641,723 read-length sequences in the training set, 949,511 in the validation set, and 632,388 in the test set (Supplementary Data 1).Architecture design and trainingRecurrent neural networks (RNNs) are a type of neural network designed to take advantage of the context dependence of sequential data (such as text, video, audio, or biological sequences), by passing information from previous items in a sequence to the current item in a sequence54. Long short-term memory networks (LSTMs)55 are an extension of RNNs, which better learn long-term dependencies by handling the RNN tendency to “forget” information farther away in a sequence56. LSTMs maintain a cell state which contains the “memory” of the information in the previous items in the sequence. LSTMs learn additional parameters which decide at each step in the sequence which information in the cell state to “forget” or “update”.LookingGlass uses a three-layer LSTM encoder model with 1152 units in each hidden layer and an embedding size of 104 based on the results of hyperparameter tuning (see below). It divides the sequence into characters using a kmer size of 1 and a stride of 1, i.e., is a character-level language model. LookingGlass is trained in a self-supervised manner to predict a masked nucleotide, given the context of the preceding nucleotides in the sequence. For each read in the training sequence, multiple training inputs are considered, shifting the nucleotide that is masked along the length of the sequence from the second position to the final position in the sequence. Because it is a character-level model, a linear decoder predicts the next nucleotide in the sequence from the possible vocabulary items “A”, “C”, “G”, and “T”, with special tokens for “beginning of read”, “unknown nucleotide” (for the case of ambiguous sequences), “end of read” (only “beginning of read” was tokenized during LookingGlass training), and a “padding” token (used for classification only).Regularization and optimization of LSTMs require special approaches to dropout and gradient descent for best performance57. The fastai library58 offers default implementations of these approaches for natural language text, and so we adopt the fastai library for all training presented in this paper. We provide the open source fastBio python package59 which extends the fastai library for use with biological sequences.LookingGlass was trained on a Pascal P100 GPU with 16GB memory on Microsoft Azure, using a batch size of 512, a back propagation through time (bptt) window of 100 base pairs, the Adam optimizer60, and utilizing a Cross Entropy loss function (Supplementary Table 1). Dropout was applied at variable rates across the model (Supplementary Table 1). LookingGlass was trained for a total of 12 days for 75 epochs, with progressively decreasing learning rates based on the results of hyperparameter optimization (see below): for 15 epochs at a learning rate of 1e−2, for 15 epochs at a learning rate of 2e−3, and for 45 epochs at a learning rate of 1e−3.Hyperparameter optimizationHyperparameters used for the final training of LookingGlass were tuned using a randomized search of hyperparameter settings. The tuned hyperparameters included kmer size, stride, number of LSTM layers, number of hidden nodes per layer, dropout rate, weight decay, momentum, embedding size, bptt size, learning rate, and batch size. An abbreviated dataset consisting of ten randomly selected read-length chunks from the GTDB representative set was created for testing many parameter settings rapidly. A language model was trained for two epochs for each randomly selected hyperparameter combination, and those conditions with the maximum performance were accepted. The hyperparameter combinations tested and the selected settings are described in the associated Github repository61.LookingGlass validation and analysis of embeddingsFunctional relevanceDataset generation
In order to assess the ability of the LookingGlass embeddings to inform the molecular function of sequences, metagenomic sequences from a diverse set of environments were downloaded from the Sequence Read Archive (SRA)62. We used MetaSeek53 to choose ten metagenomes at random from each of the environmental packages defined by the MIxS metadata standards63: built environment, host-associated, human gut, microbial mat/biofilm, miscellaneous, plant-associated, sediment, soil, wastewater/sludge, and water, for a total of 100 metagenomes. The SRA IDs used are available in (Supplementary Table 2). The raw DNA reads for these 100 metagenomes were downloaded from the SRA with the NCBI e-utilities. These 100 metagenomes were annotated with the mi-faser tool27 with the read-map option to generate predicted functional annotation labels (to the fourth digit of the Enzyme Commission (EC) number), out of 1247 possible EC labels, for each annotatable read in each metagenome. These reads were then split 80%/20% into training/validation candidate sets of reads. To ensure that there was minimal overlap in sequence similarity between the training and validation set, we compared the validation candidate sets of each EC annotation to the training set for that EC number with CD-HIT64, and filtered out any reads with >80% DNA sequence similarity to the reads of that EC number in the training set (the minimum CD-HIT DNA sequence similarity cutoff). In order to balance EC classes in the training set, overrepresented ECs in the training set were downsampled to the mean count of read annotations (52,353 reads) before filtering with CD-HIT. After CD-HIT processing, any underrepresented EC numbers in the training set were oversampled to the mean count of read annotations (52,353 reads). The validation set was left unbalanced to retain a distribution more realistic to environmental settings. The final training set contained 61,378,672 reads, while the validation set contained 2,706,869 reads. We term this set of reads and their annotations the mi-faser functional set (Table 1).
As an external test set, we used a smaller number of DNA sequences from genes with experimentally validated molecular functions. We linked the manually curated entries of Bacterial or Archaeal proteins from the Swiss-Prot database65 corresponding to the 1247 EC labels in the mi-faser functional set with their corresponding genes in the EMBL database66. We downloaded the DNA sequences, and selected ten read-length chunks at random per CDS. This resulted in 1,414,342 read-length sequences in the test set. We term this set of reads and their annotations the Swiss-Prot functional set (Table 1).

Fine-tuning procedure
We fine-tuned the LookingGlass language model to predict the functional annotation of DNA reads, to demonstrate the speed with which an accurate model can be trained using our pretrained LookingGlass language model. The architecture of the model retained the 3-layer LSTM encoder and the weights of the LookingGlass language model encoder, but replaced the language model decoder with a new multiclass classification layer with pooling (with randomly initialized weights). This pooling classification layer is a sequential model consisting of the following layers: a layer concatenating the output of the LookingGlass encoder with min, max, and average pooling of the outputs (for a total dimension of 104*3 = 312), a batch normalization67 layer with dropout, a linear layer taking the 312-dimensional output of the batch norm layer and producing a 50-dimensional output, another batch normalization layer with dropout, and finally a linear classification layer that is passed through the log(Softmax(x)) function to output the predicted functional annotation of a read as a probability distribution of the 1247 possible mi-faser EC annotation labels. We then trained the functional classifier on the mi-faser functional set described above. Because the >61 million reads in the training set were too many to fit into memory, training was done in 13 chunks of ~5-million reads each until one total epoch was completed. Hyperparameter settings for the functional classifier training are seen in Supplementary Table 1.

Encoder embeddings and MANOVA test
To test whether the LookingGlass language model embeddings (before fine-tuning, above) are distinct across functional annotations, a random subset of ten reads per functional annotation was selected from each of the 100 SRA metagenomes (or the maximum number of reads present in that metagenome for that annotation, whichever was greater). This also ensured that reads were evenly distributed across environments. The corresponding fixed-length embedding vectors for each read was produced by saving the output from the LookingGlass encoder (before the embedding vector is passed to the language model decoder) for the final nucleotide in the sequence. This vector represents a contextually relevant embedding for the overall sequence. The statistical significance of the difference between embedding vectors across all functional annotation groups was tested with a MANOVA test using the R stats package68.
Evolutionary relevance
Dataset generation
The OrthoDB database69 provides orthologous groups (OGs) of proteins at various levels of taxonomic distance. For instance, the OrthoDB group “77at2284” corresponds to proteins belonging to “Glucan 1,3-alpha-glucosidase at the Sulfolobus level”, where “2284” is the NCBI taxonomy ID for the genus Sulfolobus.
We tested whether embedding similarity of homologous sequences (sequences within the same OG) is higher than that of nonhomologous sequences (sequences from different OGs). We tested this in OGs at multiple levels of taxonomic distance—genus, family, order, class, and phylum. At each taxonomic level, ten individual taxa at that level were chosen from across the prokaryotic tree of life (e.g., for the genus level, Acinetobacter, Enterococcus, Methanosarcina, Pseudomonas, Sulfolobus, Bacillus, Lactobacillus, Mycobacterium, Streptomyces, and Thermococcus were chosen). For each taxon, 1000 randomly selected OGs corresponding to that taxon were chosen; for each of these OGs, five randomly chosen genes within this OG were chosen.
OrthoDB cross-references OGs to UniProt65 IDs of the corresponding proteins. We mapped these to the corresponding EMBL CDS IDs66 via the UniProt database API65; DNA sequences of these EMBL CDSs were downloaded via the EMBL database API. For each of these sequences, we generated LookingGlass embedding vectors.

Homologous and nonhomologous sequence pairs
To create a balanced dataset of homologous and nonhomologous sequence pairs, we compared all homologous pairs of the five sequences in an OG (total of ten homologous pairs) to an equal number of randomly selected out-of-OG comparisons for the same sequences; i.e., each of the five OG sequences was compared to 2 other randomly selected sequences from any other randomly selected OG (total of ten nonhomologous pairs). We term this set of sequences, and their corresponding LookingGlass embeddings, the OG homolog set (Table 1).

Embedding and sequence similarity
For each sequence pair, the sequence and embedding similarity were determined. The embedding similarity was calculated as the cosine similarity between embedding vectors. The sequence similarity was calculated as the Smith-Waterman alignment score using the BioPython70 pairwise2 package, with a gap open penalty of −10 and a gap extension penalty of −1. The IDs of chosen OGs, the cosine similarities of the embedding vectors, and sequence similarities of the DNA sequences are available in the associated Github repository61.

Comparison to HMM-based domain searches for distant homology detection
Distantly related homologous sequences that share, e.g., Pfam71, domains can be identified using HMM-based search methods. We used hmmscan25 (e-val threshold = 1e−10) to compare homologous (at the phylum level) sequences in the OG homolog set, for which the alignment score was less than 50 bits and the embedding similarity was greater than 0.62 (total: 21,376 gene pairs). Specifically, we identified Pfam domains in each sequence and compared whether the most significant (lowest e-value) domain for each sequence was identified in common for each homologous pair.
Environmental relevance
Encoder embeddings and MANOVA test
The LookingGlass embeddings and the environment of origin for each read in the mi-faser functional set were used to test the significance of the difference between the embedding vectors across environmental contexts. The statistical significance of this difference was evaluated with a MANOVA test using the R stats package68.
Oxidoreductase classifier
Dataset generation
The manually curated, reviewed entries of the Swiss-Prot database65 were downloaded (June 2, 2020). Of these, 23,653 entries were oxidoreductases (EC number 1.-.-.-) of Archaeal or Bacterial origin (988 unique ECs). We mapped their UniProt IDs to both their EMBL CDS IDs and their UniRef50 IDs via the UniProt database mapper API. Uniref50 IDs identify clusters of sequences with >50% amino acid identity. This cross-reference identified 28,149 EMBL CDS IDs corresponding to prokaryotic oxidoreductases, belonging to 5451 unique UniRef50 clusters. We split this data into training, validation, and test sets such that each UniRef50 cluster was contained in only one of the sets, i.e., there was no overlap in EMBL CDS IDs corresponding to the same UniRef50 cluster across sets. This ensures that the oxidoreductase sequences in the validation and test sets are dissimilar to those seen during training. The DNA sequences for each EMBL CDS ID were downloaded via the EMBL database API. These data generation process were repeated for a random selection of non-oxidoreductase UniRef50 clusters, which resulted in 28,149 non-oxidoreductase EMBL CDS IDs from 13,248 unique UniRef50 clusters.
Approximately 50 nucleotide read-length chunks (selected from the representative read-length distribution, as above) were selected from each EMBL CDS DNA sequence, with randomly selected start positions on the gene and a 50% chance of selecting the reverse complement, such that an even number of read-length sequences with “oxidoreductase” and “not oxidoreductase” labels were generated for the final dataset. This procedure produced a balanced dataset with 2,372,200 read-length sequences in the training set, 279,200 sequences in the validation set, and 141,801 sequences in the test set. We term this set of reads and their annotations the oxidoreductase model set (Table 1). In order to compare the oxidoreductase classifier performance to an HMM-based method, reads with “oxidoreductase” labels in the oxidoreductase model test set (71,451 reads) were 6-frame translated and searched against the Swiss-Prot protein database using phmmer25 (reporting e-val threshold = 0.05, using all other defaults).

Fine-tuning procedure
Since our functional annotation classifier addresses a closer classification task to the oxidoreductase classifier than LookingGlass itself, the architecture of the oxidoreductase classifier was fine-tuned starting from the functional annotation classifier, replacing the decoder with a new pooling classification layer (as described above for the functional annotation classifier) and with a final output size of 2 to predict “oxidoreductase” or “not oxidoreductase”. Fine tuning of the oxidoreductase classifier layers was done successively, training later layers in isolation and then progressively including earlier layers into training, using discriminative learning rates ranging from 1e−2 to 5e−4, as previously described72. The fine-tuned model was trained for 30 epochs, over 18 h, on a single P100 GPU node with 16GB memory.

Model performance in metagenomes
Sixteen marine metagenomes from the surface (SRF, ~5 meters) and mesopelagic (MES, 175–800 meters) from eight stations sampled as part of the TARA expedition37 were downloaded from the SRA62 (Supplementary Table 3, SRA accession numbers ERR598981, ERR599063, ERR599115, ERR599052, ERR599020, ERR599039, ERR599076, ERR598989, ERR599048, ERR599105, ERR598964, ERR598963, ERR599125, ERR599176, ERR3589593, and ERR3589586). Metagenomes were chosen from a latitudinal gradient spanning polar, temperate, and tropical regions and ranging from −62 to 76 degrees latitude. Mesopelagic depths from four out of the eight stations were sampled from oxygen minimum zones (OMZs, where oxygen More

An array of initiatives are underway to compile reference-grade genome assemblies of life on Earth. Such assemblies can shed light on many aspects of biodiversity. As Hogg says, a reference genome helps scientists determine if a sequence is a gene, to see what it encodes and assess if there is diversity at that gene. Conservation biologists might decide to move a population to improve gene flow. When one population clears a disease quicker than another, “we can move animals with the specific genetic variant that helps deal with disease.” Unfortunately, most characteristics are polygenic, she says, but “in conservation we aim to maintain and promote as much genetic diversity as we can.” Reference genomes, she says, provide a “blueprint of life” and help researchers understand how species interact with their often rapidly changing environment.A consortium has assembled the kākāpō reference genome, and Urban has been part of the team compiling one for the takahē. It involves the Takahē Recovery team, the DOC, a team at Rockefeller University and Māori members. A high-quality takahē genome can inform all the downstream conservation efforts for this species, says Urban. It was challenging to get the right kind of samples in adequate quality, she says, “but it was totally worth it because it told us a lot about the actual genomic architecture of the takahē.”Takahē genomic information has been a crucial help in developing a computational method to assemble haplotype-resolved genomes when no parental data are available, which could prove helpful in many areas of biology. The quality of this phasing, says Urban, is comparable to that of one that involved parents’ genomes. The method combines two types of genomic information: HiFi reads from Pacific Biosciences instruments and Hi-C chromatin interaction data. Pacific Biosciences introduced circular consensus sequencing a few years ago, which builds consensus reads, or HiFi reads, from multiple passes over a DNA molecule.The computational genome assembly method hifiasm has been extended. HiFi reads and Hi-C data are combined into a graph assembly that ultimate leads to haplotype-resolved assembly of diploid genomes for which parental data are lacking. Credit: Adapted with permission from ref. 5.In developing this method, Heng Li at the Dana-Farber Cancer Institute, colleagues at University of Otago in New Zealand including Lara Urban and Neil Gemmel, and several teams from other US institutions such as Rockefeller University’s Vertebrate Genome Project and the Center for Species Survival at the National Zoo, used data from the takahē and other animals, such as the critically endangered black rhinoceros.When handling diploid and polyploid genomes, many long-read assembly tools collapse differing homologous haplotypes into a ‘consensus assembly’. Some tools avoid erasing heterozygous differences and phase genomic regions with low levels of heterozygosity, and then build contiguous sequence by stitching these blocks together. The final assembly tends to include those phased blocks as an ‘alternate assembly’.With a method called trio-binning, which uses data from individuals and their parents, scientists can obtain a haplotype-resolved assembly with two sets of contiguous sequence: two haploid genomes. Other methods draw on additional data, such as chromatin interaction data from Hi-C or Strand-Seq, which applies single-cell sequencing and resolves homologs within a cell. In Strand-Seq, only the DNA template strand used during DNA replication is sequenced.Li and colleagues developed the hifiasm algorithm5 to address complications they saw in this area, such as lengthy computational pipelines. Hifiasm applies string overlap graphs, which represent different paths along the assembled genomes. In a hifiasm graph, each node is a contiguous sequence put together from ‘phased’ HiFi reads. Li and colleagues have extended hifiasm to combine HiFi reads and Hi-C data6. First, hifiasm produces a phased assembly graph onto which Hi-C reads are mapped. The graph is made up of ‘unitigs’, contiguous sequence from heterozygous and from homozygous regions. Read coverage can be used to distinguish the two. Hifiasm further processes unitigs to build a haplotype-resolved assembly of a diploid organism.The method avoids the traditional consensus assembly approach for a diploid sample, in which half of sequences are randomly discarded, and it mixes sequences from parents, which is clearly not ideal, says Li. With people, parental data can be hard to obtain and ethical approval is needed. Meanwhile, with samples obtained from animals in the wild, as in biodiversity studies, scientists usually have little or no way to locate parents. Methods exists for haplotype-resolved assembly without parent data, but they have only been tested on human samples, he says. “Making a haplotype-resolved assembler robust to various species is a lot more challenging,” says Li. An algorithm designed for species of low sequence diversity, such as humans, may not work well for species of high diversity, such as insects. “Then there are species with mixed sequence diversity, which demands an algorithm can smoothly work with all these cases without users’ intervention,“ he says. This motivated the team to extend hifiasm.There are around 440 individual South Island takahē (Porphyrio hochstetteri) left. High-quality assemblies of the species’ genome—parents and offspring—were used to benchmark a new computational tool.
Credit: I. WarrenThe takahē data from parents and chicks helped the researchers build a haplotype-resolved assembly that was a benchmark for their computational tool. “It is critical to have trio data as the ground truth,” says Li. Instead of using human ‘trios’, they wanted to develop a robust algorithm that works for various diploid samples. Says Li, “Lara’s data is invaluable.”The approach is applicable to many species, he says, but users should remember that the genomes of different species can vary dramatically in size, sequence diversity and repetitive sequence sections. “Although we have tried hard to make hifiasm work for various species, we may have overlooked cases or properties special to certain genomes,” he says. He recommends that researchers also evaluate their assemblies carefully based on what they know about the organisms they study. Users can raise a github issue or contact him and colleagues if they can’t resolve something on their own. “We are still learning how to build better assemblies,” he says, and assembly algorithms keep evolving as data quality improves.Whenua Hou, an island off New Zealand’s South Island, is a refuge for kākāpō, a critically endangered bird species.
Credit: L. Urban More

The 1986 accident at the nuclear power plant near Chernobyl in what is now Ukraine caused the largest release of radioactivity in human history. When invading Russian troops took control of the surrounding area in the province of Kyiv Oblast in February, they destroyed important research laboratories in the partially abandoned city of Chernobyl before retreating a month later.
Competing Interests
The authors declare no competing interests. More

In the current study we utilize an extensive network and data from citizen science in order to test for among taxa variation in biases and value of information (VoI) in image recognition training data. We use data from the Norwegian Species Observation Service as an example dataset due to the generic nature of this citizen science platform, where all multicellular taxa from any Norwegian region can be reported both with and without images. The platform is open to anyone willing to report under their full real name, and does not record users’ expertise or profession. The platform had 6,205 active contributors in 2021 out of its 17,655 registered users, and currently publishes almost 27 million observations through GBIF, of which 1.08 million with one or more images. Observations have been bulk-verified by experts appointed by biological societies receiving funding for this task, with particular focus on red listed species, invasive alien species, and observations out of range or season. Observations containing pictures receive additional scrutiny, as other users can alert reporters and validators to possible mistaken identifications. An advantage of this particular platform is that no image recognition model has been integrated. This ensures that the models trained in this experiment are not trained on the output resulting from the use of any model, but with identifications and taxonomic biases springing from the knowledge and interest of human observers. Moreover, the platform’s compliance with the authoritative Norwegian taxonomy allows for analyses on taxonomic coverage.In an exploration procedure we determined the taxonomic level of orders to be suitable examples of taxa with a sufficiently wide taxonomic diversity, and enough data in the dataset to be evaluated for models in this experiment. Data collection was done by acquiring taxon statistics and observation data from the Global Biodiversity Information Facility (GBIF), the largest aggregator of biodiversity observations in the world37 for the selected orders, as well as the classes used by Troudet et al.5. The authoritative taxonomy for Norway was downloaded from the Norwegian Biodiversity Information Centre38. In the experimental procedure, models were trained for 12 distinct orders (listed in Fig. 4), artificially restricting these models to different amounts of data. In the data analysis stage, model performances relative to the amount of training data were fitted for each order, allowing the estimation of a VoI. Using the number of observations per species on GBIF, and the number of species known to be present in Norway from the Norwegian Species Nomenclature Database, we calculated relative taxonomic biases.ExplorationInitial pilot runs were done on 8 taxa (see Supplementary Information), using different subset sizes of observations for each species, and training using both an Inception-ResNet-v239 as well as an EfficientNetB340 architecture for each of these subsets. These initial results indicated that the Inception-ResNet-v2 performance (F(_1)) varied less between replicate runs and was generally higher, so subsequent experiments were done using this architecture. The number of observations which still improved the accuracy of the model was found to be between 150 and 200 in the most extreme cases, so the availability of at least 220 observations with images per species was chosen as an inclusion criteria for the further experiment. This enabled us to set aside at least 20 observations per species as a test dataset for independent model analysis.From a Darwin Core Archive file of Norwegian citizen science observations from the Species Observation Service with at least one image33, a tally of the number of such observations per species was generated. We then calculated how many species, with a minimum of 220 such observations, would, at a minimum, be available per taxon if a grouping was made based on each taxon rank level with the constraint of resulting in at least 12 distinct taxa. For each taxonomic level, we calculated how many species having at least 220 such observations were available per taxon when dividing species based on that taxon level. When deciding on the appropriate taxon level to use, we limited the options to taxon levels resulting in at least 12 different taxa.A division by order was found to provide the highest minimum number of species (17) per order within these constraints, covering 12 of the 96 eligible orders. The next best alternative was the family level, which would contain 15 species per family, covering 12 of the 267 eligible families.Data collectionWe retrieved the number of species represented in the Norwegian data through the GBIF API, for all observations, all citizen science observations, and all citizen science observations with images for the 12 selected orders and the classes used by Troudet et al.5. We also downloaded the Norwegian Species Nomenclature Database38 for all kingdoms containing taxa included in these datasets. Observations with images were collected from the Darwin Core Archive file used in the exploration phase, filtering on the selected orders. For these orders, all images were downloaded and stored locally. The average number of images per observation in this dataset was 1.44, with a maximum of 17 and a median of 1.Experimental procedureFor each selected order, a list of all species with at least 220 observations with images was generated from the Darwin Core Archive file33. Then, runs were generated according to the following protocol (Fig. 5):Figure 5Data selection and subdivision. Each run is generated by selecting 17 taxonomically adjacent species per order, and randomly assigning all available images of each selected species to that run’s test-, train- or validation set. Training data are used as input during training, using the validation data to evaluate performance after each training round in order to adjust training parameters during training. The test set is used to measure model performance independently after the model is finalized28. For each subsequent model in that run, training and validation data are reduced by 25% (or slightly less than 25% if not divisible by 4). The test set is not reduced, and used for all models within a run.Full size image

1.

From a list sorted alphabetically by the full taxonomy of the species, a subset of 17 consecutive species starting from a random index was selected. If the end of the list was reached with fewer than 17 species selected, selection continued from the start of the list. The taxonomic sorting ensures that closely related species (belonging to the same family or genus), bearing more similarity, are more likely to be part of the same experimental set. This ensures that the classification task is not simplified for taxa with many eligible species.

2.

Each of the 220+ observations for each species were tagged as being either test, training or validation data. A random subset of all but 200 were assigned to the test set. The remaining 200 observations were, in a 9:1 ratio, randomly designated as training or validation data, respectively. In all cases, images from the same observation were assigned to the same subset, to keep the information in each subset independent from the others. The resulting lists of images are stored as the test set and 200-observation task.

3.

The 200 observations in the training and validation sets were then repeatedly reduced by discarding a random subset of 25% of both, maintaining a validation data proportion of (le)10%. The resulting set was saved as the next task, and this step was repeated as long as the resulting task contained a minimum of 10 observations per species. The test set remained unaltered throughout.

Following this protocol results in a single run of related training tasks with 200, 150, 113, 85, 64, 48, 36, 27, 21, 16 and 12 observations for training and validation per species. The seeds for the randomization for both the selection of the species and for the subsetting of training- and validation datasets were stored for reproducibility. The generation of runs was repeated 5 times per order to generate runs containing tasks with different species subsets and different observation subsetting.Then, a Convolutional Neural Network based on Inception-ResNet-v239 (see the Supplementary Information for model configuration) was trained using each predesignated training/validation split. When the learning rate had reached its minimum and accuracy no longer improved on the validation data, training was stopped and the best performing model was saved. Following this protocol, each of the 12 orders were trained in 5 separate runs containing 11 training tasks each, thus producing a total of 660 recognition models. After training, each model was tested on all available test images for the relevant run.Data analysisThe relative representation of species within different taxa were generated using the number of species present in the GBIF data for Norway within each taxon and the number of accepted species within that taxon present in the Norwegian Species Nomenclature Database38, in line with Troudet et al.5: (R_x = n_x – (n frac{s_x}{s})) where (R_x) is the relative representation for taxon (x), (n_x) is the number of observations for taxon (x), (n) is the total number of observations for all taxa, (s_x) is the number of species within taxon (x), and (s) is the total number of species within all taxa.As a measure of model performance, we use the F(_1) score, the harmonic mean of the model’s precision and recall, given by$$begin{aligned} F_1 = frac{tp}{tp + frac{1}{2}(fp + fn)} end{aligned}$$where (tp), (fp) and (fn) stand for true positives, false positives and false negatives, respectively. The F(_1) score is a commonly used metric for model evaluation, as it is less susceptible to data imbalance than model accuracy28.The value of information (VoI) can be generically defined as “the increase in expected value that arises from making the best choice with the benefit of a piece of information compared to the best choice without the benefit of that same information”32. In the current context, we define the VoI as the expected increase in model performance (F(_1) score) when adding one observation with at least one image. To estimate this, for every order included in the experiment, the increase in average F(_1) score over increasing training task sizes were fitted using the Von Bertalanffy Growth Function, given by$$begin{aligned} L = L_infty (1 – e^{-k(t-t_0)}) end{aligned}$$where (L) is the average F(_1) score, (L_infty) is the asymptotic maximum F(_1) score, (k) is the growth rate, (t) is the number of observations per species, and (t_0) is a hypothetical number of observations at which the F(_1) score is 0. The Von Bertalanffy curve was chosen as it contains a limited number of parameters which are intuitive to interpret, and fits the growth of model performance well.The estimated increase in performance at any given point is then given by the slope of this function, i.e. the result of the differentiation of the Von Bertalanffy Growth Curve, given41 by$$begin{aligned} frac{dL}{dt} = bke^{-kt} end{aligned}$$where$$begin{aligned} b = L_infty e^{kt_0} end{aligned}$$Using this derivative function, we can estimate the expected performance increase stemming from one additional observation with images for each of the species within the order. Filling in the average number of citizen science observations with images per Norwegian species in that order for t, and dividing the result by the total number of Norwegian species within the order, provides the VoI of one additional observation with images for that order, expressed as an average expected F(_1) increase. More

Reference data collectionIn February 2019, we involved forest experts from different regions around the world and organized a workshop to (1) discuss the variety of forest management practices that take place in various parts of the world; (2) explore what types of forest management information could be collected by visual interpretation of very high-resolution images from Google Maps and Microsoft Bing Maps, in combination with Sentinel time series and Normalized Difference Vegetation Index (NDVI) profiles derived from Google Earth Engine (GEE); (3) generalize and harmonize the definitions at global scale; (4) finalize the Geo-Wiki interface for the crowdsourcing campaigns; and (5) build a data set of control points (or the expert data set), which we used later to monitor the quality of the crowdsourced contributions by the participants. Based on the results of this analysis, we launched the crowdsourcing campaigns by involving a broader group of participants, which included people recruited from remote sensing, geography and forest research institutes and universities. After the crowdsourcing campaigns, we collected additional data with the help of experts. Hence, the final reference data consists of two parts: (1) a randomly stratified sample collected by crowdsourcing (49,982 locations); (2) a targeted sample collected by experts (176,340 locations, at those locations where the information collected from the crowdsourcing campaign was not large enough to ensure a robust classification).DefinitionsTable 1 contains the initial classification used for visual interpretation of the reference samples and the aggregated classes presented in the final reference data set. For the Geo-Wiki campaigns, we attempted to collect information (1) related to forest management practices and (2) recognizable from very high-resolution satellite imagery or time series of vegetation indices. The final reference data set and the final map contain an aggregation of classes, i.e., only those that were reliably distinguishable from visual interpretation of satellite imagery.Table 1 Forest management classes and definitions.Full size tableSampling design for the crowdsourcing campaignsInitially, we generated a random stratified sample of 110,000 sites globally. The total number of sample sites was chosen based on experiences from past Geo-Wiki campaigns12, a practical estimation of the potential number of volunteer participants that we could engage in the campaign, and the expected spatial variation in forest management. We used two spatial data sets for the stratification of the sample: World Wildlife Fund (WWF) Terrestrial Ecoregions13 and Global Forest Change14. The samples were stratified into three biomes, based on WWF Terrestrial Ecoregions (Fig. 2): boreal (25 000 sample sites), temperate (35,000 sample sites) and tropical (50,000 sample sites). Within each biome, we used Hansen’s14 Global Forest Change maps to derive areas with “forest remaining forest” 2000–2015, “forest loss or gain”, and “permanent non-forest” areas.Fig. 2Biomes for sampling stratification (1 – boreal, 2 – temperate, 3 – sub-tropical and tropical).Full size imageThe sample size was determined from previous experiences, taking into account the expected spatial variation in forest management within each biome. Tropical forests had the largest sample size because of increasing commodity-driven deforestation15, the wide spatial extent of plantations, and slash and burn agriculture. Temperate forests had a larger sample compared to boreal forests due to their higher fragmentation. Each sample site was classified by at least three different participants, thus accounting for human error and varying expertise16,17,18. At a later stage, following a preliminary analysis of the data collected, we increased the number of sample sites to meet certain accuracy thresholds for every mapped class (aiming to exceed 75% accuracy).The Geo‐Wiki applicationGeo‐Wiki.org is an online application for crowdsourcing and expert visual interpretation of satellite imagery, e.g., to classify land cover and land use. This application has been used in several data collection campaigns over the last decade16,19,20,21,22,23. Here, we implemented a new custom branch of Geo‐Wiki (‘Human impact on Forest’), which is devoted to the collection of forest management data (Fig. 3). Various map overlays (including satellite images from Google Maps, Microsoft Bing Maps and Sentinel 2), campaign statistics and tools to aid interpretation, such as time series profiles of NDVI, were provided as part of this Geo‐Wiki branch, giving users a range of options and choices to facilitate image classification and general data collection. Google Maps and Microsoft Bing Maps include mosaics of very high-resolution satellite and aerial imagery from different time periods and multiple image providers, including the Landsat satellites operated by NASA and USGS as base imagery to commercial image providers such as Digital Globe. More information on the spatial and temporal distribution of very high-resolution satellite imagery can be found in Lesiv et al.24. This collection of images was supplied as guidance for visual interpretation16,20. Participants could analyze time series profiles of NDVI from Landsat, Sentinel 2 and MODIS images, which were derived from Google Earth Engine (GEE). More information on tools can be found in Supplementary file 1.Fig. 3Screenshot of the Geo‐Wiki interface showing a very high-resolution image from Google Maps and a sample site as a 100 mx100 m blue square, which the participants classified based on the forest management classes on the right.Full size imageThe blue box in Fig. 3 corresponds to 100 m × 100 m pixels aligned with the Sentinel grid in UTM projection. It is the same geometry required for the classification workflow that is used to produce the Copernicus Land Cover product for 201511.Before starting the campaign, the participants were shown a series of slides designed to help them gain familiarity with the interface and to train them in how to visually determine and select the most appropriate type of land use and forest management classes at each given location, thereby increasing both consistency and accuracy of the labelling tasks among experts. Once completed, the participants were shown random locations (from the random stratified sample) on the Geo‐Wiki interface and were then asked to select one of the forest management classes outlined in the Definition section (see Table 1 above).Alternatively, if there was either insufficient quality in the available imagery, or if a participant was unable to determine the forest management type, they could skip such a site (Fig. 3). If a participant skipped a sample site because it was too difficult, other participants would then receive this sample site for classification, whereas in the case of the absence of high-resolution satellite imagery, i.e., Google Maps and Microsoft Bing Maps, this sample site was then removed from the pool of available sample sites. The skipped locations were less than 1% of the total amount of locations assigned for labeling. Table 2 shows the distribution of the skipped locations by countries, based on the subset of the crowdsourced data where all the participants agreed.Table 2 Distribution of the skipped locations by countries.Full size tableQuality assurance and data aggregation of the crowdsourced dataBased on the experience gained from previous crowdsourcing campaigns12,19, we invested in the training of the participants (130 persons in total) and overall quality assurance. Specifically, we provided initial guidelines for the participants in the form of a video and a presentation that were shown before the participants could start classifying in the forest management branch (Supplementary file 1). Additionally, the participants were asked to classify 20 training samples before contributing to the campaign. For each of these training samples, they received text‐based feedback regarding how each location should be classified. Summary information about the participants who filled in the survey at the end of the campaign (i.e., gender, age, level of education, and their country of residence) is provided in the Supplementary file 2. We would like to note that 130 participants is a high number, especially taking the complexity of the task into consideration.Furthermore, during the campaign, sample sites that were part of the “control” data set were randomly shown to the participants. The participants received text-based feedback regarding whether the classification had been made correctly or not, with additional information and guidance. By providing immediate feedback, our intention was that participants would learn from their mistakes, increasing the quality and classification accuracy over time. If the text‐based feedback was not sufficient to provide an understanding of the correct classification, the participants were able to submit a request (“Ask the expert”) for a more detailed explanation by email.The control set was independent of the main sample, and it was created using the same random stratified sampling procedure within each biome and the stratification by Global Forest Change maps14 (see “Sample design” section). To determine the size of the control sample, we considered two aspects: (a) the maximum number of sample sites that one person could classify during the entire campaign; (b) the frequency at which control sites would appear among the task sites (defined at 15%, which is a compromise between the classification of as many unknown locations as possible and a sufficient level of quality control, based on previous experience). Our control sample consisted of 5,000 sites. Each control sample site was classified twice by two different experts. When the two experts agreed, these sample sites were added to the final control sample. Where disagreement occurred (in 25% of cases), these sample sites were checked again by the experts and revised accordingly. During the campaign, participants had the option to disagree with the classification of the control site and submit a request with their opinion and arguments. They received an additional quality score in the situation when they were correct, but the experts were not. This procedure also ensured an increase in the quality of the control data set.To incentivize participation and high-quality classifications, we offered prizes as part of the campaign design. The ranking system for the prize competition considered both the quality of the classifications and the number of classifications provided by a participant. The quality measure was based on the control sample discussed above. The participants randomly received a control point, which was classified in advance by the experts. For every control point, a participant could receive a maximum of +30 points (fully correct classification) to a minimum of −30 points (incorrect classification). In the case where the answer was partly correct (e.g., the participant correctly classified that the forest is managed, but misclassified the regeneration type), they received points ranging from 5 to 25.The relative quality score for each participant was then calculated as the total sum of gained points divided by the maximum sum of points that this participant could have earned. For any subsequent data analysis, we excluded classifications from those participants whose relative quality score was less than 70%. This threshold corresponds to an average score of 10 points at each location (out of a maximum of 30 points), i.e., where participants were good at defining the aggregated forest management type but may have been less good at providing the more detailed classification.Unfortunately, we observed some imbalance in the proportion of participants coming from different countries, e.g. there were not so many participants from the tropics. This could have resulted in interpretation errors, even when all the participants agreed on a classification. To address this, we did an additional quality check. We selected only those sample sites where all the participants agreed and then randomly checked 100 sample sites from each class. Table 3 summarizes the results of this check and explains the selection of the final classes presented in Table 1.Table 3 Qualitative analysis of the reference sample sites with full agreement.Full size tableAs a result of the actions outlined in Table 3, we compiled the final reference data set, which consisted of 49,982 consistent sample sites.Additional expert data collectionWe used the reference data set to produce a test map of forest management (the classification algorithm used is described in the next section). By checking visually and comparing against the control data set, we found that the map was of insufficient quality for many locations, especially in the case of heterogeneous landscapes. While several reasons for such an unsatisfactory result are possible, the experts agreed that a larger sample size would likely increase the accuracy of the final map, especially in areas of high heterogeneity and for forest management classes that only cover a small spatial extent. To increase the amount of high-quality training data and hence to improve the map, we collected additional data using a targeted approach. In practice, the map was uploaded to Geo-Wiki, and using the embedded drawing tools, the experts randomly checked locations on the map, focusing on their region of expertise and added classified polygons in locations where the forest management was misclassified. To limit model overfitting and oversampling of certain classes, the experts also added points for correctly mapped classes to keep the density of the points the same. This process involved a few iterations of collecting additional points and training the classification algorithm until the map accuracy reached 75%. In total, we collected an additional 176,340 training points. With the 49,982 consistent training points from the Geo-Wiki campaigns, this resulted in 226,322 (Fig. 4). This two-pronged approach would not have been possible without the exhaustive knowledge obtained from running the initial Geo-Wiki campaigns, including numerous questions raised by the campaign participants. Figure 4 also highlights in yellow the areas of very high sampling density, I.e., those collected by the experts. The sampling intensity of these areas is much higher in comparison with the randomly distributed crowdsourced locations, and these are mainly areas with very mixed forest classes or small patches, in most cases, including plantations.Fig. 4Distribution of reference locations.Full size imageClassification algorithmTo produce the forest management map for the year 2015, we applied a workflow that was developed as part of the production of the Copernicus Global Land Services land cover at 100 m resolution (CGLS-LC100) collection 2 product11. A brief description of the workflow (Fig. 5), focusing on the implemented changes, is given below. A more thorough explanation, including detailed technical descriptions of the algorithms, the ancillary data used, and the intermediate products generated, can be found in the Algorithm Theoretical Basis Document (ATBD) of the CGLS-LC100 collection 2 product25.Fig. 5Workflow overview for the generation of the Copernicus Global Land Cover Layers. Adapted from the Algorithm Theoretical Basis Document25.Full size imageThe CGLS-LC100 collection 2 processing workflow can be applied to any satellite data, as it is unspecific to different sensors or resolutions. While the CGLS-LC100 Collection 2 product is based on PROBA-V sensor data, the workflow has already been tested with Sentinel 2 and Landsat data, thereby using it for regional/continental land cover (LC) mapping applications11,26. For generating the forest management layer, the main Earth Observation (EO) input was the PROBA-V UTM Analysis Ready Data (ARD) archive based on the complete PROBA-V L1C archive from 2014 to 2016. The ARD pre-processing included geometric transformation into a UTM coordinate system, which reduced distortions in high northern latitudes, as well as improved atmospheric correction, which converted the Top-of-Atmosphere reflectance to surface reflectance (Top-of-Canopy). In a further processing step, gaps in the 5-daily PROBA-V UTM multi-spectral image data with a Ground Sampling Distance (GSD) of ~0.001 degrees (~100 m) were filled using the PROBA-V UTM daily multi-spectral image data with a GSD of ~0.003 degrees (~300 m). This data fusion is based on a Kalman filtering approach, as in Sedano et al.27, but was further adapted to heterogonous surfaces25. Outputs from the EO pre-processing were temporally cleaned by using the internal quality flags of the PROBA-V UTM L3 data, a temporal cloud and outlier filter built on a Fourier transformation. This was done to produce consistent and dense 5-daily image stacks for all global land masses at 100 m resolution and a quality indicator, called the Data Density Indicator (DDI), used in the supervised learning process of the algorithm.Since the total time series stack for the epoch 2015 (a three-year period including the reference year 2015 +/− 1 year) would be composed of too many proxies for supervised learning, the time and spectral dimension of the data stack had to be condensed. The spectral domain was condensed by using Vegetation Indices (VIs) instead of the original reflectance values. Overall, ten VIs based on the four PROBA-V reflectance bands were generated, which included: Normalized Difference Vegetation Index (NDVI); Enhanced Vegetation Index (EVI); Structure Intensive Pigment Index (SIPI); Normalized Difference Moisture Index (NDMI); Near-Infrared reflectance of vegetation (NIRv); Angle at NIR; HUE and VALUE of the Hue Saturation Value (HSV) color system transformation. The temporal domain of the time series VI stacks was then condensed by extracting metrics, which are used as general descriptors to enable distinguishing between the different LC classes. Overall, we extracted 266 temporal, descriptive, and textual metrics from the VI times series stacks. The temporal descriptors were derived through a harmonic model, fitted through the time series of each of the VIs based on a Fourier transformation28,29. In addition to the seven parameters of the harmonic model that describe the overall level and seasonality of the VI time series, 11 descriptive statistics (mean, standard deviation, minimum, maximum, sum, median, 10th percentile, 90th percentile, 10th – 90th percentile range, time step of the first minimum appearance, and time step of the first maximum appearance) and one textural metric (median variation of the center pixel to median of the neighbours) were generated for each VI. Additionally, the elevation, slope, aspect, and purity derived at 100 m from a Digital Elevation Model (DEM) were added. Overall, 270 metrics were extracted from the PROBA-V UTM 2015 epoch.The main difference to the original CGLS-LC100 collection 2 algorithms is the use of forest management training data instead of the global LC reference data set, as well as only using the discrete classification branch of the algorithm. The dedicated regressor branch of the CGLS-LC100 collection 2 algorithm, i.e., outputting cover fraction maps for all LC classes, was not needed for generating the forest management layer.In order to adapt the classification algorithm to sub-continental and continental patterns, the classification of the data was carried out per biome cluster, with the 73 biome clusters defined by the combination of several global ecological layers, which include the ecoregions 2017 dataset30, the Geiger-Koeppen dataset31, the global FAO eco-regions dataset32, a global tree-line layer33, the Sentinel-2 tiling grid and the PROBA-V imaging extent;30,31 this, effectively, resulted in the creation of 73 classification models, each with its non-overlapping geographic extent and its own training dataset. Next, in preparation for the classification procedure, the metrics of all training points were analyzed for outliers, as well as screened via an all-relevant feature selection approach for the best metric combinations (i.e., best band selection) for each biome cluster in order to reduce redundancy between parameters used in the classification. The best metrics are defined as those that have the highest separability compared to other metrics. For each metric, the separability is calculated by comparing the metric values of one class to the metric values of another class; more details can be found in the ATBD25. The optimized training data set, together with the quality indicator of the input data (DDI data set) as a weight factor, were used in the training of the Random Forest classifier. Moreover, a 5-fold cross-validation was used to optimize the classifier parameters for each generated model (one per biome).Finally, the Random Forest classification was used to produce a hard classification, showing the discrete class for each pixel, as well as the predicted class probability. In the last step, the discrete classification results (now called the forest management map) are modified by the CGLS-LC100 collection 2 tree cover fraction layer29. Therefore, the tree cover fraction layer, showing the relative distribution of trees within one pixel, was used to remove areas with less than 10% tree cover fraction in the forest management layer, following the FAO definition of forest. Figure 6 shows the class probability layer that illustrates the model behavior, highlighting the areas of class confusion. This layer shows that there is high confusion between forest management classes in heterogeneous landscapes, e.g., in Europe and the Tropics while homogenous landscapes, such as Boreal forests, are mapped with high confidence. It is important to note that a low probability does not mean that the classification is wrong.Fig. 6The predicted class probability by the Random Forest classification.Full size image More

Data collectionThe acoustic recordings were collected during 2017 off the Changhua coast (24° 4.283 N/120° 19.102 E) (Fig. 5) by deploying a passive acoustic monitoring (PAM) device from Wildlife Acoustics, which was an SM3M recorder moored at a depth of 18–20 m. The hydrophone recorded continuously with a sampling frequency of 48 kHz and a sensitivity of −164.2 dB re:1 v/µPa. The acoustic files were recorded in the.WAV format with a duration of 60 minutes. The hydrophone setup prior to deployment is shown in Fig. 6. Table 2 contains the details for the monitoring period with the corresponding season and the number of hours of recordings each season used in this study. Studies have shown that the presence of seasonal chorusing at this monitoring site in the frequency range of 500–2500 Hz is caused by two types of chorusing15,38, with chorusing starting in early spring, reaching a peak in summer, and starting to decline late autumn, and silencing in winter6. Previous studies6,15,38 at this monitoring site have derived the details of two types of fish calls responsible for chorusing (Type 1 and Type 2); Supplementary Fig. 1 shows the spectrogram, waveform, and power spectrum density of the individual calls. Supplementary Table 1 tabulated are the acoustic features of the two call types. The monitoring region, Changhua, lies in the Eastern Taiwan Strait (ETS). The ETS is ~350 km in length and ~180 km wide64. The ETS experiences diverse oceanographic and climatic variations influenced by monsoons in summer and winter65 and extreme events caused by tropical storms, typhoons in summer, and wind/cold bursts occurring in winter66,67,68.Fig. 5: Study area located off the Taiwan Strait.Map of the Changhua coast located in Taiwan Strait, Taiwan depicting the deployed passive acoustic monitoring recorder at site A1. The map was produced in Matlab 9.11 (The Mathworks, Natick, MA; http://www.mathworks.com/) using mapping toolbox function geobasemap(). Full global basemap composed of high-resolution satellite imagery hosted by Esri (https://www.esri.com/).Full size imageFig. 6: Setup of the SM3M submersible recorder.SM3M recorder fastened to the steel frame (length and breadth = 1.22 m, height = 0.52 m) with plastic cable zip ties prior to deployment.Full size imageTable 2 Passive acoustic monitoring device specifications and monitoring duration during different seasons.Full size tableAcoustic data analysisThe acoustic data were analyzed using the PAMGuide toolbox in Matlab60. The seasonal spectrograms were computed with an FFT size of 1024 points and a 1 s time segment averaged to a 60 s resolution. The sound pressure levels (SPL) were computed in the frequency band of 500–3500 Hz and programmed to provide a single value every hour, thus resulting in 984, 1344, and 1440 data points in spring, summer, and winter, respectively (Supplementary Data 1).Determining the regularity and complexity with the complexity-entropy planeThe complexity-entropy plane was utilized in this study to quantify the structural statistical complexity and the regularity in the hourly acoustical and seasonal SPL time series data. The C-H plane is a 2D plane representation of the permutation entropy on the horizontal axis that quantifies the regularity, and the vertical axis is represented by the statistical complexity quantifying the correlation structure in the temporal series.For a given time series ({{x(t)}}_{t=1}^{N}), the N’ ≡ N − (m − 1) the values of the vectors for the length m  > 1 are ranked as$${X}_{s}=left({x}_{s-(m-1)},{x}_{s-(m-2)},ldots ,{x}_{s}right),s=1,ldots ,,{N}^{{prime} }$$
(1)
Within each vector, the values are reordered in the ascending order of their amplitude, yielding the set of ordering symbols (patterns) ({r}_{0},{r}_{1},ldots ,{r}_{m-1}) such that$${x}_{s-{r}_{0}}le {x}_{s-{r}_{1}}le ..,..le {x}_{s-{r}_{(m-1)}}$$
(2)
This symbolization scheme was introduced by Bandt and Pompe69. The scheme performs the local ordering of a time series to construct a probability mass function (PMF) of the ordinal patterns of the vector. The corresponding vectors (pi ={r}_{0},{r}_{1},ldots ,{r}_{(m-1)}) may presume any of the m! possible permutations of the set ({{{{{mathrm{0,1}}}}},ldots ,m-1}) and symbolically represent the original vector. For instance, for a given time series {9, 4, 5, 6, 1,…} with length m = 3, provides 3! different order patterns with six mutually exclusive permutation symbols are considered. The first three-dimensional vector is (9, 4, 5), following the Eq. (1), this vector corresponds to ((,{x}_{s-2},{x}_{s-1},{x}_{s})). According to Eq. (2), it yields ({x}_{s-1}le {x}_{s}le {x}_{s-2}). Then, the ordinal pattern satisfying the Eq. (2) will be (1, 0, 2). The second 3-dimensional vector is (4, 5, 6), and (2, 1, 0) will be its associated permutation, and so on.The permutation entropy (H) of order m ≥ 2 is defined as the Shannon entropy of the Brandt-Pompe probability distribution p(π)69$$Hleft(mright)=,-{mathop{sum}limits _{{pi }}}pleft(pi right){{{log }}}_{2}p(pi )$$
(3)
where ({pi }) represents the summation over all possible m! permutations of order m, (p(pi )) is the relative frequency of each permutation (pi), and the binary logarithm (base of 2) is evaluated to quantify the entropy in bits. H(m) attains the maximum ({{log }}m!) for (p(pi )=1/m!). Then the normalized Shannon entropy is given by$$0le H(m)/{{{{{rm{ln}}}}}},m!le 1$$
(4)
where the lower bound H = 0 corresponds to more predictable signals with fewer fluctuations, an either strictly increasing or decreasing series (with a single permutation), and the upper bound H = 1 corresponds to an unpredictable random series for which all the m! possible permutations are equiprobable. Thus, H quantifies the degree of disorder inherent in the time series. The choice of the pattern length m is essential for calculating the appropriate probability distribution, particularly for m, which determines the number of accessible states given by m!70,71. As a rule of thumb, the length T of the time series must satisfy the condition T (gg) m! in order to obtain reliable statistics, and for practical purposes, Bandt and Pompe suggested choosing m = 3,…,7 69.The statistical complexity measure is defined with the product form as a function of the Bandt and Pompe probability distribution P associated with the time series. (Cleft[Pright]) is represented as33$$Cleft[Pright]=frac{J[P,U]}{{J}_{{max }}}{H}_{s}[P]$$
(5)
where ({H}_{s}left[Pright]=Hleft[Pright]/{{log }}m!) is the normalized permutation entropy. (J[P,U]) is the Jensen divergence$$Jleft[P,Uright]=left{Hleft[frac{P+U}{2}right]-frac{H[P]}{2}-frac{H[U]}{2}right}$$
(6)
which quantifies the difference between the uniform distributions U and P, and ({J}_{{max }})is the maximum possible value of (Jleft[P,Uright]) that is obtained from one of the components of P = 1, with all the other components being zero:$$Jleft[P,Uright]=-frac{1}{2}left[frac{m!+1}{m!}{{log }}left(m!+1right)-2{{log }}left(2m!right)+{{log }}(m!)right]$$
(7)
For each value of the normalized permutation entropy (0le {H}_{s}[P]le 1) there is a corresponding range of possible statistical complexity (Cleft[Pright]) values. Thus, the upper (({C}_{{max }})) and lower ((C_{{min }})) complexity bounds in the C-H plane are formed. The periodic sequences such as sine and series with increasing and decreasing (with ({H}_{s}[P]=0)) and completely random series such as white noise (for which (Jleft[P,Uright]=0) and ({H}_{s}[P]=1)) will have zero complexity. Furthermore, for each given value of the (0le {H}_{s}[P]le 1), there exists a range of possible values of the statistical complexity, ({C}_{{min }}le C[P]le {C}_{{max }}). The procedure for evaluating the complexity bounds ({C}_{{min }}) and ({C}_{{max }}) is given in Martin et al.72. We utilized the R package ‘statcomp’73 to evaluate the statistical complexity (C) and the permutation entropy (H) using the command global-complexity() for the order m = 6, and the command limit_curves(m, fun = ‘min/max’) was utilized to evaluate the complexity boundaries ({C}_{{min }}) and ({C}_{{max }}). In this study, we constructed two C-H planes: (1) C and H was computed for each hourly acoustic file during different seasons. The resulting lengths of C and H during spring, summer, and autumn-winter are similar to the number of hours in the particular season (Table 2). (2) C and H was computed every 4–5 days for the seasonal SPL. The resulting length of C and H obtained during spring was 9 points (each value of C and H for every 109 h), and in summer and autumn-winter was 12 points (each value of C and H for every 112 and 120 h).Determining predictability and dynamics (linear/nonlinear) using EDMIn this study, we utilized EDM to quantify the predictability (forecasting) and distinguish between the linear stochastic and nonlinear dynamics in the seasonal soundscape SPL. EDM involves phase-space reconstruction via delay coordinate embeddings to make forecasts and to determine the ‘predictability portrait’ of time series data40. The first step in EDM is to determine the optimal embedding dimension (E), and this is obtained using a method based on simplex projection41. The simplex projection is carried out by dividing the dataset into two equal parts, of which the first part is called the library and the other part is called the target. The library set is used to build a series of non-parametric models (known as predictors) for the one step ahead predictions for the E varying between 1 and 10. Then the model’s accuracies are determined when the model is applied to the target dataset and the prediction skill (⍴) for the actual and predicted datasets is measured. The embedding dimension with the highest predictive skill is the optimal E.For the appropriate optimal E chosen, the predictability profile of the time series data is obtained by evaluating ⍴ at Tp = 1, 2, 3, … steps ahead. The flat prediction profile of the ⍴–Tp curve indicates that the time series is purely random (low ⍴) or regularly oscillating (high ⍴). In contrast, a decreasing ⍴ as Tp increases may reject the possibility of an underlying uncorrelated stochastic process and indicate the presence of low-dimensional deterministic dynamics. However, the concern with the predictability profile is that it may exhibit predictability even if time series are purely stochastic (such as autocorrelated red noise). Hence, a nonlinear test can be performed by using S-maps (sequential locally weighted global linear maps) to distinguish between linear stochastic and nonlinear dynamics in the time series dataset by fitting a local linear map. S-maps similar to simplex projects provide the forecasts in phase-space by quantifying the degree to which points are weighted when fitting the local linear map, which is given by the nonlinear localization parameter θ. When θ = 0, the entire library set will exhibit equal weights regardless of the target prediction, which mathematically resembles the model of a linear autoregressive process. In contrast, if θ  > 0, the forecasts of the library provided by the S-map depend on the local state of the target prediction, thus producing large weights, and the unique local fittings can vary in phase-space to incorporate nonlinear behavior. Consequently, if the (⍴–θ) dynamics profile shows the highest ⍴ at θ = 0 that is reduced as θ increases, it represents linear stochastic dynamics. If the ⍴ achieves the highest value at θ  > 0, then the dynamics are represented by nonlinear dynamics.In this study, the R package “rEDM”74 was used to evaluate the optimal E, prediction profile (⍴–Tp), and dynamics profile (⍴–θ) for the seasonal SPL dataset. While evaluating these entities, the data points are equally into two parts, and sequentially the first half is chosen as the library set and the other as the target set. The length of the library and the target set for spring, summer, and autumn-winter are 492, 672, and 720. The command EmbedDimension() was used to determine the forecast skill for the E ranging from 1 to 10 and the optimal E with the highest forecast skill (Supplementary Fig. 2) was chosen. In this study, we found that for all seasons, the optimal E was 2. The (⍴–Tp) was evaluated for Tp varying between 1 and 100 using the command PredictInterval() and the (⍴–θ) was evaluated using the command PredictNonlinear() for θ = 0, 0.0001, 0.0003, 0.001, 0.003, 0.01, 0.03, 0.1, 0.3, 0.5,0.75, 1.0, 1.5, 2, and 3 to 20.StatisticsThe nonparametric Kruskal–Wallis test, followed by post hoc Bonferroni’s multiple comparisons, was used to test differences in the seasonal H and C that were obtained directly from the hourly acoustic data during chorusing hours, as well as the H and C obtained for the seasonal SPL and the seasonal forecast skill. The statistical calculations were performed using the R package “agricolae” 75. More

Ellison, W. T., Southall, B. L., Clark, C. W. & Frankel, A. S. A new context-based approach to assess marine mammal behavioral responses to anthropogenic sounds. Conserv. Biol. 26, 21–28 (2012).CAS 
PubMed 
Article 

Google Scholar 
Williams, R., Clark, C. W., Ponirakis, D. & Ashe, E. Acoustic quality of critical habitats for three threatened whale populations. Anim. Conserv. 17, 174–185 (2014).Article 

Google Scholar 
Halliday, W. D., Pine, M. K. & Insley, S. J. Underwater noise and arctic marine mammals: Review and policy recommendations. Environ. Rev. 28, 438–448 (2020).Article 

Google Scholar 
Kvadsheim, P. H. et al. Impact of Anthropogenic Noise on the Marine Environment: Status of Knowledge and Management (Springer, 2020).
Google Scholar 
Weilgart, L. S. & Whitehead, H. Distinctive vocalizations from mature male sperm whales (Physeter macrocephalus). Can. J. Zool. 66, 1931–1937 (1988).Article 

Google Scholar 
Simon, M., Stafford, K. M., Beedholm, K., Lee, C. M. & Madsen, P. T. Singing behavior of fin whales in the Davis Strait with implications for mating, migration and foraging. J. Acoust. Soc. Am. 128, 3200 (2010).ADS 
PubMed 
Article 

Google Scholar 
Alves, D., Amorim, M. C. P. & Fonseca, P. J. Assessing acoustic communication active space in the Lusitanian toadfish. J. Exp. Biol. 219, 1122–1129 (2016).PubMed 

Google Scholar 
Linnenschmidt, M., Teilmann, J., Akamatsu, T., Dietz, R. & Miller, L. A. Biosonar, dive, and foraging activity of satellite tracked harbor porpoises (Phocoena phocoena). Mar. Mamm. Sci. 29, E77–E97 (2013).Article 

Google Scholar 
Giorli, G. & Goetz, K. T. Foraging activity of sperm whales (Physeter macrocephalus) off the east coast of New Zealand. Sci. Rep. 9, 12182 (2019).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Baumgartner, M. F. & Fratantoni, D. M. Diel periodicity in both sei whale vocalization rates and the vertical migration of their copepod prey observed from ocean gliders. Limnol. Oceanogr. 53, 2197–2209 (2008).ADS 
Article 

Google Scholar 
Urazghildiiev, I. R. & Van Parijs, S. M. Automatic grunt detector and recognizer for Atlantic cod (Gadus morhua). J. Acoust. Soc. Am. 139, 2532–2540 (2016).ADS 
PubMed 
Article 

Google Scholar 
Ladich, F. Ecology of sound communication in fishes. Fish Fish. 20, 552–563 (2019).Article 

Google Scholar 
Radford, C. A., Stanley, J. A., Simpson, S. D. & Jeffs, A. G. Juvenile coral reef fish use sound to locate habitats. Coral Reefs 30, 295–305 (2011).ADS 
Article 

Google Scholar 
Pierpoint, C. Harbour porpoise (Phocoena phocoena) foraging strategy at a high energy, near-shore site in south-west Wales, UK. J. Mar. Biol. Assoc. UK 88, 1167–1173 (2008).Article 

Google Scholar 
Pijanowski, B. C. et al. Soundscape ecology: The science of sound in the landscape. Bioscience 61, 203–216 (2011).Article 

Google Scholar 
Stanley, J. A., Radford, C. A. & Jeffs, A. G. Location, location, location: Finding a suitable home among the noise. Proc. R. Soc. B Biol. Sci. 279, 3622–3631 (2012).Article 

Google Scholar 
Buscaino, G. et al. Temporal patterns in the soundscape of the shallow waters of a Mediterranean marine protected area. Sci. Rep. 6, 1–13 (2016).Article 
CAS 

Google Scholar 
Gasc, A., Francomano, D., Dunning, J. B. & Pijanowski, B. C. Future directions for soundscape ecology: The importance of ornithological contributions. Auk 134, 215–228 (2017).Article 

Google Scholar 
Putland, R. L., Constantine, R. & Radford, C. A. Exploring spatial and temporal trends in the soundscape of an ecologically significant embayment. Sci. Rep. 7, 5713 (2017).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Pieretti, N., LoMartire, M., Farina, A. & Danovaro, R. Marine soundscape as an additional biodiversity monitoring tool: A case study from the Adriatic Sea (Mediterranean Sea). Ecol. Indic. 83, 13–20 (2017).Article 

Google Scholar 
Gillespie, D., Palmer, L., Macaulay, J., Sparling, C. & Hastie, G. Passive acoustic methods for tracking the 3D movements of small cetaceans around marine structures. PLoS ONE 15, e0229058 (2020).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Van Parijs, S. et al. Management and research applications of real-time and archival passive acoustic sensors over varying temporal and spatial scales. Mar. Ecol. Prog. Ser. 395, 21–36 (2009).ADS 
Article 

Google Scholar 
Warren, V. E., McPherson, C., Giorli, G., Goetz, K. T. & Radford, C. A. Marine soundscape variation reveals insights into baleen whales and their environment: a case study in central New Zealand. R. Soc. Open Sci. 8, 201503 (2021).ADS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Ahonen, H. et al. The underwater soundscape in western Fram Strait: Breeding ground of Spitsbergen’s endangered bowhead whales. Mar. Pollut. Bull. 123, 97–112 (2017).CAS 
PubMed 
Article 

Google Scholar 
Hildebrand, J. A. Anthropogenic and natural sources of ambient noise in the ocean. Mar. Ecol. Prog. Ser. 395, 5–20 (2009).ADS 
Article 

Google Scholar 
Ross, D. Ship sources of ambient noise. IEEE J. Ocean. Eng. 30, 257–261 (2005).ADS 
Article 

Google Scholar 
Popper, A. N. & Hawkins, A. The Effects of Noise on Aquatic Life Vol. 730 (Springer, 2012).Book 

Google Scholar 
Hubert, J. et al. Effects of broadband sound exposure on the interaction between foraging crab and shrimp: A field study. Environ. Pollut. 243, 1923–1929 (2018).CAS 
PubMed 
Article 

Google Scholar 
Weilgart, L. The Impact of Ocean Noise Pollution on Fish and Invertebrates (Springer, 2018).
Google Scholar 
Kvadsheim, H., Sivle, L. D., Hansen, R. R. & Karlsen, H. E. Effekter av Menneskeskapt støy på Havmiljø Rapport til Miljødirektoratet om Kunnskapsstatus FFI-RAPPORT. (2017).Parks, S. E., Johnson, M., Nowacek, D. & Tyack, P. L. Individual right whales call louder in increased environmental noise. Biol. Lett. 7, 33–35 (2011).PubMed 
Article 

Google Scholar 
Meh, F. et al. Humpback whales Megaptera novaeangliae alter calling behavior in response to natural sounds and vessel noise. Mar. Ecol. Prog. Ser. 607, 251–268 (2018).Article 

Google Scholar 
Leroy, E. C., Royer, J.-Y., Bonnel, J. & Samaran, F. Long-term and seasonal changes of large whale call frequency in the Southern Indian ocean. J. Geophys. Res. Ocean. 123, 8568–8580 (2018).ADS 
Article 

Google Scholar 
Parks, S. E., Clark, C. W. & Tyack, P. L. Short- and long-term changes in right whale calling behavior: The potential effects of noise on acoustic communication. J. Acoust. Soc. Am. 122, 3725–3731 (2007).ADS 
PubMed 
Article 

Google Scholar 
Clark, C. et al. Acoustic masking in marine ecosystems: intuitions, analysis, and implication. Mar. Ecol. Prog. Ser. 395, 201–222 (2009).ADS 
Article 

Google Scholar 
PAME. Underwater Noise in the Arctic: A State of Knowledge Report (PAME, 2019).
Google Scholar 
Beszczynska-Möller, A., Woodgate, R., Lee, C., Melling, H. & Karcher, M. A synthesis of exchanges through the main oceanic gateways to the Arctic Ocean. Oceanography 24, 82–99 (2011).Article 

Google Scholar 
Ramm, T. Hungry During Migration? Humpback Whale Movement from the Barents Sea to a Feeding Stopover in Northern Norway Revealed by Photo-ID Analysis. (MSc thesis. UiT The Arctic University of Norway, 2020).Broms, F. et al. Recent research on the migratory destinations of humpback whales (Megaptera novaeangliae) from a mid-winter feeding stop-over area in Northern Norway. in Recent research on the migratory destinations of humpback whales (Megaptera novaeangliae) from a mid-winter feeding stop-over area in Northern Norway (ed. Wenzel, F. W.) (European Cetacean Society Special Publication Series, 2015).Aniceto, A. S. et al. Arctic marine data collection using oceanic gliders: Providing ecological context to cetacean vocalizations. Front. Mar. Sci. 7, 547 (2020).Article 

Google Scholar 
Jourdain, E. & Vongraven, D. Humpback whale (Megaptera novaeangliae) and killer whale (Orcinus orca) feeding aggregations for foraging on herring (Clupea harengus) in Northern Norway. Mamm. Biol. 86, 27–32 (2017).Article 

Google Scholar 
Christiansen, J. S., Mecklenburg, C. W. & Karamushko, O. V. Arctic marine fishes and their fisheries in light of global change. Glob. Chang. Biol. 20, 352–359 (2014).ADS 
PubMed 
Article 

Google Scholar 
Rødland, E. S. & Bjørge, A. Residency and abundance of sperm whales (Physeter macrocephalus) in the Bleik Canyon, Norway. Mar. Biol. Res. 11, 974–982 (2015).Article 

Google Scholar 
Nøttestad, L. et al. Prey selection of offshore killer whales Orcinus orca in the Northeast Atlantic in late summer: spatial associations with mackerel. Mar. Ecol. Prog. Ser. 499, 275–283 (2014).ADS 
Article 

Google Scholar 
Bjørge, A., Aarefjord, H., Kaarstad, S., Kleivane, L. & Øien, N. Harbour porpoise (Phocoena phocoena) in Norwegian waters (Springer, 1991).
Google Scholar 
Gjøseter, H. et al. Fisken og Havet. https://doi.org/10.1111/maec.12075 (2010).Article 

Google Scholar 
ICES. ICES Report on Ocean Climate 2009 No.304. (2010).ICES. Report of the Working Group on Widely Distributed Stocks (WGWIDE). (2010).Rey, F. Phytoplankton: The grass of the sea. In The Norwegian Sea Ecosystem (ed. Skjoldal, H. R.) 97–136 (Academic Press, 2004).
Google Scholar 
Huse, G. et al. Effects of interactions between fish populations on ecosystem dynamics in the Norwegian Sea : Results of the INFERNO project. Mar. Biol. Res. 8, 415–419 (2012).Article 

Google Scholar 
Godø, O. R., Johnsen, S. & Torkelsen, T. The LoVe ocean observatory is in operation. Mar. Technol. Soc. J. 48, 24–30 (2014).Article 

Google Scholar 
Cooke, J. G. Balaenoptera physalus. The IUCN Red List of Threatened Species: e.T2478A50349982 (2018).Leonard, D. & Øien, N. Estimated abundances of Cetacean species in the Northeast Atlantic from Norwegian Shipboard Surveys Conducted in 2014–2018. NAMMCO Sci. Publ. 11, 4694 (2020).
Google Scholar 
Øygard, S. H. Simulations of Acoustic Transmission Loss of Fin Whale Calls Reaching the LoVe Ocean Observatory. (MSc thesis. University of Bergen, 2018).Steiner, L. et al. A link between male sperm whales, Physeter macrocephalus, of the Azores and Norway. J. Mar. Biol. Assoc. UK 92, 1751–1756 (2012).Article 

Google Scholar 
Olafsen, T., Winther, U., Olsen, Y. & Skjermo, J. Verdiskaping Basert på Produktive hav i 2050 1–76 (Springer, 2012).
Google Scholar 
Wenz, G. M. Acoustic ambient noise in the ocean: Spectra and sources. J. Acoust. Soc. Am. 34, 1936–1956 (1962).ADS 
Article 

Google Scholar 
Klinck, H. et al. Seasonal presence of cetaceans and ambient noise levels in polar waters of the North Atlantic. J. Acoust. Soc. Am. 132, 176–181 (2012).Article 

Google Scholar 
Burnham, R. E., Duffus, D. A. & Mouy, X. The presence of large whale species in Clayoquot Sound and its offshore waters. Cont. Shelf Res. 177, 15–23 (2019).ADS 
Article 

Google Scholar 
Romagosa, M. et al. Baleen whale acoustic presence and behaviour at a Mid-Atlantic migratory habitat, the Azores Archipelago. Sci. Rep. 10, 61489 (2020).
Google Scholar 
Tervo, O. Acoustic Behaviour of Bowhead Whales Balaena mysticetus in Disko Bay, Western Greenland. PhD thesis. (2011).Magnúsdóttir, E. E. & Lim, R. Subarctic singers: Humpback whale (Megaptera novaeangliae) song structure and progression from an Icelandic feeding ground during winter. PLoS ONE 14, e0210057 (2019).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Samaran, F. et al. Seasonal and geographic variation of southern blue whale subspecies in the Indian Ocean. PLoS ONE 8, e70 (2013).Article 

Google Scholar 
Norris, T. F., Dunleavy, K. J., Yack, T. M. & Ferguson, E. L. Estimation of minke whale abundance from an acoustic line transect survey of the Mariana Islands. Mar. Mammal Sci. 33, 574 (2017).Article 

Google Scholar 
Marques, T. A. et al. Estimating animal population density using passive acoustics. Biol. Rev. Camb. Philos. Soc. 88, 287–309 (2013).PubMed 
Article 

Google Scholar 
Dunlop, R. A. The effects of vessel noise on the communication network of humpback whales. R. Soc. Open Sci. 6, 190967 (2019).ADS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Christensen, I., Haug, T. & Øien, N. A review of feeding and reproduction in large baleen whales (Mysticeti) and sperm whales Physeter macrocephalus in Norwegian and adjacent waters. ICES J. Mar. Sci. 49, 341–355 (1992).Article 

Google Scholar 
Aniceto, A. S. et al. Monitoring marine mammals using unmanned aerial vehicles: quantifying detection certainty. Ecosphere 9, e02122 (2018).Article 

Google Scholar 
Pedersen, G., Storheim, E., Sivle, L. D., Godø, O. R. & Ødegaard, L. A. Concurrent passive and active acoustic observations of high-latitude shallow foraging sperm whales (Physeter macrocephalus) and mesopelagic prey layer. J. Acoust. Soc. Am. 141, 1–10 (2017).Article 

Google Scholar 
Vogel, E. F. The influence of herring (Clupea harengus) biomass and distribution on killer whale (Orcinus orca) movements on the Norwegian shelf (UiT The Arctic University of Norway, 2020).
Google Scholar 
Williams, R. et al. Impacts of anthropogenic noise on marine life: Publication patterns, new discoveries, and future directions in research and management. Ocean Coast. Manag. 115, 17–24 (2015).Article 

Google Scholar 
Garibbo, S. et al. Low-frequency ocean acoustics: Measurements from the Lofoten-Vesterålen Ocean Observatory, Norway. (2020). https://doi.org/10.1121/2.0001324.Dekeling, R. P. A. et al. Monitoring Guidance for Underwater Noise in European Seas, Part I: Executive Summary (Springer, 2014).
Google Scholar 
Erbe, C. International regulation of underwater noise. Acoust. Aust. 41, 1–10 (2013).
Google Scholar 
Halliday, W. D., Insley, S. J., Hilliard, R. C., de Jong, T. & Pine, M. K. Potential impacts of shipping noise on marine mammals in the western Canadian Arctic. Mar. Pollut. Bull. 123, 73–82 (2017).CAS 
PubMed 
Article 

Google Scholar 
Halliday, W. D. et al. Underwater sound levels in the Canadian Arctic, 2014–2019. Mar. Pollut. Bull. 168, 112437 (2021).CAS 
PubMed 
Article 

Google Scholar 
Ødegaard, L., Pedersen, G. & Johnsen, E. Underwater Noise From Wind At the High North Love Ocean Observatory. UACE 2019 Conf. Proc. 359–366 (2019).Zhang, G., Forland, T. N., Johnsen, E., Pedersen, G. & Dong, H. Measurements of underwater noise radiated by commercial ships at a cabled ocean observatory. Mar. Pollut. Bull. 153, 110948 (2020).CAS 
PubMed 
Article 

Google Scholar 
Maystrenko, Y. P., Olesen, O., Gernigon, L. & Gradmann, S. Deep structure of the Lofoten–Vesterålen segment of the Mid-Norwegian continental margin and adjacent areas derived from 3-D density modeling. J. Geophys. Res. Solid Earth 122, 1402–1433 (2017).ADS 
Article 

Google Scholar 
Gillespie, D. et al. PAMGUARD: Semiautomated, open source software for real-time acoustic detection and localization of cetaceans. J. Acoust. Soc. Am. 125, 2547–2547 (2009).ADS 
Article 

Google Scholar 
Hollander, M. & Wolfe, D. A. Nonparametric Statistical Methods (Wiley, 1973).MATH 

Google Scholar 
Vogel, E. F. et al. Killer whale movements on the Norwegian shelf are associated with herring density. Mar. Ecol. Prog. Ser. 665, 217–231 (2021).ADS 
Article 

Google Scholar 
Garcia, H. A. et al. Temporal-spatial, spectral, and source level distributions of fin whale vocalizations in the Norwegian Sea observed with a coherent hydrophone array. ICES J. Mar. Sci. 76, 268–283 (2019).Article 

Google Scholar 
Davis, G. E. et al. Exploring movement patterns and changing distributions of baleen whales in the western North Atlantic using a decade of passive acoustic data. Glob. Chang. Biol. 26, 4812–4840 (2020).ADS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Risch, D. et al. Minke whale acoustic behavior and multi-year seasonal and diel vocalization patterns in Massachusetts Bay, USA. Mar. Ecol. Prog. Ser. 489, 279–295 (2013).ADS 
Article 

Google Scholar 
Le Tixerant, M., Le Guyader, D., Gourmelon, F. & Queffelec, B. How can Automatic Identification System (AIS) data be used for maritime spatial planning?. Ocean Coast. Manag. 166, 18–30 (2018).Article 

Google Scholar 
Team, R. C. R: A Language and Environment for Statistical Computing. (R Foundation for Statistical Computing, 2020).Sumner, M. D. The Tag Location Problem. 133 (2011).Sumner, M. D., Wotherspoon, S. J. & Hindell, M. A. Bayesian estimation of animal movement from archival and satellite tags. PLoS ONE 4, e7324 (2009).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Halliday, W. D. et al. The coastal Arctic marine soundscape near Ulukhaktok, Northwest Territories, Canada. Polar Biol. 43, 623–636 (2020).Article 

Google Scholar 
Ezzet, F. & Pinheiro, J. Linear, generalized linear, and nonlinear mixed effects models. Pharm. Sci. Quant. Pharmacol. 1, 103–135. https://doi.org/10.1002/9780470087978.ch4 (2006).Article 

Google Scholar 
Mazerolle, M. J. AICcmodavg: Model Selection and Multimodel Inference Based on (Q)AIC(c). (2020).Wickham, H. ggplot2: Elegant Graphics for Data Analysis (Springer, 2016).MATH 
Book 

Google Scholar 
Pante, E. & Simon-Bouhet, B. marmap: A package for importing, plotting and analyzing bathymetric and topographic data in R. PLoS ONE 8, e73051 (2013).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Wessel, P. & Walter, H. F. S. A global self-consistent, hierarchical, high-resolution shoreline database. J. Geophys. Res. 101, 8741–8743 (1996).ADS 
Article 

Google Scholar 
Sueur, J., Aubin, T. & Simonis, C. Equipment review: Seewave, a free modular tool for sound analysis and synthesis. Bioacoustics 18, 213–226 (2008).Article 

Google Scholar 
The Mathworks Inc. MATLAB (R2019a). (MathWorks Inc., 2019).Merchant, N. D. et al. Measuring acoustic habitats. Methods Ecol. Evol. 6, 257–265 (2015).PubMed 
PubMed Central 
Article 

Google Scholar  More

IPCC (2013) Climate Change: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change [Stocker TF, Qin D Qin, Plattner G-K, Tignor M, Allen SK, Boschung J, Nauels A, Xia Y, Bex V & Midgley PM (eds.)]. Cambridge University Press, Cambridge, United Kingdom and New York, NY, USA, 1535 pp (2013).Graham, L.P. Projections of Future Anthropogenic Climate Change [in:] Assessment of Climate Change for the Baltic Sea Basin. Regional Climate Studies. Bolle H.J., Menenti M., Rasool I. Series Editors Springer-Verlag Berlin Heidelberg s.133–220 (2008).Früchtenich, E., Bock, J., Feucht, V., Früchtenich W. Reactions of three European oak species ( Q. robur, Q. petraea and Q. ilex ) to repetitive summer drought in sandy soil. Trees, Forests and People 5: 100093 (2021).Gray, S. B. & Brady, S. M. Plant developmental responses to climate change. Dev. Biol. 419, 64–77 (2016).CAS 
Article 

Google Scholar 
Willliams, A. & De Vries, F. T. Plant root exudation under drought: Implications for ecosystem functioning. New Phytol. 225, 1899–1905 (2019).Article 

Google Scholar 
Canarini, A., Merchant, A. & Dijkstra, F. A. Drought effects on Helianthus annuus and Glycine max metabolites: From phloem to root exudates. Rhizosphere 2, 85–97 (2016).Article 

Google Scholar 
De Vries, F. T. et al. Changes in root-exudate-induced respiration reveal a novel mechanism through which drought affects ecosystem carbon cycling. New Phytol. 224, 132–145 (2019).Article 

Google Scholar 
Phillips, R. P., Finzi, A. C. & Bernhardt, E. S. Enhanced root exudation indu ces microbial feedbacks to N cycling in a pine forest under long-term CO2 fumigation. Ecol. Lett. 14, 187–194 (2011).Article 

Google Scholar 
Meier, I. C. et al. Root exudation of mature beech forests across a nutrient availability gradient: The role of root morphology and fungal activity. New Phytol. 226, 583–594 (2020).CAS 
Article 

Google Scholar 
Gianfreda, L. Enzymes of importance to rhizosphere processes. J. Soil Sci. Plant Nutr. 15, 283–306 (2015).
Google Scholar 
Małek, S., Ważny, R., Błońska, E., Jasik, M. & Lasota, J. Soil fungal diversity and biological activity as indicators of fertilization strategies in a forest ecosystem after spruce disintegration in the Karpaty Mountains. Sci. Total Environ. 751, 142335 (2021).ADS 
Article 

Google Scholar 
Zuccarini, P., Asensio, D., Ogaya, R., Sardans, J. & Penuelas, J. Effects of seasonal and decadal warming on soil enzymatic activity in a P-deficient Mediterranean shrubland. Glob. Change Biol. 26, 3698–3714 (2019).ADS 
Article 

Google Scholar 
Sing, S. et al. Soil organic carbon cycling in response to simulated soil moisture variation under field conditions. Sci. Rep. 11, 10841 (2021).ADS 
Article 

Google Scholar 
Sardans, J. & Penuelas, J. Drought decreases soil enzyme activity in a Mediterranean Quercus ilex L. forest. Soil Biol. Biochem. 37, 455–461 (2005).CAS 
Article 

Google Scholar 
Czúcz, B., Gálhidy, L. & Mátyás, C. Present and forecasted xeric climatic limits of beech and sessile oak distribution at low altitudes in Central Europe. Ann. For. Sci. 68, 99–108. https://doi.org/10.1007/s13595-011-0011-4 (2011).Article 

Google Scholar 
Sáenz-Romero, C. et al. Adaptive and plastic responses of Quercus petraea populations to climate across Europe. Glob. Change Biol. 23, 2831–2847 (2018).ADS 
Article 

Google Scholar 
Regulation of the Minister of the Environment. Detailed requirements for the forest reprudactive material (Dz. U. Nr 31, poz. 272) (in Polish) (2004).Phillips, R. P., Erlitz, Y., Bier, R. & Bernhardt, E. S. New approach for capturing soluble root exudates in forest soils. Funct. Ecol. 22, 990–999. https://doi.org/10.1111/j.1365-2435.2008.01495.x (2008).Article 

Google Scholar 
Ostonen, I., Lõhmus, K. & Lasn, R. The role of soil conditions in fine root ecomorphology in Norway spruce (Picea abies (L.) Karst.). Plant Soil 208, 283–292 (1999).CAS 
Article 

Google Scholar 
Pritsch, K. et al. A rapid and highly sensitive method for measuring enzyme activities in single mycorrhizal tips using 4-methylumbelliferone-labelled fluorogenic substrates in a microplate system. J. Microbiol. Methods 58, 233–241 (2004).CAS 
Article 

Google Scholar 
Sanaullah, M., Razavi, B. S., Blagodatskaya, E. & Kuzyakov, Y. Spatial distribution and catalytic mechanisms of β-glucosidase activity at the root-soil interface. Biol. Fert. Soils 52, 505–514 (2016).CAS 
Article 

Google Scholar 
R Core Team (2021). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL http://www.R-project.org/.Hartmann, H. Will a 385million year-struggle for light become a struggle for water and for carbon?–how trees may cope with more frequent climate change-type drought events. Glob. Change Biol. 17, 642–655 (2011).ADS 
Article 

Google Scholar 
Brunner, I., Herzog, C., Dawes, M. A., Arend, M. & Sperisen, C. How tree roots respond to drought. Front. Plant Sci. 6, 547 (2015).Article 

Google Scholar 
Markesteijn, L. & Poorter, L. Seedling root morphology and biomass allocation of 62 tropical tree species in relation to drought- and shade-tolerance. J. Ecol. 97, 311–325 (2009).Article 

Google Scholar 
Poorter, L. & Markesteijn, L. Seedling Traits Determine Drought Tolerance of Tropical Tree Species. Biotropica 40, 321–331 (2008).Article 

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

Google Scholar 
Lozano, Y. M., Aguilar-Triqueros, C. A., Flaig, I. C. & Rillig, M. C. Root trait responses to drought are more heterogeneous than leaf trait responses. Funct. Ecol. 34, 2224–2235 (2020).Article 

Google Scholar 
De Vries, F. T., Brown, C. & Stevens, C. J. Grassland species root response to drought: consequences for soil carbon and nitrogen availability. Plant Soil 409, 297–312 (2016).Article 

Google Scholar 
Sell, M. et al. Responses of fine root exudation, respiration and morphology in three early successional ree species to increased air humidity and different soil nitrogen sources. Tree Physiol. 42, 557–569. https://doi.org/10.1093/treephys/tpab118 (2021).Article 

Google Scholar 
Karlowsky, S. et al. Drought-induced accumulation of root exudates supports post-drought recovery of microbes in mountain grassland. Front. Plant Sci. 9, 1593. https://doi.org/10.3389/fpls.2018.01593 (2018).Article 
PubMed 
PubMed Central 

Google Scholar 
Fuchslueger, L., Bahn, M., Fritz, K., Hasibeder, R. & Richter, A. Experimental drought reduces the transfer of recently fixed plant carbon to soil microbes and alters the bacterial community composition in a mountain meadow. New Phytol. 201, 916–927. https://doi.org/10.1111/nph.12569 (2014).CAS 
Article 
PubMed 

Google Scholar 
Gargallo-Garriga, A. et al. Root exudate metabolomes change under drought and show limited capacity for recovery. Sci. Rep. 8, 12696. https://doi.org/10.1038/s41598-018-30150-0 (2018).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Zhang, X., Dippold, M. A., Kuzyakov, Y. & Razavi, B. S. Spatial pattern of enzyme activities depends on root exudate composition. Soil Biol. Biochem. 133, 83–93. https://doi.org/10.1016/j.soilbio.2019.02.010 (2019).CAS 
Article 

Google Scholar 
Hommel, R. et al. Impact of interspecific competition and drought on the allocation of new assimilates in trees. Plant Biol. 18, 785–796. https://doi.org/10.1111/plb.12461 (2016).CAS 
Article 
PubMed 

Google Scholar  More