Philippa Kaur

The chemical components analysis of different extracts of A. argyi
In our preliminary study, we accidentally found that A. argyi powder significantly inhibited the germination and reduced the varieties and biomass of weeds in the field, when it was applied as a fertilizer originally. Therefore, we speculated that certain allelochemicals present in A. argyi might inhibit the growth of weeds. To investigated the possible allelochemicals in A. argyi, three solvents (water, 50% ethanol and pure ethanol) were used to extract the metabolites in A. argyi leaves. The three type of extracts were analysed by UPLC-Q-TOF-MS and the components were confirmed by comparison with synthetic standards and MS data in literatures9,10,11. As shown in Table 1 and supplement Fig. 1, we have identified a total of 29 components in A. argyi. Six main compound mass signals were identified in the water extract: caffeic acid, schaftoside, 4-caffeoylquinic acid, 5-caffeoylquinic acid, 3,5-dicaffeoylquinic acid and 3-caffeoylquinic acid. The main compounds of the 50% ethanol extract were 4,5-dicaffeoylquinic acid, 3-caffeoylquinic acid, schaftoside, rutin, kaempferol 3-rutinoside, 3,4-dicaffeoylquinic acid, 3,5-dicaffeoylquinic acid, 3-caffeoy,1-p-coumaroylquinic acid, 1,3,4-tri-caffeoylquinic acid and eupatilin. The metabolites with higher contents in the pure ethanol extract were eupatilin, jaceosidin and casticin. Among these compounds, caffeic acid is very unique in water extract. Higher contents of schaftoside, 4-caffeoylquinic acid and 3-caffeoylquinic acid were observed in water extract and 50% ethanol extract, but very low concentrations were detected in the pure ethanol extract. 3,4-dicaffeoylquinic acid, jaceosidin, eupatilin and casticin were present at higher concentrations in the 50% ethanol extract and pure ethanol extract, but were detected at very low concentrations or were absent in the water-soluble extract. In a word, we have preliminarily identified the chemical components of different extracts.
Table 1 The chemical composition of different solvent extracts of A. argyi.
Full size table

Comparison of the allelopathic effects of different extracts of A. argyi
To explore the allelopathic effects of three different extracts of A. argyi, seed germination and seedling growth of B. pekinensis, L. sativa and O. sativa were investigated after treatment of A. argyi powder extracts. The results showed that the allelopathic inhibition increased in a concentration dependent manner. When seeds were incubated with extracts in a range of concentrations, the water-soluble extract of A. argyi powder exerted an extremely significant inhibitory effect on the germination index of all the three plants (Fig. 1a,b). While the 50% ethanol extract also showed striking allelopathic inhibitory effects on the germination index of B. pekinensis and L. sativa, but moderately inhibitory effects on O. sativa (Fig. 1c,d). Similarly, the pure ethanol extract only showed powerful inhibitory effects on the germination index of B. pekinensis and L. sativa, but no effects on O. sativa (Fig. 1e,f). Additionally, the water-soluble extract of A. argyi powder displayed extremely inhibition of the biomass of the three plants (Fig. 2a), while the 50% ethanol extract also exerted extremely significant allelopathic inhibitory effects on the biomass of B. pekinensis and L. sativa but inhibited O. sativa moderately (Fig. 2b). However, the pure ethanol extract exerted inhibitory effects on the biomass of these three plants only in high concentrations (Fig. 2c). Based on these results, the allelopathic intensity of the three different extracts of A. argyi was in the order of water-soluble extract  > 50% ethanol extract  > pure ethanol extract.
Figure 1

The different solvent exracts of A. argyi: (a,b) the water-soluble extract, (c,d) the 50% ethanol extract,and (e,f) the pure ethanol extract exert allelopathic effects on germination index of different plants. (n = 3,*P  germination index  > biomass  > germination rate  > root length  > stem length. All allelopathy response indexes reached -1.00 when plants were treated with 150 mg/ml extract. For O. sativa, the six physiological indexes also could be inhibited by a low concentration of extract (50 mg/ml), but the changes were not as obvious as the changes in B. pekinensis and L. sativa. The intensity of inhibition on the six indexes was root length  > stem length  > biomass  > germination index  > germination speed index  > germination rate. When O. sativa seeds treated with 100 mg/ml of extract, the allelopathic response index of root length and stem length were -1.00. The germination rate, germination speed index, germination index and biomass were -0.79, -0.91, -0.91 and -0.84, respectively, under the treatment with 150 mg/ml of extract. In brief, according to the comprehensive allelopathy index of the 6 indicators , the order in which they were sensitive to water-soluble extract of A. argyi were B. pekinensis (Cruciferae)  > L. sativa (Compositae)  > O. sativa (Gramineae).
Figure 3

The water-soluble extract of A. argyi inhibits the germination and growth of Brassica pekinensis, Lactuca sativa and Oryza sativa. Specific performance is in a series of indicators: (a) the germination rate, (b) the germination speed index, (c) the root length, (d) the stem length. (n = 3,*P  More

Viromes outperform total metagenomes in the recovery of viral sequences from complex soil communities
To determine the extent to which viral sequences were enriched and bacterial and archaeal sequences were depleted in viromes, relative to total metagenomes, we performed a series of analyses to compare these two approaches. After quality filtering, total metagenomes yielded an average of 8,741,015 paired reads per library for April samples and 14,551,631 paired reads for August samples, while viromes yielded an average of 9,519,518 and 5,770,419 paired reads in April and August, respectively (Fig. 1A and Supplementary Table 2). Viromes displayed a significant depletion of bacterial and archaeal sequences, as evidenced by fewer reads classified as 16S rRNA gene fragments: 0.006% of virome reads, compared to 0.042% of reads in total metagenomes (Fig. 1B). Moreover, taxonomic classification of the recovered 16S rRNA gene reads revealed clear differences in the microbial profiles associated with each approach: total metagenomes were significantly enriched in Acidobacteria, Actinobacteria, Firmicutes, and Thaumarchaeota, whereas viromes were significantly enriched in Armatimonadetes, Saccharibacteria, and Parcubacteria (Supplementary Fig. 2A, B). These last two taxa belong to the candidate phyla radiation and are typified by small cells [59,60,61], which would be more likely to pass through the 0.22-um filter that we used for viral particle purification [37, 62]. Although we acknowledge that taxon-specific differences in 16S rRNA gene copy numbers could theoretically account for some of the observed differences in absolute numbers of reads assigned to 16S rRNA genes between viromes and total metagenomes [63], in the context of subsequent analyses (see below), the most parsimonious interpretation is that both the abundances and types of bacterial and archaeal genomic content differed between the two datasets.
Fig. 1: Differences in sequence composition and assembly performance between total metagenomes and viromes.

A Sequencing depth distribution across profiling methods and time points. The y-axis displays the number of paired reads in each library after quality trimming and adapter removal. Boxes display the median and interquartile range (IQR), and data points further than 1.5x IQR from box hinges are plotted as outliers. B Percent of reads classified as 16S rRNA gene fragments in the set of quality trimmed reads; the distribution of data within boxes, whiskers, and outliers is as in A. C Sequence complexity as measured by the frequency distribution of a representative set of k-mers (k = 31) detected in each library. The x-axis displays occurrence, i.e., the number of times a particular k-mer was found in a library, while the y-axis shows the number of k-mers that exhibited a specific occurrence. D Length distribution of contigs assembled from each library (min. length = 2Kbp). White dots represent the N50 of each assembly, and green squares display the viral enrichment, as measured by the percent of contigs classified as putatively viral by DeepVirFinder and/or VirSorter. Total MG = total metagenome.

Full size image

To assess differences in sequence complexity between the two profiling methods, we calculated the k-mer frequency spectrum for each library (Fig. 1C). Relative to viromes, total metagenomes displayed an increased number of singletons (k-mers observed only once) and an overall tendency toward lower k-mer occurrences, indicating that size-fractionating our soil communities reduced sequence complexity. These differences in sequence complexity translated into notable contrasts in the quality of de novo assemblies obtained from individual libraries (Fig. 1D), while viromes yielded 800 Mbp of assembled sequences across 169,421 contigs (250 Mbp assembled in ≥10 Kbp contigs), total metagenomes produced only 65 Mbp across 22,951 contigs (1.5 Mbp assembled in ≥10 Kbp contigs). The improved assembly quality from the viromes was despite lower sequencing throughput relative to total metagenomes, particularly for the August samples (Fig. 1A). Using DeepVirFinder [48] and VirSorter [47] to mine assemblies for viral contigs, we found that 52.4% of virome contigs and only 2.2% of total metagenome contigs were identified as viral. Together, these results show that our laboratory methods for removing contamination from cells and free DNA reduced genomic signatures from cellular organisms, substantially improved sequence assembly, and successfully enriched the viral signal in soil viromes relative to total metagenomes.
Viromes facilitate exploration of the rare virosphere
To remove redundancy in our assemblies, we clustered all 192,372 contigs into a set of 105,909 representative contigs (global identity threshold = 0.95). Following current standards to define viral populations (vOTUs) [51, 52], we then screened all nonredundant ≥10 Kbp contigs for viral signatures. We identified 4065 vOTUs with a median sequence length of 17,870 bp (max = 259,025 bp) and a median gene content of 27 predicted ORFs (max = 421 ORFs). To profile the viral communities in our samples, we mapped reads against this database of nonredundant vOTU sequences (≥90% average nucleotide identity, ≥75% coverage over the length of the contig). On average, 0.04% of total metagenomic reads and 23.4% of viromic reads were mapped to vOTUs (Supplementary Fig. 3A). One August virome sample (CS-H) had particularly low sequencing throughput and low vOTU recovery (Fig. 1A and Supplementary Fig. 3B) and was discarded from downstream analyses.
In total, 2961 vOTUs were detected through read mapping in at least one sample. Of these, 2864 were exclusively found in viromes, 94 in both viromes and total metagenomes, and three in total metagenomes alone. Thus, viromes were able to recover 30 times as many viral populations as total metagenomes, even when vOTUs assembled from viromes were part of the reference set for read mapping. Notably, the three vOTUs exclusively detected in total metagenomes were only present in one metagenome from April that did not have a successful paired virome (Supplementary Fig. 4). Considering that all other vOTUs detected in total metagenomes were detected in at least one virome, it seems possible that the corresponding virome could have contained these vOTUs if sequencing had been successful. Consistent with capturing a representative amount of viral diversity from the viromes but not total metagenomes, our sampling effort was sufficient to approach a richness asymptote in vOTU accumulation curves derived from viromes but not total metagenomes (Fig. 2A).
Fig. 2: Viral richness, abundance, and occupancy patterns captured by viromes compared to total metagenomes.

A Accumulation curves of vOTUs in total metagenomes (red, n = 16) and viromes (blue, n = 14). Dots represent cumulative richness at each sampling effort across 100 permutations of sample order; the overlaid line displays the mean cumulative richness. The right graph includes the same total metagenomic data as the left graph, zoomed in along the y-axis. B Abundance-occupancy data based on vOTU profiles derived from viromes. Data in blue are from vOTUs detected only in viromes, and data in red are from vOTUs detected in both viromes and total metagenomes. Bottom left: dots represent the mean relative abundance (x-axis) and occupancy (percent of samples in which a given vOTU was detected, y-axis) that individual vOTUs displayed in viromes within a collection time point (April or August). Thus, vOTUs detected in both time points are represented twice. Red dots highlight the set of vOTUs shared between total metagenomes and viromes. Top: density curves showing the distribution of relative abundances for all vOTUs detected in viromes (blue) or the subset of vOTUs detected in viromes and total metagenomes (red). Bottom right: percent of vOTUs (x-axis) found at each occupancy level (y-axis). Red bars highlight the percent of vOTUs detected in both profiling methods. C Euler diagram displaying the overlap in detection for each vOTU (n = 2961) across profiling methods. Red vOTUs were detected by both profiling methods, and three vOTUs were detected exclusively in total metagenomes. Total MG = total metagenome.

Full size image

To examine the distribution of vOTUs along the abundance-occupancy spectrum, we compared mean relative abundances of vOTUs against the number of samples in which each vOTU was detected. Given the contrasting experimental conditions between the April and August collections, we performed this analysis within each time point. In viromes, highly abundant vOTUs tended to be recovered in the majority of samples (i.e., they displayed high occupancies), while rare vOTUs were typically recovered in only a few samples (Fig. 2B, C and Supplementary Fig. 5A), a trend usually observed in microbial communities [64]. Furthermore, more than 30% of vOTUs were found in all sampled plots, indicating the presence of a sizable core virosphere distributed throughout the field. In contrast, the distribution of 16S rRNA gene OTUs in viromes leaned toward lower occupancies (Supplementary Fig. 5B) as expected from the significant depletion of cellular genomes upon size fractionation (Fig. 1B). On the other hand, more than 80% of vOTUs in total metagenomes were detected only once (Supplementary Fig. 5C), despite the widespread distribution displayed by the 16S rRNA gene OTUs identified in the same samples (Supplementary Fig. 5D), suggesting a sparse recovery of viral diversity compared to a more complete recovery of bacterial and archaeal diversity in total metagenomes.
Inspecting the abundance-occupancy patterns for the 94 vOTUs detected in both viromes and total metagenomes revealed that vOTUs recovered from total metagenomes were among the most abundant and ubiquitous in virome profiles (Fig. 2B), indicating that total soil metagenomes were more likely to miss the rare virosphere. Notably, comparing the relative abundances of vOTUs across paired total metagenomes and viromes showed that their abundance-based ranks were not always preserved (Supplementary Fig. 4). While this discrepancy could stem from methodological challenges associated with virome preparations (e.g., differential adsorption of viruses to the soil matrix could have impacted their resuspension and recovery, therefore affecting the relative abundances of the associated vOTUs), it is also likely that total metagenomes were more susceptible to subsampling biases as evidenced by the sparse and inconsistent recovery of vOTUs exhibited by this profiling method (Supplementary Fig. 5C).
Viromes reveal a diverse taxonomic landscape
To examine the taxonomic spread covered by our vOTUs, we compared them against the RefSeq prokaryotic virus database using vConTACT2, a network-based method to classify viral contigs [49]. Under this approach, vOTUs are grouped by shared predicted protein content into taxonomically informative viral clusters (VCs) that approximate viral genera. Of the 2961 vOTUs, 1712 were confidently assigned to VCs, while the rest were only weakly connected to other clusters (outliers, 784 vOTUs) or shared no genus-level predicted protein content with any other contigs (singletons, 465 vOTUs) (Supplementary Table 3). Only 130 vOTUs were grouped with RefSeq genomes, indicating that this dataset has substantially expanded known viral taxonomic diversity (Fig. 3A). Subsetting the vOTUs detected by each profiling method showed that viromes captured a more taxonomically diverse set of viruses: 1711 vOTUs detected in viromes were assigned to 533 VCs, while 68 vOTUs detected in total metagenomes (67 of which were also detected in viromes) were assigned to 54 VCs (53 of which were also detected in viromes). Thus, any potential biases in the types of viruses recovered through soil viromics (e.g., through preferential recovery of certain viral taxonomic groups from the soil matrix) were not immediately obvious. Any such biases were either eclipsed by the much greater taxonomic diversity of viruses recovered in viromes, relative to total metagenomes, and/or they also apply to total metagenomes, at least for the viruses and soils examined here.
Fig. 3: Taxonomic diversity and predicted hosts of viral populations (vOTUs) identified in viromes compared to total metagenomes.

A Gene-sharing network of vOTUs detected in viromes alone (blue nodes), total metagenomes (red nodes; nodes outlined in white were also detected in viromes, nodes outlined in black were detected exclusively in total metagenomes), and RefSeq prokaryotic virus genomes (gray nodes). Edges connect contigs or genomes with a significant overlap in predicted protein content. Only vOTUs and genomes assigned to a viral cluster (VC) are shown. Accompanying bar plots indicate the number of distinct VCs detected in total metagenomes and viromes (VCs detected in both profiling methods are counted twice, once per bar plot). B, C Subnetwork of all RefSeq genomes and co-clustered vOTUs. Colored nodes indicate the virus family (B) or the associated host phylum (C) of each RefSeq genome. Bar plots display the number of vOTUs classified as each predicted family (B) or host phylum (C) across total metagenomes and viromes (vOTUs detected in both profiling methods are counted twice, once per bar plot). Total MG = total metagenome.

Full size image

Most of the 130 vOTUs clustered with RefSeq viral genomes could be taxonomically classified at the family level (Fig. 3B). Podoviridae was the most highly represented family, followed by Siphoviridae and Myoviridae. Myoviridae were only detected in viromes, not total metagenomes, further confirming that viromes do not seem to exclude viral groups relative to total metagenomes—if anything, the opposite seems to be true. Among the Siphoviridae clusters, we could further identify three vOTUs as belonging to the genus Decurrovirus, which are phages of Arthrobacter, a genus of Actinobacteria common in soil [65,66,67]. Because the genome network was highly structured by host taxonomy (Fig. 3A, [68]), we used consistent host signatures among RefSeq viruses in the same VC to assign putative hosts to vOTUs in VCs shared with RefSeq genomes. Most such vOTUs were putatively assigned to Proteobacteria, Actinobacteria, or Bacteroidetes hosts, and a few were linked to Firmicutes. Interestingly, these bacterial phyla were among the most abundant taxa in the 16S rRNA gene profiles derived from the total metagenomes from these soils (Supplementary Fig. 2A).
Although soil viromes and total metagenomes have been compared [62] and their presumed advantages and disadvantages have been reviewed [10], here a comprehensive comparison of results from both profiling approaches applied to the same samples showed that soil viromes recover richer (Fig. 2A) and more taxonomically diverse (Fig. 3) soil viral communities than total metagenomes.
Compositional patterns of agricultural soil viral communities and their ecological drivers
Since viromes vastly outperformed total metagenomes in capturing the viral diversity in our samples, we focused on viromes to explore the compositional relationships among viral communities. To assess the impact of each individual experimental factor on beta diversity, we performed separate permutational multivariate analyses of variance (PERMANOVA) on Bray–Curtis dissimilarities (Supplementary Table 4). Collection time point had a significant effect (R2 = 0.50, p = 0.001), but biochar (R2 = 0.19, p = 0.58) and nitrogen (R2 = 0.12, p = 0.75) treatments did not (only samples from the August time points, after nitrogen amendments, were considered for the nitrogen analysis). Additionally, to determine if the location of each sampled plot had an impact on community composition, we tested the effect of plot position along the West–East (W–E) and South–North (S–N) axes of the field (Supplementary Table 4). Viral communities displayed a significant spatial gradient along the W–E axis (R2 = 0.17, p = 0.046) but not the S–N axis (R2 = 0.10, p = 0.20). Given the significant spatiotemporal structuring in our samples, we performed an additional PERMANOVA to examine the effect of biochar while accounting for these factors (Supplementary Table 5) and detected a significant effect (R2 = 0.19, p = 0.012) only when both collection time point and W–E gradient were part of the model. We did not detect a significant effect of nitrogen treatment, even after accounting for the W–E gradient in the August samples (Supplementary Table 5).
To assess whether the bacterial and archaeal communities displayed similar compositional trends and could therefore potentially explain patterns in viral community composition, we attempted to generate metagenome assembled genomes (MAGs) from our total metagenomes. However, the low quality of total metagenomic assemblies (Fig. 1D) precluded MAG reconstruction (19 MAGs with a median completeness of 30.3), so instead we used 16S rRNA gene profiles recovered from total metagenomes (Supplementary Fig. 2A). Although 16S rRNA genes accounted for More

1.
Saunders JK, Fuchsman CA, McKay C, Rocap G. Complete arsenic-based respiratory cycle in the marine microbial communities of pelagic oxygen-deficient zones. Proc Natl Acad Sci USA. 2019;116:9925–30.
CAS  PubMed  Article  PubMed Central  Google Scholar 
2.
Callbeck CM, Lavik G, Ferdelman TG, Fuchs B, Gruber-Vodicka HR, Hach PF, et al. Oxygen minimum zone cryptic sulfur cycling sustained by offshore transport of key sulfur oxidizing bacteria. Nat Commun. 2018;9:1729.
PubMed  PubMed Central  Article  CAS  Google Scholar 

3.
Garcia-Robledo E, Padilla CC, Aldunate M, Stewart FJ, Ulloa O, Paulmier A, et al. Cryptic oxygen cycling in anoxic marine zones. Proc Natl Acad Sci USA. 2017;114:8319–24.
CAS  PubMed  Article  PubMed Central  Google Scholar 

4.
Sun X, Kop LFM, Lau MCY, Frank J, Jayakumar A, Lücker S, et al. Uncultured Nitrospina-like species are major nitrite oxidizing bacteria in oxygen minimum zones. ISME J. 2019;13:2391–402.
CAS  PubMed  PubMed Central  Article  Google Scholar 

5.
Tsementzi D, Wu J, Deutsch S, Nath S, Rodriguez-R LM, Burns AS, et al. SAR11 bacteria linked to ocean anoxia and nitrogen loss. Nature. 2016;536:179–83.
CAS  PubMed  PubMed Central  Article  Google Scholar 

6.
Thamdrup B, Steinsdóttir HGR, Bertagnolli AD, Padilla CC, Patin NV, Garcia-Robledo E, et al. Anaerobic methane oxidation is an important sink for methane in the ocean’s largest oxygen minimum zone. Limnol Oceanogr. 2019;64:2569–85.
CAS  Article  Google Scholar 

7.
Breitburg D, Levin LA, Oschlies A, Grégoire M, Chavez FP, Conley DJ, et al. Declining oxygen in the global ocean and coastal waters. Science. 2018;359:eaam7240.
PubMed  Article  CAS  PubMed Central  Google Scholar 

8.
Brown CT, Hug LA, Thomas BC, Sharon I, Castelle CJ, Singh A, et al. Unusual biology across a group comprising more than 15% of domain Bacteria. Nature. 2015;523:208–U173.
CAS  PubMed  Article  PubMed Central  Google Scholar 

9.
Castelle CJ, Wrighton KC, Thomas BC, Hug LA, Brown CT, Wilkins MJ, et al. Genomic expansion of domain archaea highlights roles for organisms from new phyla in anaerobic carbon cycling. Curr Biol. 2015;25:690–701.
CAS  PubMed  Article  PubMed Central  Google Scholar 

10.
Hug LA, Baker BJ, Anantharaman K, Brown CT, Probst AJ, Castelle CJ, et al. A new view of the tree of life. Nat Microbiol. 2016;1:16048.

11.
Mylroie JE, Carew JL, Moore AI. Blue holes: definition and genesis. Carbonates Evaporates. 1995;10:225–33.
CAS  Article  Google Scholar 

12.
Canganella F, Bianconi G, Kato C, Gonzalez J. Microbial ecology of submerged marine caves and holes characterized by high levels of hydrogen sulphide. Rev Environ Sci Biotechnol. 2007;6:61–70.

13.
Gischler E, Shinn EA, Oschmann W, Fiebig J, Buster NA. A 1500-year holocene caribbean climate archive from the blue hole, Lighthouse Reef, Belize. J Coast Res. 2008;246:1495–505.
Article  CAS  Google Scholar 

14.
Pohlman JW. The biogeochemistry of anchialine caves: progress and possibilities. Hydrobiologia. 2011;677:33–51.
CAS  Article  Google Scholar 

15.
Davis MC, Garey JR. Microbial function and hydrochemistry within a stratified anchialine sinkhole: A window into coastal aquifer interactions. Water. 2018;10:972–972.
Article  CAS  Google Scholar 

16.
Garman KM, Rubelmann H, Karlen DJ, Wu T, Garey JR. Comparison of an inactive submarine spring with an active nearshore anchialine spring in Florida. Hydrobiologia. 2011;677:65–87.

17.
Gonzalez BC, Iliffe TM, Macalady JL, Schaperdoth I, Kakuk B. Microbial hotspots in anchialine blue holes: Initial discoveries from the Bahamas. Hydrobiologia. 2011;677:149–56.
CAS  Article  Google Scholar 

18.
Seymour JR, Humphreys WF, Mitchell JG. Stratification of the microbial community inhabiting an anchialine sinkhole. Aquat Microb Ecol. 2007;50:11–24.

19.
Yao P, Wang XC, Bianchi TS, Yang ZS, Fu L, Zhang XH, et al. Carbon cycling in the world’s deepest blue hole. J Geophys Res. 2020;125:e2019JG005307.

20.
He H, Fu L, Liu Q, Fu L, Bi N, Yang Z, et al. Community Structure, abundance and potential functions of bacteria and archaea in the Sansha Yongle blue hole, Xisha, South China Sea. Front Microbiol. 2019;10:2404–2404.
PubMed  PubMed Central  Article  Google Scholar 

21.
He P, Xie L, Zhang X, Li J, Lin X, Pu X, et al. Microbial diversity and metabolic potential in the stratified Sansha Yongle Blue Hole in the South China Sea. Sci Rep. 2020;10:5949–5949.
CAS  PubMed  PubMed Central  Article  Google Scholar 

22.
DeWitt D. Submarine springs and other Karst features in offshore waters of the Gulf of Mexico and Tampa Bay, Southwest Florida Water Management District. 2003.

23.
Hu C, Muller-Karger FE, Swarzenski PW. Hurricanes, submarine groundwater discharge, and Florida’s red tides. Geophys Res Lett. 2006;33:L11601.
Google Scholar 

24.
Smith CG, Swarzenski PW. An investigation of submarine groundwater-borne nutrient fluxes to the west Florida shelf and recurrent harmful algal blooms. Limnol Oceanogr. 2012;57:471–85.
CAS  Article  Google Scholar 

25.
Vargo GA, Heil CA, Fanning KA, Dixon LK, Neely MB, Lester K, et al. Nutrient availability in support of Karenia brevis blooms on the central West Florida Shelf: What keeps Karenia blooming? Continental Shelf Res. 2008;28:73–98.
Article  Google Scholar 

26.
Walsh DA, Zaikova E, Howes CG, Song YC, Wright JJ, Tringe SG, et al. Metagenome of a versatile chemolithoautotroph from expanding oceanic dead zones. Science. 2009;326:578–82.
CAS  PubMed  Article  PubMed Central  Google Scholar 

27.
Weisberg RH, Liu YG, Lembke C, Hu CM, Hubbard K, Garrett M. The coastal ocean circulation influence on the 2018 West Florida Shelf K. brevis Red Tide Bloom. J Geophys Res Oceans. 2019;124:2501–12.
Article  Google Scholar 

28.
Quast C, Pruesse E, Yilmaz P, Gerken J, Schweer T, Yarza P, et al. The SILVA ribosomal RNA gene database project: improved data processing and web-based tools. Nucleic Acids Res. 2013;41:D590–6.
CAS  PubMed  Article  PubMed Central  Google Scholar 

29.
Parks DH, Chuvochina M, Waite DW, Rinke C, Skarshewski A, Chaumeil PA, et al. A standardized bacterial taxonomy based on genome phylogeny substantially revises the tree of life. Nat Biotechnol. 2018;36:996–1004.
CAS  PubMed  Article  PubMed Central  Google Scholar 

30.
Rodriguez RLM, Gunturu S, Tiedje JM, Cole JR, Konstantinidis KT. Nonpareil 3: fast estimation of metagenomic coverage and sequence diversity. Msystems. 2018;3: e00039-18.

31.
Baker BJ, De Anda V, Seitz KW, Dombrowski N, Santoro AE, Lloyd KG. Diversity, ecology and evolution of Archaea. Nat Microbiol. 2020;5:887–900.
PubMed  Article  CAS  PubMed Central  Google Scholar 

32.
Thiel V, Costas AMG, Fortney NW, Martinez JN, Tank M, Roden EE, et al. “Candidatus Thermonerobacter thiotrophicus,” a non-phototrophic member of the Bacteroidetes/Chlorobi with dissimilatory sulfur metabolism in hot spring mat communities. Front Microbiol. 2019;9:3159–3159.
PubMed  PubMed Central  Article  Google Scholar 

33.
Helly JJ, Levin LA. Global distribution of naturally occurring marine hypoxia on continental margins. Deep-Sea Res Part I. 2004;51:1159–68.
CAS  Article  Google Scholar 

34.
Xie LP, Wang BD, Pu XM, Xin M, He PQ, Li CX, et al. Hydrochemical properties and chemocline of the Sansha Yongle blue hole in the South China Sea. Sci Total Environ. 2019;649:1281–92.
CAS  PubMed  Article  PubMed Central  Google Scholar 

35.
Thamdrup B, Dalsgaard T, Revsbech NP. Widespread functional anoxia in the oxygen minimum zone of the Eastern South Pacific. Deep-Sea Res Part I. 2012;65:36–45.
CAS  Article  Google Scholar 

36.
Wyrtki K. The oxygen minima in relation to ocean circulation. Deep-Sea Res Oceanographic Abstr. 1962;9:11–23.
CAS  Article  Google Scholar 

37.
Ghosh W, Dam B. Biochemistry and molecular biology of lithotrophic sulfur oxidation by taxonomically and ecologically diverse bacteria and archaea. FEMS Microbiol Rev. 2009;33:999–1043.
CAS  PubMed  Article  PubMed Central  Google Scholar 

38.
Luther GW, Findlay AJ, MacDonald DJ, Owings SM, Hanson TE, Beinart RA, et al. Thermodynamics and kinetics of sulfide oxidation by oxygen: a look at inorganically controlled reactions and biologically mediated processes in the environment. Front Microbiol. 2011;2:1–9.
Article  CAS  Google Scholar 

39.
Houghton JL, Foustoukos DI, Flynn TM, Vetriani C, Bradley AS, Fike DA. Thiosulfate oxidation by Thiomicrospira thermophila: metabolic flexibility in response to ambient geochemistry. Environ Microbiol. 2016;18:3057–72.
CAS  PubMed  PubMed Central  Article  Google Scholar 

40.
Kelly DP, Shergill JK, Lu WP, Wood AP. Oxidative metabolism of inorganic sulfur compounds by bacteria. Antonie Van Leeuwenhoek Int J Gen Mol Microbiol. 1997;71:95–107.
CAS  Article  Google Scholar 

41.
Grimm F, Franz B, Dahl C. Thiosulfate and sulfur oxidation in purple sulfur bacteria. In: Dahl C, Friedrich C, editors. Microbial sulfur metabolism. Springer: Heidelberg, Germany; 2008. p. 101–16.

42.
Zopfi J, Ferdelman TG, Fossing H. Distribution and fate of sulfur intermediates – sulfite, tetrathionate, thiosulfate, and elemental sulfur – in marine sediments. In: Amend JP, Edwards KJ, Lyons TW, editors. Sulfur biogeochemistry: past and present. The Geological Society of America: Boulder, Colorado; 2004. p. 97–116.

43.
Wright JJ, Konwar KM, Hallam SJ. Microbial ecology of expanding oxygen minimum zones. Nat Rev Microbiol. 2012;10:381–94.
CAS  PubMed  Article  PubMed Central  Google Scholar 

44.
Bertagnolli AD, Stewart FJ. Microbial niches in marine oxygen minimum zones. Nat Rev Microbiol. 2018;16:723–9.
CAS  PubMed  Article  PubMed Central  Google Scholar 

45.
Hawley AK, Brewer HM, Norbeck AD, Pasǎ-Tolić L, Hallam SJ. Metaproteomics reveals differential modes of metabolic coupling among ubiquitous oxygen minimum zone microbes. Proc Natl Acad Sci USA. 2014;111:11395–400.
CAS  PubMed  Article  PubMed Central  Google Scholar 

46.
Anantharaman K, Hausmann B, Jungbluth SP, Kantor RS, Lavy A, Warren LA, et al. Expanded diversity of microbial groups that shape the dissimilatory sulfur cycle. ISME J. 2018;12:1715–28.
CAS  PubMed  PubMed Central  Article  Google Scholar 

47.
Murillo AA, Ramírez-Flandes S, DeLong EF, Ulloa O. Enhanced metabolic versatility of planktonic sulfur-oxidizing gamma-proteobacteria in an oxygen-deficient coastal ecosystem. Front Mar Sci. 2014;1:1–13.

48.
Shah V, Chang BX, Morris RM. Cultivation of a chemoautotroph from the SUP05 clade of marine bacteria that produces nitrite and consumes ammonium. ISME J. 2017;11:263–71.
CAS  PubMed  Article  PubMed Central  Google Scholar 

49.
Wirsen CO, Sievert SM, Cavanaugh CM, Molyneaux SJ, Ahmad A, Taylor LT, et al. Characterization of an autotrophic sulfide-oxidizing marine Arcobacter sp. that produces filamentous sulfur. Appl Environ Microbiol. 2002;68:316–25.
CAS  PubMed  PubMed Central  Article  Google Scholar 

50.
Luther GW, Glazer BT, Hohmann L, Popp JI, Tailefert M, Rozan TF, et al. Sulfur speciation monitored in situ with solid state gold amalgam voltammetric microelectrodes: polysulfides as a special case in sediments, microbial mats and hydrothermal vent waters. J Environ Monit. 2001;3:61–6.
CAS  PubMed  Article  PubMed Central  Google Scholar 

51.
Rozan TF, Theberge SM, Luther G. Quantifying elemental sulfur (S0), bisulfide (HS-) and polysulfides (S(x)2-) using a voltammetric method. Analyt Chim Acta. 2000;415:175–84.
CAS  Article  Google Scholar 

52.
Sievert SM, Wieringa EBA, Wirsen CO, Taylor CD. Growth and mechanism of filamentous-sulfur formation by Candidatus Arcobacter sulfidicus in opposing oxygen-sulfide gradients. Environ Microbiol. 2007;9:271–6.
CAS  PubMed  Article  PubMed Central  Google Scholar 

53.
Moussard H, Corre E, Cambon-Bonavita MA, Fouquet Y, Jeanthon C. Novel uncultured Epsilonproteobacteria dominate a filamentous sulphur mat from the 13 degrees N hydrothermal vent field, East Pacific Rise. FEMS Microbiol Ecol. 2006;58:449–63.
CAS  PubMed  Article  PubMed Central  Google Scholar 

54.
Heylen K, Vanparys B, Wittebolle L, Verstraete W, Boon N, De PV. Cultivation of denitrifying bacteria: optimization of isolation conditions and diversity study. Appl Environ Microbiol. 2006;72:2637–43.
CAS  PubMed  PubMed Central  Article  Google Scholar 

55.
Taillefert M, Bono AB, Luther GW. Reactivity of freshly formed Fe(III) in synthetic solutions and (pore)waters: voltammetric evidence of an aging process. Environ Sci Technol. 2000;34:2169–77.
CAS  Article  Google Scholar 

56.
Barco RA, Emerson D, Sylvan JB, Orcutt BN, Jacobson Meyers ME, Ramírez GA, et al. New insight into microbial iron oxidation as revealed by the proteomic profile of an obligate iron-oxidizing chemolithoautotroph. Appl Environ Microbiol. 2015;81:5927–37.
CAS  PubMed  PubMed Central  Article  Google Scholar 

57.
Garber AI, Nealson KH, Okamoto A, McAllister SM, Chan CS, Barco RA, et al. FeGenie: a comprehensive tool for the identification of iron genes and iron gene neighborhoods in genome and metagenome assemblies. Front Microbiol. 2020;11:37.
PubMed  PubMed Central  Article  Google Scholar 

58.
Canfield DE, Stewart FJ, Thamdrup B, De Brabandere L, Dalsgaard T, Delong EF, et al. A cryptic sulfur cycle in oxygen-minimum-zone waters off the Chilean Coast. Science. 2010;330:1375–8.
CAS  PubMed  Article  PubMed Central  Google Scholar 

59.
Ding J, Zhang Y, Wang H, Jian H, Leng H, Xiao X. Microbial community structure of deep-sea hydrothermal vents on the ultraslow spreading Southwest Indian Ridge. Front Microbiol. 2017;8:1012.

60.
Leon-Zayas R, Peoples L, Biddle JF, Podell S, Novotny M, Cameron J, et al. The metabolic potential of the single cell genomes obtained from the Challenger Deep, Mariana Trench within the candidate superphylum Parcubacteria (OD1). Environ Microbiol. 2017;19:2769–84.
CAS  PubMed  PubMed Central  Article  Google Scholar 

61.
Liu X, Li M, Castelle CJ, Probst AJ, Zhou Z, Pan J, et al. Insights into the ecology, evolution, and metabolism of the widespread Woesearchaeotal lineages. Microbiome. 2018;6:102–102.
PubMed  PubMed Central  Article  Google Scholar 

62.
Ortiz-Alvarez R, Casamayor EO. High occurrence of Pacearchaeota and Woesearchaeota (Archaea superphylum DPANN) in the surface waters of oligotrophic high-altitude lakes. Environ Microbiol Rep. 2016;8:210–7.
CAS  PubMed  Article  PubMed Central  Google Scholar 

63.
Suominen S, Dombrowski N, Damste JSS, Villanueva L. A diverse uncultivated microbial community is responsible for organic matter degradation in the Black Sea sulphidic zone. Environ Microbiol. 2021. https://doi.org/10.1111/1462-2920.14902.

64.
Castelle CJ, Banfield JF. Major new microbial groups expand diversity and alter our understanding of the tree of life. Cell. 2018;172:1181–97.
CAS  PubMed  Article  PubMed Central  Google Scholar 

65.
Dombrowski N, Lee JH, Williams TA, Offre P, Spang A. Genomic diversity, lifestyles and evolutionary origins of DPANN archaea. FEMS Microbiol Lett. 2019;366:fnz008.

66.
Tian RM, Ning DL, He ZL, Zhang P, Spencer SJ, Gao SH, et al. Small and mighty: adaptation of superphylum Patescibacteria to groundwater environment drives their genome simplicity. Microbiome. 2020;8:51.

67.
Vigneron A, Cruaud P, Langlois V, Lovejoy C, Culley AI, Vincent WF. Ultra-small and abundant: Candidate phyla radiation bacteria are potential catalysts of carbon transformation in a thermokarst lake ecosystem. Limnol Oceanogr Lett. 2020;5:212–20.
Article  Google Scholar 

68.
Beam JP, Becraft ED, Brown JM, Schulz F, Jarett JK, Bezuidt O, et al. Ancestral absence of electron transport chains in Patescibacteria and DPANN. Front Microbiol. 2020;11:1848.

69.
Luef B, Frischkorn KR, Wrighton KC, Holman HYN, Birarda G, Thomas BC, et al. Diverse uncultivated ultra-small bacterial cells in groundwater. Nat Commun. 2015;6:6372.

70.
Wrighton KC, Castelle CJ, Wilkins MJ, Hug LA, Sharon I, Thomas BC, et al. Metabolic interdependencies between phylogenetically novel fermenters and respiratory organisms in an unconfined aquifer. ISME J. 2014;8:1452–63.
CAS  PubMed  PubMed Central  Article  Google Scholar 

71.
Konstantinidis KT, Tiedje JM. Trends between gene content and genome size in prokaryotic species with larger genomes. Proc Natl Acad Sci USA. 2004;101:3160–5.
CAS  PubMed  Article  PubMed Central  Google Scholar 

72.
Moya A, Pereto J, Gil R, Latorre A. Learning how to live together: genomic insights into prokaryote-animal symbioses. Nat Rev Genet. 2008;9:218–29.
CAS  PubMed  Article  PubMed Central  Google Scholar 

73.
Moran NA, Plague GR. Genomic changes following host restriction in bacteria. Curr Opin Genet Dev. 2004;14:627–33.
CAS  PubMed  Article  PubMed Central  Google Scholar 

74.
Chaudhury P, Quax TEF, Albers SV. Versatile cell surface structures of archaea. Mol Microbiol. 2018;107:298–311.
CAS  PubMed  Article  PubMed Central  Google Scholar 

75.
Pohlschroder M, Esquivel RN. Archaeal type IV pili and their involvement in biofilm formation. Front Microbiol. 2015;6:190.

76.
Aylward FO, Santoro AE. Heterotrophic Thaumarchaea with small genomes are widespread in the dark ocean. mSystems. 2020;5:e00415–20.
CAS  PubMed  PubMed Central  Article  Google Scholar 

77.
Reji L, Francis CA. Metagenome-assembled genomes reveal unique metabolic adaptations of a basal marine Thaumarchaeota lineage. ISME J. 2020;14:2105–15.
CAS  PubMed  Article  PubMed Central  Google Scholar 

78.
Santoro AE, Richter RA, Dupont CL. Planktonic marine Archaea. Annu Rev Mar Sci. 2019;11:131–58.
Article  Google Scholar 

79.
Rinke C, Rubino F, Messer LF, Youssef N, Parks DH, Chuvochina M, et al. A phylogenomic and ecological analysis of the globally abundant Marine Group II archaea (Ca. Poseidoniales ord. nov.). ISME J. 2019;13:663–75.
CAS  PubMed  Article  PubMed Central  Google Scholar 

80.
Pereira O, Hochart C, Auguet JC, Debroas D, Galand PE. Genomic ecology of Marine Group II, the most common marine planktonic Archaea across the surface ocean. Microbiol Open. 2019;8:e00852.
Google Scholar 

81.
Martin-Cuadrado AB, Garcia-Heredia I, Moltó AG, López-Úbeda R, Kimes N, López-García P, et al. A new class of marine Euryarchaeota group II from the mediterranean deep chlorophyll maximum. ISME J. 2015;9:1619–34.
CAS  PubMed  Article  PubMed Central  Google Scholar 

82.
Martin-Cuadrado AB, Rodriguez-Valera F, Moreira D, Alba JC, Ivars-Martínez E, Henn MR, et al. Hindsight in the relative abundance, metabolic potential and genome dynamics of uncultivated marine archaea from comparative metagenomic analyses of bathypelagic plankton of different oceanic regions. ISME J. 2008;2:865–86.
CAS  PubMed  Article  PubMed Central  Google Scholar 

83.
Moreira D, Rodríguez-Valera F, López-García P. Analysis of a genome fragment of a deep-sea uncultivated Group II euryarchaeote containing 16S rDNA, a spectinomycin-like operon and several energy metabolism genes. Environ Microbiol. 2004;6:959–69.
CAS  PubMed  Article  PubMed Central  Google Scholar 

84.
Sforna MC, Philippot P, Somogyi A, Van Zuilen MA, Medjoubi K, Schoepp-Cothenet B, et al. Evidence for arsenic metabolism and cycling by microorganisms 2.7 billion years ago. Nat Geosci. 2014;7:811–5.
CAS  Article  Google Scholar 

85.
Meheust R, Burstein D, Castelle CJ, Banfield JF. The distinction of CPR bacteria from other bacteria based on protein family content. Nat Commun. 2019;10:4173.

86.
Luther GW, Glazer BT, Ma S, Trouwborst RE, Moore TS, Metzger E, et al. Use of voltammetric solid-state (micro)electrodes for studying biogeochemical processes: laboratory measurements to real time measurements with an in situ electrochemical analyzer (ISEA). Mar Chem. 2008;108:221–35.
CAS  Article  Google Scholar 

87.
Brendel PJ, Luther GW. Development of a gold amalgam voltammetric microelectrode for the determination of dissolved Fe, Mn, O2, and S(-II) in porewaters of marine and freshwater sediments. Environ Sci Technol. 1995;29:751–61.
CAS  PubMed  Article  PubMed Central  Google Scholar 

88.
Arar EJ, Collins GB. Method 445.0 in vitro determination of chlorophyll a and pheophytin a in marine and freshwater algae by fluorescence. Washington, DC: U.S. Environmental Protection Agency; 1997.

89.
Bran+Luebbe/Seal. Ammonia in water and seawater, in Method No G-171-96. 2005. Norderstedt, Germany.

90.
Bran+Luebbe/Seal. Nitrate and nitrite in water and seawater; total nitrogen in persulfate digests, in Metho No G-172-96. 2010. Norderstedt, Germany.

91.
Solórzano L, Sharp JH. Determination of total dissolved nitrogen in natural waters. Limnol Oceanogr. 1980;25:751–4.
Article  Google Scholar 

92.
Solórzano L, Sharp JH. Determination of total dissolved phosphorus and particulate phosphorus in natural waters. Limnol Oceanogr. 1980;25:754–8.
Article  Google Scholar 

93.
Dickson AG, Sabine CL, Christian JR. Guide to best practices for ocean CO2 measurements. PICES Special Publication 3. 2007.

94.
Murphy J, Riley JP. A modified single solution method for the determination of phosphate in natural waters. Analy Chim Acta. 1962;27:31–6.
CAS  Article  Google Scholar 

95.
Stookey LL. Ferrozine—a new spectrophotometric reagent for iron. Anal Chem. 1970;42:779–81.
CAS  Article  Google Scholar 

96.
Schindelin J, Arganda-Carreras I, Frise E, Kaynig V, Longair M, Pietzsch T, et al. Fiji: an open-source platform for biological-image analysis. Nat Methods. 2012;9:676–82.
CAS  PubMed  Article  PubMed Central  Google Scholar 

97.
Padilla CC, Bertagnolli AD, Bristow LA, Sarode N, Glass JB, Thamdrup B, et al. Metagenomic binning recovers a transcriptionally active gammaproteobacterium linking methanotrophy to partial denitrification in an anoxic oxygen minimum zone. Front Mar Sci. 2017;4:23–23.
Article  Google Scholar 

98.
Kozich JJ, Westcott SL, Baxter NT, Highlander SK, Schloss PD. Development of a dual-index sequencing strategy and curation pipeline for analyzing amplicon sequence data on the MiSeq illumina sequencing platform. Appl Environ Microbiol. 2013;79:5112–20.
CAS  PubMed  PubMed Central  Article  Google Scholar 

99.
Callahan BJ, McMurdie PJ, Rosen MJ, Han AW, Johnson AJA, Holmes SP. DADA2: High-resolution sample inference from Illumina amplicon data. Nat Methods. 2016;13:581–3.
CAS  PubMed  PubMed Central  Article  Google Scholar 

100.
Bolyen E, Rideout JR, Dillon MR, Bokulich NA, Abnet CC, Al-Ghalith GA, et al. Reproducible, interactive, scalable and extensible microbiome data science using QIIME 2. Nat Biotechnol. 2019;37:852–7.
CAS  PubMed  PubMed Central  Article  Google Scholar 

101.
Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15:550–550.
PubMed  PubMed Central  Article  CAS  Google Scholar 

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

103.
McMurdie PJ, Holmes S. Waste not, want not: why rarefying microbiome data is inadmissible. PLoS Comput Biol. 2014;10:e1003531.

104.
Willis AD, Martin BD. DivNet: estimating diversity in networked communities. bioRxiv. 2018. Available from https://www.biorxiv.org/content/10.1101/305045v1.

105.
Wickham H. Elegant graphics for data analysis. New York: Springer-Verlag; 2016.
Google Scholar 

106.
Oksanen J, Blanchet F, Friendly M, Kindt R, Legendre P, Mcglinn D, et al. vegan: Community Ecology package, in R package version 2.5-5. 2019. https://cran.r-project.org/package=vegan.

107.
Nayfach S, Pollard KS. Average genome size estimation improves comparative metagenomics and sheds light on the functional ecology of the human microbiome. Genome Biol. 2015;16:51–51.
PubMed  PubMed Central  Article  CAS  Google Scholar 

108.
Nayfach S, Bradley PH, Wyman SK, Laurent TJ, Williams A, Eisen JA, et al. Automated and accurate estimation of gene family abundance from shotgun metagenomes. PLoS Comput Biol. 2015;11:e1004573.

109.
Nayfach S, Pollard KS. Toward accurate and quantitative comparative metagenomics. Cell. 2016;166:1103–16.
CAS  PubMed  PubMed Central  Article  Google Scholar 

110.
Nurk S, Meleshko D, Korobeynikov A, Pevzner PA. MetaSPAdes: a new versatile metagenomic assembler. Genome Res. 2017;27:824–34.
CAS  PubMed  PubMed Central  Article  Google Scholar 

111.
Li D, Liu CM, Luo R, Sadakane K, Lam TW. MEGAHIT: an ultra-fast single-node solution for large and complex metagenomics assembly via succinct de Bruijn graph. Bioinformatics. 2015;31:1674–6.
CAS  PubMed  Article  PubMed Central  Google Scholar 

112.
Mikheenko A, Saveliev V, Gurevich A. MetaQUAST: evaluation of metagenome assemblies. Bioinformatics. 2016;32:1088–90.
CAS  PubMed  Article  PubMed Central  Google Scholar 

113.
Seemann T. Prokka: rapid prokaryotic genome annotation. Bioinformatics. 2014;30:2068–9.
CAS  PubMed  Article  PubMed Central  Google Scholar 

114.
Hyatt D, Chen GL, LoCascio PF, Land ML, Larimer FW, Hauser LJ. Prodigal: prokaryotic gene recognition and translation initiation site identification. BMC Bioinformatics. 2010;11:119–119.
PubMed  PubMed Central  Article  CAS  Google Scholar 

115.
James BT, Luczak BB, Girgis HZ. MeShClust: an intelligent tool for clustering DNA sequences. Nucleic Acids Res. 2018;46:e83.
PubMed  PubMed Central  Article  CAS  Google Scholar 

116.
Aramaki T, Blanc-Mathieu R, Endo H, Ohkubo K, Kanehisa M, Goto S, et al. KofamKOALA: KEGG ortholog assignment based on profile HMM and adaptive score threshold. Bioinformatics. 2019;36:2251–2.
PubMed Central  Article  CAS  Google Scholar 

117.
Boratyn GM, Thierry-Mieg J, Thierry-Mieg D, Busby B, Madden TL. Magic-BLAST, an accurate RNA-seq aligner for long and short reads. BMC Bioinformatics. 2019;20:405–405.
PubMed  PubMed Central  Article  CAS  Google Scholar 

118.
Dunivin TK, Yeh SY, Shade A. A global survey of arsenic-related genes in soil microbiomes. BMC Biol. 2019;17:45–45.
CAS  PubMed  PubMed Central  Article  Google Scholar 

119.
Wood DE, Lu J, Langmead B. Improved metagenomic analysis with Kraken 2. Genome Biol. 2019;20:257.

120.
Lu J, Breitwieser FP, Thielen P, Salzberg SL. Bracken: estimating species abundance in metagenomics data. PeerJ Comput Sci. 2017;3:e104.

121.
Parks DH, Imelfort M, Skennerton CT, Hugenholtz P, Tyson GW. CheckM: assessing the quality of microbial genomes recovered from isolates, single cells, and metagenomes. Genome Res. 2015;25:1043–55.
CAS  PubMed  PubMed Central  Article  Google Scholar 

122.
Olm MR, Brown CT, Brooks B, Banfield JF. DRep: a tool for fast and accurate genomic comparisons that enables improved genome recovery from metagenomes through de-replication. ISME J. 2017;11:2864–8.
CAS  PubMed  PubMed Central  Article  Google Scholar 

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

124.
Stamatakis A. RAxML version 8: a tool for phylogenetic analysis and post-analysis of large phylogenies. Bioinformatics. 2014;30:1312–3.
CAS  PubMed  PubMed Central  Article  Google Scholar 

125.
Letunic I, Bork P. Interactive Tree of Life (iTOL) v4: Recent updates and new developments. Nucleic Acids Res. 2019;47:p. W256–9.
Article  CAS  Google Scholar  More

Distribution, frequency, and functional implications of mutations during laboratory evolution of obligate syntrophy
We evaluated whether the selection of mutations in the same genes (i.e., “parallel evolution” [17]) had contributed to improvements in syntrophic growth of Dv and Mm across independent evolution lines, all of which started with the same ancestral clone of each organism. The goal of this analysis was to focus on generalized strategies for adaptation to syntrophy, irrespective of the culturing condition so we investigated parallelism across both U and H lines. Based on the number of mutations (normalized to gene length and genome size) in Dv and Mm across 13 evolved lines (six lines designated U for “uniform” conditions with continuous shaking and seven H lines for “heterogenous” conditions without shaking), we calculated a G-score [18] (“goodness-of-fit”, see “Methods” section [18]) to assess if the observed parallel evolution rate was higher than expected by chance. The “observed G-score” was calculated as the sum of G-scores for all genes in the genome of each organism; mean and standard deviation of “expected G-scores” were calculated through 1000 simulations of randomizing locations of observed numbers of mutations across the genome of each organism. The observed total G-score for Dv (1092.617) and Mm (805.02) was significantly larger than the expected mean G-score (Dv: 798.19 ± 14.99, Z = 19.63 and Mm: 564.83 ± 15.95, Z = 15.06), demonstrating significant parallel evolution across lines.
With the exception of five high G-score genes (DVU0597, DVU1862, DVU0436, DVU0013, and DVU2394), which were mutated during long term salt adaptation of Dv [19], mutations in other high G-score genes appeared to be putatively specific to syntrophic interactions. Altogether, 24 genes in Dv and 16 genes in Mm associated with core processes had accumulated function modulating mutations across at least 2 or more evolution lines (Fig. 2 and Supplementary Table S1). Signal transduction and regulatory gene mutations (seven in Dv and six in Mm) represented 19.9% and 27.2% of all mutations in Dv and Mm, respectively, similar to long term laboratory evolution of E. coli [18], potentially because their influence on the functions of many genes [20, 21]. We also observed missense and nonsense mutations in outer membrane and transport functions (four genes in Dv and three genes in Mm). For example, the highest G-score gene in Dv, DVU0799—an abundant outer membrane porin for the uptake of sulfate and other solutes in low-sulfate conditions [22], was mutated early across all lines, with at least two missense mutations in UE3 (S223Y) and UA3 (T242P). Notably, the regulator of the archaellum operon (MMP1718) had the highest G-score with frameshift (11 lines) and nonsynonymous coding (2 lines) mutations [23]. Similarly, two motility-associated genes of Dv (DVU1862 and DVU3227) also accumulated frameshift, nonsense and nonsynonymous mutations across 4 H and 3 U lines. Together, these observations were consistent with other laboratory evolution experiments performed in liquid media [24], suggesting that retaining motility has a fitness cost during syntrophy [25, 26].
Fig. 2: Frequency and location of high G-score mutations in Dv (A) and Mm (B) across 13 independent evolution lines.

SnpEff predicted impact of mutations* are indicated as moderate (orange circles) or high (red circles) with the frequency of mutations indicated by node size. Expected number of mutations for each gene was calculated based on the gene length and the total number of mutations in a given evolution line. Genes with parallel changes were ranked by calculating a G (goodness of fit) score between observed and expected values and indicated inside each panel. Mutations for each gene are plotted along their genomic coordinates (vertical axes) across 13 evolution lines (horizontal axes). Total number of mutations for a given gene is shown as horizontal bar plots. [*HIGH impact mutations: gain or loss of start and stop codons and frameshift mutations; MODERATE impact mutations: codon deletion, nonsynonymous in coding sequence, change or insertion of codon; low impact mutations: synonymous coding and nonsynonymous start codon].

Full size image

Consistent with our previous observation that obligate mutual interdependence drove the erosion of metabolic independence of Dv [5, 27], mutations in SR genes were among the top contributors to the total G-score in Dv (DVU2776 (74.7), DVU1295 (46.5), DVU0846 (42.9), and DVU0847 (22.3)). However, it was intriguing that DVU2776 (DsrC), which catalyzes the conversion of sulfite to sulfide, the final step in SR, accumulated function modulating but not loss-of-function mutations across six lines. The functional impact of these mutations is not clear but it is possible that these changes might alter previously suggested alternative roles for this gene, including electron confurcation for the oxidation of lactate [28], sulfite reduction, 2-thiouridine biosynthesis and possibly gene regulation [29].
Analysis of temporal appearance and combinations of mutations across evolution lines
Growth characteristics of all evolution lines improved by the 300th generation [4], and in some lines even before the appearance of SR− mutations, indicating that mutations in other genes had also contributed to improvements in syntrophy. Each evolution line had at least 8 and up to 13 out of 24 high G-score mutations in Dv, while Mm had mutations in at least 5 and up to 10 out of 16 high G-score genes. We interrogated the temporal order in which high G-score mutations were selected and the combinations in which they co-existed in each evolution line to uncover evidence for epistatic interactions in improving obligate syntrophy. Indeed, missense mutations in DsrC (DVU2776) were fixed simultaneously with the appearance of loss of function mutations in one of two sigma 54 type regulators (DVU2894, DVU2394) in lines HA2, and UR1 (P = 5.40 × 10−5). In rare instances, we also observed that some high G-score mutations co-occurred across evolution lines, e.g., two U- and one H-line consistently showed for at least two time points a mutation in DVU1283 (GalU) coexisting with mutations in DVU2394 (P = 5.04 × 10−3). More commonly, the combinations of high G-score gene mutations varied across multiple lines. In fact, no two lines possessed identical combination of high G-score gene mutations (Fig. 3A, B), and many high-frequency mutations were uniquely present or absent in different lines (Fig. 3C, D).
Fig. 3: Frequency and time of appearance of mutations through 1 K generations of laboratory evolution lines of Dv and Mm cocultures.

The heat maps display frequency of mutations in genes (rows) in Dv (A) and Mm (B) in each evolution line, ordered from early to later generations (horizontal axis). High G-score genes are shown in red font and their G-score rank is shown to the left in gray shaded box, also in red font. Bar plots above heat maps indicate total number of mutations in each generation and the color indicates impact of mutation. Use “Frequency”, “Generations”, and “Mutation impact” key below the heat maps for interpretation. Mutations that were unique to each evolution line is shown in (C, D) for Dv and Mm, respectively. E The heatmap illustrates a selective sweep across both organisms in line HS3.

Full size image

Mutations in high G-score genes appeared consistently in all evolution lines (P 80% EPD-03 vs, More

Miniature flight recorders
Honey bees regularly carry a payload of 55–65 mg6, but a mounted flight recorder should consume only a small fraction of this allowable weight to avoid significantly affecting bee behavior. The dimensions of the recorder must also be small, as the available mounting area on the bee thorax is limited. We propose a flight recorder consisting of a (2times 2times 0.3,text {mm}^3) ASIC mounted on a (3times 3times 0.4,text {mm}^3) printed circuit board, which is similar in size to 3 mm diameter, 1.5 mm tall conventional bee tags (“Queen number set,” Betterbee) and is on par with previous studies using honey bee tagging methods30. The ASIC provides most core functionality including signal detection, memory, power harvesting, and communications circuitry, and the PCB provides a magnetic backscatter coil for near-field wireless communication. The combined weight of our chip-PCB assembly is expected to be at most 10 mg, which is a small fraction of the honey bee payload (Fig. 1c). Based on previous studies, we expect this may slightly reduce foraging trip time but not significantly impact flight characteristics, which is the focus of our system31. Furthermore, future iterations of our flight recorder will have smaller size and weight, minimizing the overall impact on honey bees. Power will be harvested from sunlight, which can provide intensity greater than 1 mW/mm(^2). On-chip photovoltaics can offer power conversion efficiency on the order of 5%, supplying 50 μW of electrical power for the chip. This power budget, while low, is sufficient for the chip, since solar angle measurement and storage of data in memory are not energy-intensive operations and need only occur a few times per second. Furthermore, wireless communication for data upload is only used when the recorder is at the base station, thus allowing the base station to fully power the near-field wireless link. Additionally, IC technology is generally robust to the environmental factors likely to be encountered by honey bees. For instance, the variations in humidity level and temperature experienced by honey bees are not expected to affect chip operation. Our proposed design is a fully power-autonomous, environmentally-robust, miniature flight recorder well-suited to the task of recording honey bee activity.
The orientation of a bee during flight can be described by the yaw ((gamma)), which represents the absolute heading relative to the sun, and the angle-of-incidence (AOI) ((psi)), which represents the overhead angle between the sun and the sensor (Fig. 1b). To record flights, the chip uses ASPs to measure the AOI of sunlight and stores these measurements in on-chip memory. ASPs achieve AOI-sensitivity via a pair of diffraction gratings stacked over a photodiode (Fig. 1b), wherein the first grating induces a diffraction pattern that shifts laterally across the second grating as AOI is swept, thus passing a periodically-varying intensity of light to the photodiode. The stored AOI measurements can be downloaded to a base station upon return to the hive, and from these data the heading throughout the flight can be extracted. Assuming a constant speed of 6.5 m/s32, we can use the sequence of recorded headings to reconstruct the honey bee’s trajectory in post-processing.
The data taken by the flight recorder will be subject to measurement errors, and these errors will manifest in the reconstructed trajectory. Here, we identify and explore methods to mitigate the primary sources of error. We posit that errors will stem primarily from finite heading measurement resolution, finite sampling rate, and random fluctuations in sampling rate (jitter). Each of these error sources can be suppressed through careful chip design, but improvements in these performance variables can only be made at the expense of larger chip area. For instance, the heading measurement resolution will increase if more ASPs are used to measure AOI, but each pixel consumes significant silicon area and contributes additional data that must be stored in memory. Increasing the sampling rate and decreasing sampling rate jitter requires that more measurements be stored if flight time is unchanged, thus increasing the chip area required for memory. The size of the chip is thus inversely related to the severity of the expected measurement error, and trade-offs between chip size and achievable precision should be examined. A smaller sensor is feasible if trajectory reconstruction performance requirements are relaxed; more stringent requirements will necessitate a larger chip that will increase the burden on the bee. If the relationship between final uncertainty in the reconstructed position of the bee and the core sensor specifications is understood, then chip-level performance goals can be formed based on trajectory-level precision requirements.
To determine required heading resolution, sampling rate, and sampling rate jitter, we first describe the procedure to reconstruct trajectories. We define the timestep estimate (hat{Delta t}) as the inverse of the sampling rate, and for known flight speed v the sequence of measured heading estimates ({hat{gamma }_0, hat{gamma }_1, …, hat{gamma }_{n-1}}) can be mapped to position estimate (hat{mathbf {p}}_n = [hat{x}_n, , hat{y}_n]^T) as a function of discrete time index n via the motion model

$$begin{aligned} hat{mathbf {p}}_n = sum _{i=0}^{n-1} v, mathbf {h}(hat{gamma }_i) hat{Delta t} end{aligned}$$
(1)

where h is a unit vector pointing along the heading of the bee. Each sensor estimate of heading and timestep will be subject to errors, and by modeling these errors as additive white noise, we can evaluate a confidence region for each position estimate (hat{mathbf {p}}_n) (Fig. 2a,b). Each analog heading measurement (hat{gamma }_i) must be digitized for storage on-chip and will therefore suffer from quantization error, a form of rounding. We denote this error (epsilon _{gamma ,i}) and model it as a uniform random variable with variance (sigma ^2_gamma) on the interval (pm frac{Delta gamma }{2}), where (Delta gamma) is the heading bin width. Furthermore, the timestep estimate (hat{Delta t}) will be subject to random clock jitter that can be modeled as a Gaussian random variable (epsilon _{Delta t,i}) with variance (sigma ^2_{Delta t}). In this model, we posit for simplicity that the random error in timing scales linearly in proportion to oscillator frequency, thus maintaining a fixed ratio of (sigma _{Delta t}/hat{Delta t}). These measurement errors contribute random error to position estimate ({hat{mathbf {p}}}_n), and thus each position estimate should be viewed as a random variable (mathbf {p}_n). The confidence region surrounding ({hat{mathbf {p}}}_n) depends on the covariance matrix of ({mathbf {p}_n}), and we evaluate these terms by first using the small-angle approximation to linearize the motion model with respect to (epsilon _{gamma ,i}):

$$begin{aligned} mathbf {p}_n = sum _{i=0}^{n-1} v , mathbf {h}(hat{gamma }_i + epsilon _{gamma ,i}) (hat{Delta t}+epsilon _{Delta t, i}) &approx sum _{i=0}^{n-1} v , left( mathbf {h}(hat{gamma }_i) + mathbf {h}^perp (hat{gamma }_i)epsilon _{gamma ,i}right) (hat{Delta t}+epsilon _{Delta t, i}) = sum _{i=0}^{n-1} v , mathbf {h}(hat{gamma }_i)(hat{Delta t}+epsilon _{Delta t,i}) + sum _{i=0}^{n-1} v , mathbf {h}^perp (hat{gamma }_i)epsilon _{gamma ,i} (hat{Delta t} + epsilon _{Delta t, i}) end{aligned}$$
(2)

Figure 2

(a) Sensor measurement errors produce confidence regions surrounding each estimated position shown in an example circular trajectory. (b) Zoomed-in view of confidence region at end of flight, with principle components shown. (c) The two directed standard deviations grow throughout the duration of the trajectory, and are bounded above and below by the directed standard deviations computed from the case of the straight-line trajectory. (d–f) Both directional standard deviations characterizing the final error region will depend on all three of the core sensor specs ({Delta gamma , hat{Delta t}, sigma _{Delta t}}), and the max confidence region dimension will grow if these specs are relaxed. Plots were created in MATLAB33.

Full size image

This approximation is valid if (epsilon _{gamma ,i}) is kept small, which can be guaranteed by keeping heading bin width (Delta gamma) small. The covariance matrix of (mathbf {p}_n) is then given by

$$begin{aligned} mathrm {Cov}(mathbf {p}_n, mathbf {p}_n) = varvec{ Sigma }_n approx sum _{i=0}^{n-1} R(hat{{gamma }}_{i}) begin{bmatrix} v^2sigma ^2_{Delta t} &{} 0 \ 0 &{} v^2sigma ^2_{gamma }(sigma _{Delta t}^2 + (hat{Delta t})^2) end{bmatrix} {R^T({hat{gamma }}_i)} end{aligned}$$
(3)

where R is the standard 2 × 2 rotation matrix (Appendix: Derivation of Trajectory Precision Equation). By the Central Limit Theorem, after a sufficient number of timesteps (mathbf {p}_n) will become Gaussian distributed. Thus, the confidence region will be an ellipse with major and minor axes spanned by the eigenvectors ({mathbf {v}_1, mathbf {v}_2}) of ({varvec{Sigma }}_n). The covariances of the confidence region along each of these axes are given by the eigenvalues ({lambda _1, lambda _2}) of (varvec{Sigma }_n), and directed standard deviations can then be defined as ({sigma _1,sigma _2}). A 99% confidence region for position estimate ({mathbf {hat{p}}}_{n}) is given by an ellipse with major and minor axes lengths ({3sigma _1mathbf {v}_1, 3sigma _2mathbf {v}_2}). We therefore conclude that (3sigma _1) and (3sigma _2) are critical values defining achievable trajectory reconstruction precision.
These (3sigma)-bounds can be computed for any measured sequence of headings, but general upper and lower bounds for (3sigma _1) and (3sigma _2) across all possible trajectories can be derived from the (3sigma _1) and (3sigma _2) given by the case in which the bee flies in a straight line. In the straight-line case, the eigenvalues of (varvec{Sigma }_n) are given by n times the diagonal entries of the diagonal matrix in Eq. (3). The directed (3sigma)-bounds are then

$$begin{aligned} 3sigma _{1,n}^*= & {} 3sqrt{n} v sigma _{Delta t} end{aligned}$$
(4)

$$begin{aligned} 3sigma _{2,n}^*= & {} 3sqrt{n} v sigma _{gamma }sqrt{sigma _{Delta t}^2 + (hat{Delta t})^2} end{aligned}$$
(5)

For some values of ({sigma _{gamma }, hat{Delta t}, sigma _{Delta t}}), (3sigma _{1,n}^*) will be larger than (3sigma _{2,n}^*); for others, the converse will be true. These equations provide simple bounds on achievable reconstruction precision that are valid for any trajectory and can be computed from sensor characteristics. Since heading error (epsilon _{gamma _i}) is uniformly distributed, heading variance (sigma ^2_gamma) is defined by heading bin width (Delta gamma), and thus the reconstruction precision is defined by a core suite of sensor specifications: ({Delta gamma , hat{Delta t}, sigma _{Delta t}}). An illustration of the trajectory reconstruction process, along with confidence regions, is shown in Fig. 2, as well as the relationship between reconstruction precision and each of the core chip specs.
For a maximum specified chip sensing area, trajectory precision should be optimized through balanced allocation of area to solar AOI detection and to memory (Fig. 3). We evaluate the optimal area allocation by defining the standard deviation upper bound (3sigma ^*_{max} = max (3sigma _{1,f}^*, 3sigma _{2,f}^*)), where (3sigma _{1,f}^*) and (3sigma _{2,f}^*) are the directed (3sigma)-values computed from a straight-line trajectory that is long enough to completely fill the memory. If more area is spent on pixels for AOI detection, heading resolution can be increased, thus causing (sigma _gamma) to be reduced and correspondingly lowering (3sigma ^*_{2,f}). Conversely, if more area is spent on memory, measurements can be taken more frequently, and timestep and clock jitter can be reduced, thus reducing (hat{Delta t}) and (sigma _{Delta t}). This will cause (3sigma ^*_{1,f}) to decrease, but may cause an increase in (3sigma ^*_{2,f}) since less area is now available for heading sensors. As shown in Fig. 3, the upper bound (3sigma ^*_{max}) minimizes when the two counteracting variables (3sigma ^*_{1,f}) and (3sigma ^*_{2,f}) are equal, and this intersection point prescribes the optimal allocation of sensor area. Our proposed flight recorder features a 4 (mathrm {mm}^2) chip, which can offer sensing area of approximately 3 (mathrm {mm}^2). When this area is allocated optimally, the heading resolution is (2^circ) and timestep is approximately 240 ms at timestep jitter of (3sigma _{Delta t}/Delta t = 0.03). With these specifications, the maximum recordable trajectory length is approximately 4 km, with (3sigma) uncertainty of (pm 2.4) m.
Figure 3

(a) Area usage is optimized where (3sigma ^*_{1,f}) and (3sigma ^*_{2,f}) cross, as this design point co-minimizes the two directional standard deviations characterizing trajectory precision. (b) This design point specifies the heading resolution and number of data words in memory that minimize trajectory uncertainty given a fixed sensor area constraint. Plots were created in MATLAB33.

Full size image

Sensor calibration and noise modeling
We next examined the output from existing ASP array sensors in order to create a model for future flight recorders. These particular sensors have 96 pixels, or 24 sets of 4 oriented in 90° angles. Specifically, we designed a calibration apparatus that consists of a platform holding the ASP array driven by a custom microcontroller PCB and an arm with a light emitting diode (LED) to imitate the sun. The platform rotates to mimic a change in yaw, while the arm rotates to mimic different AOI solar light at different times of the day (Fig. 4a). Using this apparatus, we measured light input in 0.9° increments across the entire hemisphere and record the ASP array response (Fig. 4a inset) for a total resolution of 40,000 measurements in a single sweep with 200 AOI angles and 200 yaw angles. Measurements sampled by the microcontroller were transmitted to a desktop computer for logging and processing via a custom MATLAB33 script. Each ASP array response consists of a 48-bit sequence. To interpret the output, we created a lookup table from the unique 48-bit sequence that is stored for each yaw-AOI angle pair. Future data was then compared to these stored sequences in parallel using an XOR operation and the pair with the least difference in bit values was returned.
Figure 4

(a) Calibration apparatus with inset example of a measured ASP quadrature response. (b) ASP array repeatability over a single sensor (left) and multiple sensors (right). The magenta lines denote the expected operating region. (c) Expected system operating region (45°–75°) shown in magenta given the peak foraging hours (10 a.m.–4 p.m.) shown in grey and the AOI solar light during the Summer in Ithaca, NY. Plots in (a-inset), (b) and (c) were created in MATLAB33.

Full size image

We characterized sensor repeatability by repeating a sweep three times with a single ASP array and, similarly, characterized precision by comparing sweeps from two additional ASP arrays. A sample curve from the calibration sweep seen in the inset in Fig. 4a shows the characteristic angle dependence. The sensor exhibits poor response uniqueness when the sun approaches zenith and when it nears the horizon; upwards of 100(^circ) error in yaw near 90(^circ) AOI (zenith) and up to 180(^circ) error in yaw at 0°–25° AOI (horizon) (Fig. 4b). The former occurs because the position of a light source directly overhead is ambiguous to the sensor across yaw and consequently, indeterminate. We find that the sensor simply does not operate well in the latter region where light is arriving nearly parallel to the surface of the chip. Furthermore, as is expected, the difference in response is generally greater when comparing different sensors. The horizontal dark bars in the right graph in Fig. 4b are examples of this increased error. We expect our system to operate under favorable foraging conditions. Based on previously published data, we estimate this operating region to be May through September at an example location of the authors’ hometown of Ithaca, New York, USA, with the most active foraging hours being from 10 a.m. to 4 p.m.34. During this time, the AOI spans (45^circ) to (75^circ) (Fig. 4c). Within this region, we see a significantly reduced same-chip error in yaw with a mean and standard deviation of (1.52^circ pm 1.23^circ). We compensate for the remaining error within the operating region as discussed in the following sections.
In order to realistically simulate sensor output, we create a lookup table with an error model for each individual AOI value. Similar to our theoretical model, we fit a normal distribution to error in the yaw angle measured at each AOI. We use this error model to inform reconstruction of recorded bee paths as described in the following sections. For this work, we assume that we have access to calibration data for each particular sensor, however, given the low discrepancy between sensors (Fig. 4b right), we believe that it is possible to avoid individual calibration with more sophisticated data processing. We leave this aspect for future work.
Honey bee foraging simulation
To properly develop our methodology for using instrumented bees to monitor the state of pollination and bloom, we designed a colony foraging simulator with an example apple orchard. Central to our approach is an understanding of the behavior and environmental conditions surrounding honey bee foraging, summarized in Fig. 1c. The following subsections detail orchard, honey bee motion, and colony foraging models.
Orchard model
We modelled the orchard based on common characteristics seen in real orchards (Fig. 5a) as well as those reported by the University of Vermont Cooperative Extension for Growing Fruit Trees35. Specifically, these include a tree trunk radius of 0.15 m, separation between individual trees in a given planted row as 2.4 m, and separation between rows as 5 m. To make the model realistic to a variety of orchards, we add randomness to the trunk radius (0.15–0.30 m) and to the tree locations (up to 0.5 m in any direction). We further use a 60 × 60 m2 area with 200 trees, as is representative of the common grower practice utilizing a single colony per acre36. We account for the fact that trees can be in different stages of bloom by assigning each a randomly generated quality factor between 1 and 10; this quality factor affects the number of feeding events in a flight.
Colony foraging model
A high-quality honey bee colony for commercial apple pollination contains a laying queen, developing brood, and 20,000–40,000 worker bees, of which approximately 25% are “foragers”, or those that leave the hive to collect pollen, nectar, resin, and water29,37,38. Since resin and water foragers are a small proportion of the forager workforce, we expect our flight dataset to be largely from bees that are visiting flowers39. For this work, we assume favorable foraging conditions as previously described and a colony size of 35,000, 25% of which are foragers for a total of 8750 foragers conducting ~ 36,000 flights per day (an average of 4 flights per forager)37,40. In our model, a forager can perform either a learning-, return-, or scout flight. Note that we exclude orientation flights, which are conducted by new foragers, under the assumption that these can be easily classified given their tortuous nature23. A learning flight is when a bee orients to a feeding site it has not previously visited after learning the bearing and distance from one of its sisters in the hive41,42. A return flight occurs when a bee orients to a feeding site it has previously visited, and can be thought of as an optimized version of the learning flight in terms of distance flown43. A scout flight occurs when a forager leaves the hive to search independently for new feeding sites. In our model, we make the assumption based on published behavioral research that 20% of foragers are acting as scout foragers and the remaining 80% perform an initial learning flight followed by return flights to the same source41,42,43.
We incorporate that return flights will frequent the same feeding sites and that neighboring trees are likely to bloom together by randomly assigning initial goal locations (trees that bees advertise in the colony as high quality food sources) to 5 neighboring trees. Bees will randomly choose between these 5, then continue feeding on neighboring trees until they have visited trees with quality factors accumulating to at least 10 before returning to the hive. Scout foragers randomly visit trees in the orchard. Realistically, not all bees will be tagged and some tags will be lost. Here, we consider a conservative estimate that at least 430 or 5% of all foragers will be tagged, leading us to 1750 recorded flights per day, and use accumulated data to overcome the loss of tagged bees, which we expect to occur as a result of predation, senescence, stress, and other factors. Note that honey bees have a pronounced division of labor associated with worker age41, making it easy to tag a cohort and wait for them to become foragers, or to identify foragers and tag them specifically. Tagging 430 bees would take our honey bee technician approximately a day; speeding up this process is an area of future investigation.
Honey bee motion model
To better illustrate the characteristics of foraging flights, we recorded activity between a queenright colony with about 10,000 workers and a nearby feeder station (Fig. 5b–d). Three distinct phases of the bee flights were recorded: an initial orientation flight upon leaving the hive, flights between the hive and the feeder, and search flights near the feeder. Flights near the entrance and the feeder were characterized by rapid turning, whereas flights in between the hive and the feeder were nearly straight “bee lines”. While at the feeder, bees crawl around at a significantly reduced velocity.
Figure 5

(a) Photo of a honey bee in a conventional apple orchard. (b–d) Setup to showcase different types of honey bee flights. Still images from videos recorded at the hive entrance, in between, and at the feeder station, with tracked paths overlaid. (e–f) Recorded heading over the course of a simulated foraging flight and feeding event. Straight line flights are marked in grey, turns in blue, and feeding events in green. Image overlays in (b), (c) and (d) were created in MATLAB33. Plots in (e) and (f) were created in Python 3.7.

Full size image

We use this study to inform our foraging simulation. We assume generally straight paths in obstacle-free environments, slow turns when avoiding obstacles, and rapid turns in AOI and yaw when nearing and crawling on a food source. We further base our flight model on the following assumptions, summarized in Fig. 1c. (1) We assume the starting location is well known since the flight path will always originate and terminate at the hive entrance. (2) Based on past studies and the fact that our simulation takes place in a dense apple orchard, we assume most flights will be within the 4 km range of our flight recorder44,45,46. When leaving the hive, bees will fly an average of ~ 7.5 m/s, but once loaded with nectar, flight speed is reduced to ~ 6.5 m/s32. Here, we assume a constant velocity of ~ 6.5 m/s. (3) Based on prior honey bee tracking studies23,47, we represent flight in only two dimensions. The apple trees in the orchards we model are not tall and bees will therefore experience much greater motion in the horizontal plane than the vertical. Regardless of whether a bee in reality will fly over or under the canopy, we can model this issue in two dimensions. Turns around tree trunks represent the largest source of error for flight reconstruction, therefore by modeling flight under the canopy, we model the “worst case” scenario. (4) We estimate that the AOI of sunlight with respect to the orchard will remain within a quantifiable margin throughout the duration of the simulated flights, as bees have been found to spend an average of 20–45 min on foraging flights48. (5) We model the yaw of a bee as constant during bee-line flights, i.e. given no nearby obstacles. When the bee changes its heading to circumvent obstacles, this causes a change in yaw. (6) Once implemented on bees, we expect to be able to add sensors near the hive which would help us acquire current temperature and weather patterns as well as other dynamic factors specific to a particular environment for calibrating our model.
To simulate scouting, learning, and return foraging flights, we combine the honey bee motion model previously described, with the Bug2 algorithm and grid-based path planning49. Grid based path planning uses a discrete grid of points over which an agent searches to find obstacle-free path segments. We compute scout and learning flights as follows. Using the Bug2 algorithm, honey bee paths are generated by first assuming direct flight along a known heading from the hive. Once an obstacle is encountered, the bee searches for a path around it, until it can once more move unhindered toward the goal. The process is repeated until all goals are reached and the bee has returned to the hive. Since no two paths are identical in nature, we plan obstacle navigation with a randomly generated grid. Return flights are found by forming a graph of all the points visited during a learning flight and using a Dijkstra’s search algorithm49 to find the shortest path through these points from the hive to the goal. Once paths are generated, we compute the sequence of headings given the 240 ms sensor sampling frequency reported earlier. An example flight is shown in Fig. 5e,f. These headings are then discretized based on the ASP calibration data discussed earlier, and noise is added given the error shown in Fig. 4c left.
When bees land at a feeding site, they tend to crawl on and among flowers to gather nectar and pollen. We simulate this by generating random motion centered around the feeding site. The average feeding time was reported to be 1–2 min per feeding site6. To make the simulation more realistic, we randomly generate a feeding time between 60 and 120 s for foragers, and between 20 and 130 s for scouts. We furthermore assume that the AOI of sunlight changes as the honey bee tilts up and down while crawling on flowers.
Path reconstruction and generation of foraging activity maps
Path reconstruction inherently depends on the accuracy with which our sensor is able to describe the motion of an instrumented bee. Beyond limited angle resolution, errors related to the sampling rate accumulate when turns occur, at worst (v_{loaded} Delta t) = 1.6 m. To increase the accuracy of our foraging activity maps, we use models of sensor noise and flight speed, and leverage all recorded flights. The full workflow is shown in Fig. 6, where the foraging simulation portion generates the data we expect if our sensors are placed on actual bees, and the remaining flow is the data processing portion of our methodology.
Figure 6

Flow chart combining the colony foraging simulation, simulated flight recorder, path reconstruction, and data processing to output the final foraging activity map. Recorded heading plot made in Python 3.7. Other plots were created in MATLAB33.

Full size image

We first identify feeding and turn features in our path data that stand out above the noise floor. Feeding features consist of crawling behavior in which the bee is moving at much lower speed compared to flying, but with rapidly varying yaw and AOI, inducing an elevated rate of change in measured AOI. A turn is marked by a significant change in yaw. Given the time of day and location, we can find the expected AOI on the orchard and use this to find the yaw from our lookup table. Our algorithm then indexes the sensor noise lookup table to find the related mean and standard deviation and uses this to estimate turns and feeding status. Our method marks a turn for a change greater than three standard deviations in yaw ((6.06^circ)). A feeding site is marked if the detected AOI deviates more than three standard deviations ((3.72^circ)) from the one expected, or if two consecutive AOI samples deviate more than 3 standard deviations, adjustable depending on the total flight time. The average detection accuracy of a turn is 99% with a standard deviation of 0.28% and average detection accuracy of a feeding site is 99% with a standard deviation of 0.29%.
We explore three methods to generate activity maps: from the raw data, we classify feeding features using the aforementioned statistical approach and produce maps based on accumulated path reconstructions; in “iteration 1” and “iteration 4” we take a particle filter approach to improve path localization based on knowledge of the hive and frequented feeding sites respectively. Particle filters are used to track a variable of interest over time by creating many representative particles, generating predictions according to dynamics and error models, and then updating them according to observation models50. In this case, each particle forms a candidate trajectory and predictions are based on speed, sensor readings, and the sensor error model. Assuming that we start without knowledge of the orchard, we build up an observation model by reconstructing and accumulating the raw flight paths (Fig. 6, raw data). To account for the fact that bees may turn at any point between sensor readings, we then upsample our sensor readings by a factor of 8, essentially producing 8 guesses for where the bee actually turned. Based on the artificially upsampled data and a random sample from the sensor noise distribution at a given AOI, we then generate the displacement in each path segment as follows:

$$begin{aligned} begin{bmatrix} Delta x_{t} \ Delta y_{t} end{bmatrix} = v Delta t begin{bmatrix} cos {(gamma _{t} + gamma _{noise})} \ sin {(gamma _{t} + gamma _{noise})} end{bmatrix} end{aligned}$$
(6)

The final path is found as a cumulative sum of these displacements. We repeat this process to generate 5000 particles for each flight. The choice of 5000 is guided by our variable dimensions; the upsampling rate of 8 was the highest we could handle on a quadcore desktop computer with 16 GB RAM—to truly represent all potential turns we would need a number of particles equal to 8 to the power of the number of turns per flight. For reference, the average number of turns per flight is 13, thus the true representation in our sampling approach would require (8^{13}) particles.
In iteration 1, we choose ten of these particles according to proximity to the hive upon return, and use the feeding features from these to form an initial foraging activity map, represented by a discrete grid with computed visit numbers. Note that in a typical localization approach a single particle, or average of several particles, is chosen as the final reconstruction. In our approach, we retain 10 different reconstructions for each individual path in order to better represent the distribution of points in a path due to sensor noise.
After this first pass, we repeat the process, but now filter particles by using the initial activity map as an observation model for the particle filter (Fig. 6, iteration 2–4). Specifically, we update particle probability at each detected feeding site by assigning the probability of the nearest grid cell in the map to the particle, and then sampling the particles by weight. At the end of the process, the ten particles with the highest probability product are used to construct an updated activity map. More

Chemistry
General
Racemic mixture of cis/trans (35%:65%) isomers of phytol (1) (PYT) (97% purity), N-bromosuccinimide (NBS, 99% purity) and N-chlorosuccinimide (NCS, 98% purity) were purchased from Sigma-Aldrich Chemical Co. (St. Louis, MO, USA), while trimethylortoacetate was purchased from Fluka. Analytical grade acetic acid, sodium hydrogen carbonate, acetone, hexane, diethyl ether, tetrahydrofuran (THF), anhydrous magnesium sulfate, sodium chloride were purchased from Chempur (Poland).
Analytical Thin Layer Chromatography (TLC) was carried out on silica gel coated aluminium plates (DC-Alufolien Kieselgel 60 F254, Merck, Darmstadt, Germany) with a mixture of hexane, acetone and diethyl ether in various ratios as the developing systems. Compounds were visualized by spraying the plates with solution of 1% Ce(SO4)2 and 2% H3[P(Mo3O10)4] (2 g) in 10% H2SO4, followed by heating to 120–200 °C.
The products of chemical synthesis were purified by column chromatography on silica gel (Kieselgel 60, 230–400 mesh ASTM, 40–63 μm, Merck) using a mixture of hexane, acetone, and diethyl ether (in various ratios) as eluents.
Gas chromatography (GC) analysis was carried out on an Agilent Technologies 6890 N Network GC instrument (Santa Clara, CA, USA) equipped with autosampler, split injection (20:1) and FID detector using a DB-5HT column (Agilent, Santa Clara, USA) (polyimide-coated fused silica tubing, 30 m × 0.25 mm × 0.1 µm) with hydrogen as the carrier gas. Products of the chemical reactions were analysed using the following temperature programme: injector 250 °C, detector (FID) 250 °C, initial column temperature: 100 °C, 100–300 °C (rate 30 °C/min), final column temperature 300 °C (hold 2 min).
Nuclear magnetic resonance spectra 1H NMR, 13C NMR, DEPT 135, HSQC, 1H–1H COSY and NOESY were recorded in CDCl3 solutions with signals of residual solvent (δH = 7.26 δC = 77) on a Brüker Avance II 600 MHz (Brüker, Rheinstetten, Germany) spectrometer.
High-resolution mass spectra (HRMS) were recorded using electron spray ionization (ESI) technique on spectrometer Waters ESI-Q-TOF Premier XE (Waters Corp., Milford, MA, USA).
General procedure for the synthesis of compounds (2–7)
The preparation of ester 2 and acid 3 has been illustrated in detail in our previous work39, and so the synthesis method would not be listed here.
To a solution of acid 3 (7.8 mmol) in THF (30 mL) the N-bromosuccinimide (7.8 mmol) or N-chlorosuccinimide (7.8 mmol) was added. The mixture was stirred at room temperature for 48–96 h. When the substrate reacted completely (TLC, GC) the mixture was diluted with diethyl ether and washed with saturated NaHCO3 solution and brine. Organic layer of ether extract was separated and dired over anhydrous magnesium sulfate and evaporated on a rotary evaporator. New δ-halogeno-γ-lactones (4–7) were separated by silica gel column using for elution hexane/diethyl eter in gradient system. Bromo- and chlorolactonization afforded products with the following physical and spectral data presented below:
trans-5-Bromomethyl-4-methyl-4-(4′,8′,12′-trimethyltridecyl)dihydrofuran-2-one ( 4 )
(25% yield); 1H NMR (600 MHz, CDCl3): δ 0.85 (four t, J = 6.4 Hz, 12H, CH3-4′, CH3-8′, (CH3)2–12′), 1.05–1.55 (m, 21H, CH2-1′, CH2-2′, CH2-3′, CH2-5′, CH2-6′, CH2-7′, CH2-9′, CH2-10′, CH2-11′, H-4′, H-8′, H-12′), 1.24 (s, 3H, CH3-4), 2.30 and 2.60 (two d, J = 17.2 Hz, 2H, CH2-3), 3.47 (dd, J = 11.3, 7.3 Hz, 1H, one of CH2-Br), 3.55 (dd, J = 11.3, 4.5 Hz, 1H, one of CH2-Br), 4.39 (dd, J = 7.2, 4.5 Hz, 1H, H-5); 13C NMR (150 MHz, CDCl3): δ 19.61, 19.69 (CH3)2–12′), 22.65, 22.75 (CH3-4′, CH3-8′), 24.50 (CH3-4), 29.21 (CH2-Br), 41.49 (CH2-3), 42.86 (C-4), 22.00, 24.47, 24.82, 34.08, 37.28, 37.41, 37.61, 37.70, 39.38 (CH2-1′, CH2-2′, CH2-3′, CH2-5′, CH2-6′, CH2-7′, CH2-9′, CH2-10′, CH2-11′), 28.00, 32.69, 32.80 (H-4′, H-8′, H-12′), 87.66 (H-5), 174.92 (C-2); HRMS (ESI): m/z calcd. for C22H41BrO2 [M + Na]+ 439.2188; found 439.2182.
cis-5-Bromomethyl-4-methyl-4-(4′,8′,12′-trimethyltridecyl)dihydrofuran-2-one ( 5 )
(46% yield); 1H NMR (600 MHz, CDCl3): δ 0.85 (four t, J = 6.4 Hz, 12H, CH3-4′, CH3-8′, (CH3)2–12′), 1.04–1.55 (m, 21H, CH2-1′, CH2-2′, CH2-3′, CH2-5′, CH2-6′, CH2-7′, CH2-9′, CH2-10′, CH2-11′, H-4′, H-8′, H-12′), 1.08 (s, 3H, CH3-4), 2.38 and 2.49 (two d, J = 17.2 Hz, 2H, CH2-3), 3.48 (m, 2H, one CH2-Br), 4.41 (dd, J = 7.5, 4.5 Hz, 1H, H-5); 13C NMR (150 MHz, CDCl3): δ 18.96 (CH3-4), 19.69, 19.76 (CH3)2–12′), 22.65, 22.75 (CH3-4′, CH3-8′), 29.17 (CH2-Br), 42.70 (CH2-3), 43.09 (C-4), 22.20, 24.46, 24.82, 37.28, 37.39, 37.41, 37.49, 39.38, 39.96 (CH2-1′, CH2-2′, CH2-3′, CH2-5′, CH2-6′, CH2-7′, CH2-9′, CH2-10′, CH2-11′), 28.00, 30.95, 32.72 (H-4′, H-8′, H-12′), 86.46 (H-5), 174.76 (C-2); HRMS (ESI): m/z calcd. for C22H41BrO2 [M + Na]+ 439.2188; found 439.2183.
trans-5-Chloromethyl-4-methyl-4-(4′,8′,12′-trimethyltridecyl)dihydrofuran-2-one ( 6 )
(21% yield); 1H NMR (600 MHz, CDCl3): δ 0.84 (four t, J = 6.4 Hz, 12H, CH3-4′, CH3-8′, (CH3)2–12′), 1.03–1.54 (m, 21H, CH2-1′, CH2-2′, CH2-3′, CH2-5′, CH2-6′, CH2-7′, CH2-9′, CH2-10′, CH2-11′, H-4′, H-8′, H-12′), 1.23 (s, 3H, CH3-4), 2.23 and 2.60 (two d, J = 17.2 Hz, 2H, CH2-3), 3.67 (dd, J = 12.1, 6.1 Hz, 1H, one of CH2-Cl), 3.73 (dd, J = 12.1, 4.6 Hz, 1H, one of CH2-Cl), 4.33 (dd, J = 7.2, 4.6 Hz, 1H, H-5); 13C NMR (150 MHz, CDCl3): δ 19.61, 19.67 (CH3)2–12′), 22.65, 22.75 (CH3-4′, CH3-8′), 24.67 (CH3-4), 42.30 (CH2-Cl), 41.48 (CH2-3), 42.38 (C-4), 22.09, 24.46, 24.82, 34.21, 37.29, 37.39, 37.61, 37.70, 39.38 (CH2-1′, CH2-2′, CH2-3′, CH2-5′, CH2-6′, CH2-7′, CH2-9′, CH2-10′, CH2-11′), 28.00, 32.70, 32.80 (H-4′, H-8′, H-12′), 87.45 (H-5), 175.21 (C-2); HRMS (ESI): m/z calcd. for C22H41ClO2 [M + Na]+ 395.2693; found 395.2698.
cis-5-chloromethyl-4-methyl-4-(4′,8′,12′-trimethyltridecyl)dihydrofuran-2-one ( 7 )
(39% yield); 1H NMR (600 MHz, CDCl3): δ 0.85 (four t, J = 6.6 Hz, 12H, CH3-4′, CH3-8′, (CH3)2–12′), 1.04–1.56 (m, 21H, CH2-1′, CH2-2′, CH2-3′, CH2-5′, CH2-6′, CH2-7′, CH2-9′, CH2-10′, CH2-11′, H-4′, H-8′, H-12′), 1.06 (s, 3H, CH3-4), 2.38 and 2.47 (two d, J = 17.2 Hz, 2H, CH2-3), 3.67 (m, 2H, one CH2-Cl), 4.35 (dd, J = 6.4, 4.9 Hz, 1H, H-5); 13C NMR (150 MHz, CDCl3): δ 19.06 (CH3-4), 19.69, 19.76 (CH3)2–12′), 22.65, 22.74 (CH3-4′, CH3-8′), 42.52 (CH2-Cl), 42.39 (CH2-3), 42.54 (C-4), 22.14, 24.46, 24.83, 37.28, 37.37, 37.41, 37.49, 39.38, 40.16 (CH2-1′, CH2-2′, CH2-3′, CH2-5′, CH2-6′, CH2-7′, CH2-9′, CH2-10′, CH2-11′), 28.00, 32.64, 32.80 (H-4′, H-8′, H-12′), 86.24 (H-5), 175.04 (C-2); HRMS (ESI): m/z calcd. for C22H41ClO2 [M + Na]+ 395.2693; found 395.2697.
Deterrent activity of phytol and its derivatives
Aphids, plants and compound application
The peach potato aphids Myzus persicae (Sulzer) and the Chinese cabbage Brassica rapa subsp. pekinensis (Lour.) Hanelt were reared in laboratory at 20 °C, 65% r.h., and L16:8D photoperiod. One to seven days old apterous females of M. persicae and 3-week old plants with 4–5 fully developed leaves were used for experiments. M. persicae were obtained from the laboratory culture maintained at the Department of Botany and Ecology for many generations since 2000. All experiments were carried out under the same conditions of temperature, relative humidity, and photoperiod. The bioassays were started at 10–11.a.m. Each compound was dissolved in 70% ethanol to obtain the recommended 0.1% solution40. All compounds were applied on the adaxial and abaxial leaf surfaces by immersing a leaf in the ethanolic solution of a given compound for 30 s.20. Control leaves of similar size were immersed in 70% ethanol that was used as a solvent for the studied compounds. Experiments were performed 1 h after the compounds application to allow the evaporation of the solvent. Every plant and aphid were used only once.
Aphid settling (choice test)
This bioassay allows the study of aphid host preferences under semi-natural conditions41. In the present study, aphids were given free choice between control and treated excised leaves that were placed in a Petri dish. Aphids were placed in the dish equidistance from treated and untreated leaves, so that aphids could choose between treated (on one half of a Petri dish) and control leaves (on the other half of the dish). Aphids that settled, i.e. they did not move, and the position of their antennae indicated feeding, on each leaf were counted at 1 h, 2 h, and 24 h intervals after access to the leaf. Each experiment was replicated 8 times (n = 8 replicates, 20 viviparous apterous females/replicate). Aphids that were moving or not on any of the leaves were not counted.
Behavioral responses of aphids Myzus persicae during probing and feeding (no-choice test)
Aphid probing and the phloem sap uptake by M. persicae was monitored using the technique of electronic registration of aphid probing in plant tissues, known as EPG (= Electrical Penetration Graph), that is frequently employed in insect–plant relationship studies considering insects with sucking-piercing mouthparts42,43,44. In this experimental set-up, aphid and plant are connected to electrodes and thus made parts of an electric circuit, which is completed when the aphid inserts its stylets into the plant. Weak voltage is supplied in the circuit, and all changing electric properties are recorded as EPG waveforms that can be correlated with aphid activities and stylet position in plant tissues45,46. The parameters describing aphid behaviour during probing and feeding, such as total time of probing, proportion of phloem patterns E1 and E2, number of probes, etc., are good indicators of plant suitability or interference of probing by chemical or physical factors in individual plant tissues44,45,46. In the present study, aphids were attached to a golden wire electrode with conductive silver paint (epgsystems. eu) and starved for 1 h prior to the experiment. Probing behaviour of 12 apterous females/studied compound and control was monitored for 8 h continuously with Giga-4 and Giga-8 DC EPG with 1 GΩ of input resistance recording equipment (EPG Systems, Wageningen, The Netherlands). Each aphid was given access to a freshly prepared plant and each aphid/plant combination was used only once. Various behavioural phases were labelled manually using the Stylet + software (www.epgsystems.eu). The following aphid behaviours were distinguished: no penetration (waveform ‘np’ – aphid stylets outside the plant), pathway phase—penetration of non-phloem tissues (waveforms ‘ABC’), phloem phase (salivation into sieve elements, waveform ‘E1’ and ingestion of phloem sap, waveform ‘E2’), and xylem phase (ingestion of xylem sap, waveform ‘G’). Waveform ‘G’ occurred rarely irrespective of the treatment. Therefore, in all calculations the xylem phase was added to the pathway phase and termed as probing in non-phloem tissues. The E1/E2 transition patterns were included in E2. Waveform patterns that were not terminated before the end of the experimental period (8 h) were not excluded from the calculations. The parameters derived from EPG recordings were analyzed according to their frequency and duration in configuration related to activities in peripheral and vascular tissues.
Statistical analysis
The data of the choice-test were analyzed using Student’s t-test (STATISTICA 13.1. package). If aphids showed clear preference for the leaf treated with the tested compound (p  More

1.
Leroux, S. J. & Loreau, M. Subsidy hypothesis and strength of trophic cascades across ecosystems. Ecol. Lett. 11, 1147–1156 (2008).
PubMed  Article  PubMed Central  Google Scholar 
2.
Moore, J. C. et al. Detritus, trophic dynamics and biodiversity. Ecol. Lett. 7, 584–600 (2004).
ADS  Article  Google Scholar 

3.
Nowlin, W. H., Vanni, M. J. & Yang, L. H. Comparing resource pulses in aquatic and terrestrial ecosystems. Ecology 89, 647–659 (2008).
PubMed  Article  PubMed Central  Google Scholar 

4.
Wilson, E. E. & Wolkovich, E. M. Scavenging: how carnivores and carrion structure communities. Trends Ecol. Evol. 26, 129–135 (2011).
PubMed  Article  PubMed Central  Google Scholar 

5.
Margalida, A., Donázar, J. A., Carrete, M. & Sánchez-Zapata, J. A. Sanitary versus environmental policies: fitting together two pieces of the puzzle of European vulture conservation. J. Appl. Ecol. 47, 931–935 (2010).
Article  Google Scholar 

6.
Margalida, A., Colomer, M. À. & Oro, D. Man-induced activities modify demographic parameters in a long-lived species: effects of poisoning and health policies. Ecol. Appl. 24, 436–444 (2014).
PubMed  Article  PubMed Central  Google Scholar 

7.
Moreno-Opo, R. & Margalida, A. Carcasses provide resources not exclusively to scavengers: patterns of carrion exploitation by passerine birds. Ecosphere 4, art105 (2013).

8.
DeVault, T. L., Rhodes, O. E. Jr. & Shivik, J. A. Scavenging by vertebrates: behavioral, ecological, and evolutionary perspectives on an important energy transfer pathway in terrestrial ecosystems. Oikos 102, 225–234 (2003).
Article  Google Scholar 

9.
Barton, P. S., Cunningham, S. A., Lindenmayer, D. B. & Manning, A. D. The role of carrion in maintaining biodiversity and ecological processes in terrestrial ecosystems. Oecologia 171, 761–772 (2013).
ADS  PubMed  Article  PubMed Central  Google Scholar 

10.
Bump, J. K. et al. Ungulate carcasses perforate ecological filters and create biogeochemical hotspots in forest herbaceous layers allowing trees a competitive advantage. Ecosystems 12, 996–1007 (2009).
Article  Google Scholar 

11.
Danell, K., Berteaux, D. & Bråthen, K. A. Effect of muskox carcasses on nitrogen concentration in tundra vegetation. Arctic 55, 389–392 (2002).
Article  Google Scholar 

12.
Klink, R., Laar-Wiersma, J., Vorst, O. & Smit, C. Rewilding with large herbivores: positive direct and delayed effects of carrion on plant and arthropod communities. PLoS ONE 15, e0226946 (2020).
PubMed  PubMed Central  Article  CAS  Google Scholar 

13.
Turner, W. C. et al. Fatal attraction: vegetation responses to nutrient inputs attract herbivores to infectious anthrax carcass sites. Proc. R. Soc. Lond. B Biol. Sci. 281, e20141785 (2014).

14.
Mateo-Tomás, P. et al. From regional to global patterns in vertebrate scavenger communities subsidized by big game hunting. Divers. Distrib. 21, 913–924 (2015).
Article  Google Scholar 

15.
Markandya, A. et al. Counting the cost of vulture decline—an appraisal of the human health and other benefits of vultures in India. Ecol. Econ. 67, 194–204 (2008).
Article  Google Scholar 

16.
Selva, N., Jędrzejewska, B., Jędrzejewski, W. & Wajrak, A. Factors affecting carcass use by a guild of scavengers in European temperate woodland. Can. J. Zool. 83, 1590–1601 (2005).
Article  Google Scholar 

17.
DeVault, T. L., Brisbin, J., Lehr, I., Rhodes, J. & Olin, E. Factors influencing the acquisition of rodent carrion by vertebrate scavengers and decomposers. Can. J. Zool. 82, 502–509 (2004).
Article  Google Scholar 

18.
Arrondo, E. et al. Rewilding traditional grazing areas affects scavenger assemblages and carcass consumption patterns. Basic Appl. Ecol. 41, 56–66 (2019).
Article  Google Scholar 

19.
Morales-Reyes, Z. et al. Scavenging efficiency and red fox abundance in Mediterranean mountains with and without vultures. Acta Oecologica 79, 81–88 (2017).
ADS  Article  Google Scholar 

20.
Ruzicka, R. E. & Conover, M. R. Does weather or site characteristics influence the ability of scavengers to locate food? Ethology 118, 187–196 (2012).
Article  Google Scholar 

21.
Moleón, M., Sánchez-Zapata, J., Sebastián-González, E. & Owen-Smith, N. Carcass size shapes the structure and functioning of an African scavenging assemblage. Oikos 124, 1391–1403 (2015).

22.
Cornwell, W. K. et al. Plant species traits are the predominant control on litter decomposition rates within biomes worldwide. Ecol. Lett. 11, 1065–1071 (2008).
PubMed  Article  PubMed Central  Google Scholar 

23.
Ogada, D. L., Torchin, M. E., Kinnaird, M. F. & Ezenwa, V. O. Effects of vulture declines on facultative scavengers and potential implications for mammalian disease transmission. Conserv. Biol. 26, 453–460 (2012).
CAS  PubMed  Article  PubMed Central  Google Scholar 

24.
Sekercioglu, Ç. H., Wenny, D. G. & Whelan, C. J. Why Birds Matter: Avian Ecological Function and Ecosystem Services (University of Chicago Press, 2016).

25.
Pereira, L. M., Owen-Smith, N. & Moleón, M. Facultative predation and scavenging by mammalian carnivores: seasonal, regional and intra-guild comparisons. Mammal Rev. 44, 44–55 (2014).
Article  Google Scholar 

26.
Selva, N. & Fortuna, M. A. The nested structure of a scavenger community. Proc. R. Soc. B Biol. Sci. 274, 1101–1108 (2007).
Article  Google Scholar 

27.
Wolf, C. & Ripple, W. J. Range contractions of the world’s large carnivores. R. Soc. Open Sci. 4, 170052 (2017).
ADS  PubMed  PubMed Central  Article  Google Scholar 

28.
Grimm, N. B. et al. The impacts of climate change on ecosystem structure and function. Front. Ecol. Environ. 11, 474–482 (2013).
Article  Google Scholar 

29.
Lauenroth, W. et al. Potential effects of climate change on the temperate zones of North and South America. Rev. Chil. Hist. Nat. 77, 439–453 (2004).
Article  Google Scholar 

30.
Shanley, C. S. et al. Climate change implications in the northern coastal temperate rainforest of North America. Clim. Change 130, 155–170 (2015).
ADS  CAS  Article  Google Scholar 

31.
Wilmers, C. C. & Getz, W. M. Gray wolves as climate change buffers in yellowstone. PLOS Biol. 3, e92 (2005).
PubMed  PubMed Central  Article  CAS  Google Scholar 

32.
Sebastián-González, E. et al. Network structure of vertebrate scavenger assemblages at the global scale: drivers and ecosystem functioning implications. Ecography 43, 1143–1155 (2020).
Article  Google Scholar 

33.
Pardo-Barquín, E., Mateo-Tomás, P. & Olea, P. P. Habitat characteristics from local to landscape scales combine to shape vertebrate scavenging communities. Basic Appl. Ecol. 34, 126–139 (2019).
Article  Google Scholar 

34.
Sebastián-González, E. et al. Scavenging in the Anthropocene: human impact drives vertebrate scavenger species richness at a global scale. Glob. Change Biol. 25, 3005–3017 (2019).
ADS  Article  Google Scholar 

35.
Turner, K. L., Abernethy, E. F., Conner, L. M., Rhodes, O. E. & Beasley, J. C. Abiotic and biotic factors modulate carrion fate and vertebrate scavenging communities. Ecology 98, 2413–2424 (2017).
PubMed  Article  PubMed Central  Google Scholar 

36.
Janßen, F., Treude, T. & Witte, U. Scavenger assemblages under differing trophic conditions: a case study in the deep Arabian Sea. Deep Sea Res. Part II Top. Stud. Oceanogr. 47, 2999–3026 (2000).
ADS  Article  Google Scholar 

37.
Houston, D. C. To the vultures belong the spoils. Nat. Hist. 103, 34–41 (1994).
Google Scholar 

38.
Houston, D. C. Scavenging efficiency of turkey vultures in tropical forest. The Condor 88, 318–323 (1986).
Article  Google Scholar 

39.
Sauer, J. et al. The North American breeding bird survey, results and analysis 1966–2015. (2017).

40.
Hill, J. E., DeVault, T. L., Beasley, J. C., Rhodes, O. E. & Belant, J. L. Effects of vulture exclusion on carrion consumption by facultative scavengers. Ecol. Evol. 8, 2518–2526 (2018).
PubMed  PubMed Central  Article  Google Scholar 

41.
Heinrich, B. Winter foraging at carcasses by three sympatric corvids, with emphasis on recruitment by the raven, Corvus corax. Behav. Ecol. Sociobiol. 23, 141–156 (1988).
Article  Google Scholar 

42.
Bellan, S. E., Turnbull, P. C. B., Beyer, W. & Getz, W. M. Effects of experimental exclusion of scavengers from carcasses of anthrax-infected herbivores on bacillus anthracis sporulation, survival, and distribution. Appl. Environ. Microbiol. 79, 3756–3761 (2013).
CAS  PubMed  PubMed Central  Article  Google Scholar 

43.
The IUCN Red List of Threatened Species. IUCN Red List of Threatened Species https://www.iucnredlist.org/en.

44.
Kiff, L. F. The current status of North American vultures. In Raptors at Risk 175–189 (World Working Group on Birds of Prey/Hancock House, 2000).

45.
Prasad, A. M., Iverson, L. R., Peters, M. P. & Matthews, S. N. Climate change tree atlas (Northern Research Station, US Forest Service, Delaware, OH, 2014).
Google Scholar 

46.
Kiff, L. The current status of North American vultures. in 175–189 (2000).

47.
Houston, D. C. Competition for food between Neotropical vultures in forest. Ibis 130, 402–417 (1988).
Article  Google Scholar 

48.
Gomez, L. G., Houston, D. C., Cotton, P. & Tye, A. The role of greater yellow-headed vultures Cathartes melambrotus as scavengers in neotropical forest. Ibis 136, 193–196 (1994).
Article  Google Scholar 

49.
Ripple, W. J. et al. Status and ecological effects of the world’s largest carnivores. Science 343, e1241484 (2014).

50.
Tomberlin, J. K., Barton, B. T., Lashley, M. A. & Jordan, H. R. Mass mortality events and the role of necrophagous invertebrates. Curr. Opin. Insect Sci. 23, 7–12 (2017).
PubMed  Article  PubMed Central  Google Scholar 

51.
Fey, S. B. et al. Recent shifts in the occurrence, cause, and magnitude of animal mass mortality events. Proc. Natl. Acad. Sci. 112, 1083–1088 (2015).
ADS  CAS  PubMed  Article  PubMed Central  Google Scholar 

52.
Wikenros, C., Sand, H., Ahlqvist, P. & Liberg, O. Biomass flow and scavengers use of carcasses after re-colonization of an apex predator. PLoS ONE 8, e77373 (2013).
ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

53.
Kočárek, P. Decomposition and Coleoptera succession on exposed carrion of small mammal in Opava, the Czech Republic. Eur. J. Soil Biol. 39, 31–45 (2003).
Article  Google Scholar 

54.
Matuszewski, S., Bajerlein, D., Konwerski, S. & Szpila, K. Insect succession and carrion decomposition in selected forests of Central Europe. Part 1: pattern and rate of decomposition. Forensic Sci. Int. 194, 85–93 (2010).
PubMed  Article  PubMed Central  Google Scholar 

55.
Reed, H. B. A study of dog carcass communities in tennessee, with special reference to the insects. Am. Midl. Nat. 59, 213–245 (1958).
Article  Google Scholar 

56.
Bauer, J. W., Logan, K. A., Sweanor, L. L. & Boyce, W. M. Scavenging behavior in Puma. Southwest. Nat. 50, 466–471 (2005).
Article  Google Scholar 

57.
Burkepile, D. E. et al. Chemically mediated competition between microbes and animals: microbes as consumers in food webs. Ecology 87, 2821–2831 (2006).
PubMed  Article  PubMed Central  Google Scholar 

58.
Janzen, D. H. Why fruits rot, seeds mold, and meat spoils. Am. Nat. 111, 691–713 (1977).
CAS  Article  Google Scholar 

59.
DeVault, T. L. & Rhodes, O. E. Identification of vertebrate scavengers of small mammal carcasses in a forested landscape. Acta Theriol. (Warsz.) 47, 185–192 (2002).
Article  Google Scholar 

60.
Parker, K. L., Robbins, C. T. & Hanley, T. A. Energy expenditures for locomotion by Mule Deer and Elk. J. Wildl. Manag. 48, 474–488 (1984).
Article  Google Scholar 

61.
Crête, M. & Larivière, S. Estimating the costs of locomotion in snow for coyotes. Can. J. Zool. 81, 1808–1814 (2003).
Article  Google Scholar 

62.
Droghini, A. & Boutin, S. The calm during the storm: snowfall events decrease the movement rates of grey wolves (Canis lupus). PLoS ONE 13, e0205742 (2018).

63.
Green, G. I., Mattson, D. J. & Peek, J. M. Spring feeding on ungulate carcasses by grizzly bears in Yellowstone National Park. J. Wildl. Manag. 61, 1040–1055 (1997).
Article  Google Scholar 

64.
De Jong, G. D. & Chadwick, J. W. Decomposition and arthropod succession on exposed rabbit carrion during summer at high altitudes in colorado, USA. J. Med. Entomol. 36, 833–845 (1999).
PubMed  Article  PubMed Central  Google Scholar 

65.
Sun, S.-J. et al. Climate-mediated cooperation promotes niche expansion in burying beetles. Elife 3, e02440 (2014).
PubMed  PubMed Central  Article  Google Scholar 

66.
Krofel, M. Monitoring of facultative avian scavengers on large mammal carcasses in Dinaric forest of Slovenia. Acrocephalus 32, 45–51 (2011).
Article  Google Scholar 

67.
DeVault, T. L., Seamans, T. W., Linnell, K. E., Sparks, D. W. & Beasley, J. C. Scavenger removal of bird carcasses at simulated wind turbines: Does carcass type matter?. Ecosphere 8, e01994 (2017).
Article  Google Scholar 

68.
Turner, K. L., Conner, L. M. & Beasley, J. C. Effect of mammalian mesopredator exclusion on vertebrate scavenging communities. Sci. Rep. 10, 2644 (2020).
ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

69.
Abernethy, E. F., Turner, K. L., Beasley, J. C. & Rhodes, O. E. Scavenging along an ecological interface: utilization of amphibian and reptile carcasses around isolated wetlands. Ecosphere 8, e01989 (2017).
Article  Google Scholar 

70.
Olson, Z. H., Beasley, J. C. & Rhodes, O. E. Carcass type affects local scavenger guilds more than habitat connectivity. PLoS ONE 11, (2016).

71.
Ragg, J., Mackintosh, C. & Moller, H. The scavenging behaviour of ferrets (Mustela furo), feral cats (Felis domesticus), possums (Trichosurus vulpecula), hedgehogs (Erinaceus europaeus) and harrier hawks (Circus approximans) on pastoral farmland in New Zealand: Implications for bovine tuberculosis transmission. N. Z. Vet. J. 48, 166–175 (2001).
Article  Google Scholar 

72.
Laundré, J. W., Hernández, L. & Altendorf, K. B. Wolves, elk, and bison: reestablishing the” landscape of fear” in Yellowstone National Park, USA. Can. J. Zool. 79, 1401–1409 (2001).
Article  Google Scholar 

73.
Ripple, W. J. & Beschta, R. L. Linking wolves to willows via risk-sensitive foraging by ungulates in the northern Yellowstone ecosystem. For. Ecol. Manag. 230, 96–106 (2006).
Article  Google Scholar 

74.
Ripple, W. J. & Beschta, R. L. Trophic cascades in Yellowstone: the first 15 years after wolf reintroduction. Biol. Conserv. 145, 205–213 (2012).
Article  Google Scholar 

75.
Smith, D. W., Peterson, R. O. & Houston, D. B. Yellowstone after wolves. Bioscience 53, 330–340 (2003).
Article  Google Scholar 

76.
White, P. J. & Garrott, R. A. Northern Yellowstone elk after wolf restoration. Wildl. Soc. Bull. 33, 942–955 (2005).
Article  Google Scholar 

77.
Clark, P. J. & Evans, F. C. Distance to nearest neighbor as a measure of spatial relationships in populations. Ecology 35, 445–453 (1954).
Article  Google Scholar 

78.
Cook, R. C., Cook, J. G. & Irwin, L. L. Estimating elk body mass using chest-girth circumference. Wildl. Soc. Bull. 1973-2006 31, 536–543 (2003).
Google Scholar 

79.
Craine, J. M., Towne, E. G. & Elmore, A. Intra-annual bison body mass trajectories in a tallgrass prairie. Mammal Res. 60, 263–270 (2015).
Article  Google Scholar 

80.
Lott, D. F. & Galland, J. C. Body mass as a factor influencing dominance status in American Bison Cows. J. Mammal. 68, 683–685 (1987).
Article  Google Scholar 

81.
Fox, J. & Weisberg, S. An R Companion to Applied Regression. (Sage Publications, 2018).

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

83.
Pan, Y. & Jackson, R. T. Ethnic difference in the relationship between acute inflammation and serum ferritin in US adult males. Epidemiol. Infect. 136, 421–431 (2008).
CAS  PubMed  Article  PubMed Central  Google Scholar 

84.
Brewer, M. J., Butler, A. & Cooksley, S. L. The relative performance of AIC, AICC and BIC in the presence of unobserved heterogeneity. Methods Ecol. Evol. 679, 692. https://doi.org/10.1111/2041-210X.12541 (2016).
Article  Google Scholar 

85.
Burnham, K. P. & Anderson, D. R. Multimodel inference: understanding AIC and BIC in model selection. Sociol. Methods Res. 33, 261–304 (2004).
MathSciNet  Article  Google Scholar 

86.
Franklin, J. Mapping Species Distributions: Spatial Inference and Prediction (Cambridge University Press, Cambridge, 2010).
Google Scholar 

87.
Kleiber, C. & Zeileis, A. Applied Econometrics with R (Springer, Berlin, 2008).
Google Scholar 

88.
Cameron, A. C. & Trivedi, P. K. Regression-based tests for overdispersion in the Poisson model. J. Econom. 46, 347–364 (1990).
MathSciNet  Article  Google Scholar  More

1.
Fischer, H. et al. Palaeoclimate constraints on the impact of 2 °C anthropogenic warming and beyond. Nat. Geosci. 11, 474–485 (2018).
ADS  CAS  Article  Google Scholar 
2.
Diffenbaugh, N. S., Singh, D. & Mankin, J. S. Unprecedented climate events: historical changes, aspirational targets, and national commitments. Sci. Adv. 4, eaao3354 (2018).
ADS  PubMed  PubMed Central  Article  CAS  Google Scholar 

3.
Kim, J.-S. et al. Reduced North American terrestrial primary productivity linked to anomalous Arctic warming. Nat. Geosci. 10, 572–576 (2017).
ADS  CAS  Article  Google Scholar 

4.
United Nations Framework Convention on Climate Change (UNFCCC). Adoption of the Paris Agreement (2015).

5.
Intergovernmental Panel on Climate Change (IPCC). Fifth Assessment Report: Climate Change (AR5) (2014).

6.
Nabuurs, G. J. et al. Forestry. In Climate Change 2007: Mitigation (Cambridge University Press, 2007).

7.
Smyth, C. E. et al. Quantifying the biophysical climate change mitigation potential of Canada’s forest sector. Biogeosciences 11, 3515–3529 (2014).
ADS  Article  Google Scholar 

8.
Xu, Z., Smyth, C. E., Lemprière, T. C., Rampley, G. J. & Kurz, W. A. Climate change mitigation strategies in the forest sector: biophysical impacts and economic implications in British Columbia. Can. Mitig. Adapt. Strateg. Glob. Change 23, 257–290 (2018).
Article  Google Scholar 

9.
Bastin, J.-F. et al. The global tree restoration potential. Science 365, 76–79 (2019).
ADS  CAS  PubMed  Article  Google Scholar 

10.
Peterson St-Laurent, G., Hagerman, S., Kozak, R. & Hoberg, G. Public perceptions about climate change mitigation in British Columbia’s forest sector. PLoS ONE 13, e0195999 (2018).
PubMed  PubMed Central  Article  CAS  Google Scholar 

11.
Ma, Z. et al. Regional drought-induced reduction in the biomass carbon sink of Canada’s boreal forests. Proc. Natl Acad. Sci. USA 109, 2423–2427 (2012).
ADS  CAS  PubMed  Article  Google Scholar 

12.
Charney, N. D. et al. Observed forest sensitivity to climate implies large changes in 21st century North American forest growth. Ecol. Lett. 19, 1119–1128 (2016).
PubMed  Article  Google Scholar 

13.
Girardin, M. P. et al. No growth stimulation of Canada’s boreal forest under half-century of combined warming and CO2 fertilization. Proc. Natl Acad. Sci. USA 113, E8406–E8414 (2016).
CAS  PubMed  Article  Google Scholar 

14.
Marchand, W. et al. Untangling methodological and scale considerations in growth and productivity trend estimates of Canada’s forests. Environ. Res. Lett. 13, 093001 (2018).
ADS  Article  Google Scholar 

15.
Browne, L., Wright, J. W., Fitz-Gibbon, S., Gugger, P. F. & Sork, V. L. Adaptational lag to temperature in valley oak (Quercus lobata) can be mitigated by genome-informed assisted gene flow. Proc. Natl Acad. Sci. USA. 116, 25179–25185 (2019).
CAS  PubMed  Article  Google Scholar 

16.
Sally, N. A. & Whitlock, M. C. Assisted gene flow to facilitate local adaptation to climate change. Annu Rev. Ecol. Evol. Syst. 44, 367–388 (2013).
Article  Google Scholar 

17.
Lemprière, T. C. et al. Canadian boreal forests and climate change mitigation. Environ. Rev. 21, 293–321 (2013).
Article  Google Scholar 

18.
Winder, R., Nelson, E. & Beardmore, T. Ecological implications for assisted migration in Canadian forests. For. Chron. 87, 731–744 (2011).
Article  Google Scholar 

19.
Teskey, R. et al. Responses of tree species to heat waves and extreme heat events. Plant Cell Environ. 38, 1699–1712 (2015).
PubMed  Article  Google Scholar 

20.
Isaac-Renton, M. et al. Northern forest tree populations are physiologically maladapted to drought. Nat. Commun. 9, 5254 (2018).
ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

21.
Field, E., Schönrogge, K., Barsoum, N., Hector, A. & Gibbs, M. Individual tree traits shape insect and disease damage on oak in a climate‐matching tree diversity experiment. Ecol. Evol. 9, 8524–8540 (2019).
PubMed  PubMed Central  Article  Google Scholar 

22.
Depardieu, C. et al. Adaptive genetic variation to drought in a widely distributed conifer suggests a potential for increasing forest resilience in a drying climate. N. Phytol. 227, 427–439 (2020).
CAS  Article  Google Scholar 

23.
Montwé, D., Isaac-Renton, M., Hamann, A. & Spiecker, H. Cold adaptation recorded in tree rings highlights risks associated with climate change and assisted migration. Nat. Commun. 9, 1574 (2018).
ADS  PubMed  PubMed Central  Article  CAS  Google Scholar 

24.
Girardin, M. P. et al. Negative impacts of high temperatures on growth of black spruce forests intensify with the anticipated climate warming. Glob. Change Biol. 22, 627–643 (2016).
ADS  Article  Google Scholar 

25.
Hember, R. A., Kurz, W. A. & Coops, N. C. Increasing net ecosystem biomass production of Canada’s boreal and temperate forests despite decline in dry climates. Glob. Biogeochem. Cycles 31, 134–158 (2017).
CAS  Article  Google Scholar 

26.
Boucher, D. et al. Current and projected cumulative impacts of fire, drought, and insects on timber volumes across Canada. Ecol. Appl. 28, 1245–1259 (2018).
PubMed  Article  Google Scholar 

27.
Klein, R. J. T., Schipper, E. L. F. & Dessai, S. Integrating mitigation and adaptation into climate and development policy: three research questions. Environ. Sci. Policy 8, 579–588 (2005).
Article  Google Scholar 

28.
Zamudio, K. R., Bell, R. C. & Mason, N. A. Phenotypes in phylogeography: species’ traits, environmental variation, and vertebrate diversification. Proc. Natl Acad. Sci. USA 113, 8041–8048 (2016).
CAS  PubMed  Article  Google Scholar 

29.
Massatti, R. et al. Population history provides foundational knowledge for utilizing and developing native plant restoration materials. Evol. Appl. 11, 2025–2039 (2018).
PubMed  PubMed Central  Article  Google Scholar 

30.
Sork, V. L. et al. Putting the landscape into the genomics of trees: approaches for understanding local adaptation and population responses to changing climate. Tree Genet. Genomes 9, 901–911 (2013).
Article  Google Scholar 

31.
Morgenstern, E. K. & Mullin, T. J. Growth and survival of black spruce in the range-wide provenance study. Can. J. Res. 20, 130–143 (1990).
Article  Google Scholar 

32.
Rehfeldt, G. E., Wykoff, W. R. & Ying, C. C. Physiologic plasticity, evolution, and impacts of a changing climate on Pinus contorta. Clim. Change 50, 355–376 (2001).
Article  Google Scholar 

33.
Morgenstern, E. K. Range-wide genetic variation of black spruce. Can. J. Res. 8, 463–473 (1978).
Article  Google Scholar 

34.
Thomson, A. M., Riddell, C. L. & Parker, W. H. Boreal forest provenance tests used to predict optimal growth and response to climate change: 2. Black spruce. Can. J. Res. 39, 143–153 (2009).
Article  Google Scholar 

35.
Pedlar, J. H. & McKenney, D. W. Assessing the anticipated growth response of northern conifer populations to a warming climate. Sci. Rep. 7, 43881 (2017).
ADS  PubMed  PubMed Central  Article  Google Scholar 

36.
Mahony, C. R. et al. Evaluating genomic data for management of local adaptation in a changing climate: a lodgepole pine case study. Evol. Appl. 13, 116–131 (2020).
PubMed  Article  Google Scholar 

37.
Housset, J. M. et al. Tree rings provide a new class of phenotypes for genetic associations that foster insights into adaptation of conifers to climate change. N. Phytol. 218, 630–645 (2018).
Article  Google Scholar 

38.
Heer, K. et al. Linking dendroecology and association genetics in natural populations: stress responses archived in tree rings associate with SNP genotypes in silver fir (Abies alba Mill.). Mol. Ecol. 27, 1428–1438 (2018).
CAS  PubMed  Article  Google Scholar 

39.
Bouriaud, O., Teodosiu, M., Kirdyanov, A. V. & Wirth, C. Influence of wood density in tree-ring-based annual productivity assessments and its errors in Norway spruce. Biogeosciences 12, 6205–6217 (2015).
ADS  CAS  Article  Google Scholar 

40.
Babst, F. et al. When tree rings go global: challenges and opportunities for retro- and prospective insight. Quat. Sci. Rev. 197, 1–20 (2018).
ADS  Article  Google Scholar 

41.
Jaramillo-Correa, J. P., Beaulieu, J. & Bousquet, J. Variation in mitochondrial DNA reveals multiple distant glacial refugia in black spruce (Picea mariana), a transcontinental North American conifer. Mol. Ecol. 13, 2735–2747 (2004).
CAS  PubMed  Article  Google Scholar 

42.
Gérardi, S., Jaramillo-Correa, J. P., Beaulieu, J. & Bousquet, J. From glacial refugia to modern populations: new assemblages of organelle genomes generated by differential cytoplasmic gene flow in transcontinental black spruce: assemblages of organelle genomes. Mol. Ecol. 19, 5265–5280 (2010).
PubMed  Article  CAS  Google Scholar 

43.
Rehfeldt, G. E., Leites, L. P., Joyce, D. G. & Weiskittel, A. R. Role of population genetics in guiding ecological responses to climate. Glob. Change Biol. 24, 858–868 (2018).
ADS  Article  Google Scholar 

44.
Beaulieu, J., Corriveau, A. & Daoust, G. Phenotypic Stability and Delineation of Black Spruce Breeding Zones in Quebec. Vol. LAU-X-85E (Forestry Canada, Quebec Region, Sainte-Foy, Quebec, 1989).

45.
Perrin, M., Rossi, S. & Isabel, N. Synchronisms between bud and cambium phenology in black spruce: early-flushing provenances exhibit early xylem formation. Tree Physiol. 37, 593–603 (2017).
PubMed  Article  Google Scholar 

46.
Sniderhan, A. E., McNickle, G. G. & Baltzer, J. L. Assessing local adaptation vs. plasticity under different resource conditions in seedlings of a dominant boreal tree species. AoB Plants 10, ply004 (2018).

47.
Newton, P. F. Systematic review of yield responses of four North American conifers to forest tree improvement practices. Ecol. Manag. 172, 29–51 (2003).
Article  Google Scholar 

48.
Marchand, W. et al. Strong overestimation of water‐use efficiency responses to rising CO2 in tree‐ring studies. Glob. Change Biol. https://doi.org/10.1111/gcb.15166 (2020).

49.
Metsaranta, J. M. Long-term tree-ring derived carbon dynamics of an experimental plantation in relation to species and density in Northwestern Ontario. Can. Ecol. Manag. 441, 229–241 (2019).
Article  Google Scholar 

50.
Büntgen, U. et al. Limited capacity of tree growth to mitigate the global greenhouse effect under predicted warming. Nat. Commun. 10, 2171 (2019).
ADS  PubMed  PubMed Central  Article  CAS  Google Scholar 

51.
DOrangeville, L. et al. Northeastern North America as a potential refugium for boreal forests in a warming climate. Science 352, 1452–1455 (2016).
ADS  CAS  Article  Google Scholar 

52.
Babst, F. et al. Twentieth century redistribution in climatic drivers of global tree growth. Sci. Adv. 5, eaat4313 (2019).
ADS  PubMed  PubMed Central  Article  Google Scholar 

53.
Rossi, S. Bud break responds more strongly to daytime than night‐time temperature under asymmetric experimental warming. Glob. Change Biol. 9 (2016).

54.
Frechette, E., Ensminger, I., Bergeron, Y., Gessler, A. & Berninger, F. Will changes in root-zone temperature in boreal spring affect recovery of photosynthesis in Picea mariana and Populus tremuloides in a future climate? Tree Physiol. 31, 1204–1216 (2011).
CAS  PubMed  Article  Google Scholar 

55.
Verbyla, D. Remote sensing of interannual boreal forest NDVI in relation to climatic conditions in interior Alaska. Environ. Res. Lett. 10, 125016 (2015).
ADS  Article  Google Scholar 

56.
Trujillo, E. Elevation-dependent influence of snow accumulation on forest greening. Nat. Geosci. 5, 5 (2012).
Article  CAS  Google Scholar 

57.
Vaganov, E. A., Hughes, M. K., Kirdyanov, A. V., Schweingruber, F. H. & Silkin, P. P. Influence of snowfall and melt timing on tree growth in subarctic Eurasia. Nature 400, 149–151 (1999).
ADS  CAS  Article  Google Scholar 

58.
Ols, C., Girardin, M. P., Hofgaard, A., Bergeron, Y. & Drobyshev, I. Monitoring climate sensitivity shifts in tree-rings of eastern boreal North America using model-data comparison: shifts in tree growth sensivity to climate. Ecosystems 21, 1042–1057 (2018).
CAS  Article  Google Scholar 

59.
Prunier, J., Gérardi, S., Laroche, J., Beaulieu, J. & Bousquet, J. Parallel and lineage-specific molecular adaptation to climate in boreal black spruce. Mol. Ecol. 21, 4270–4286 (2012).
CAS  PubMed  Article  Google Scholar 

60.
Capblancq, T. et al. Climate-associated genetic variation in Fagus sylvatica and potential responses to climate change in the French Alps. J. Evol. Biol. https://doi.org/10.1111/jeb.13610 (2020).

61.
Liepe, K. J., Hamann, A., Smets, P., Fitzpatrick, C. R. & Aitken, S. N. Adaptation of lodgepole pine and interior spruce to climate: implications for reforestation in a warming world. Evol. Appl. 9, 409–419 (2016).
PubMed  PubMed Central  Article  Google Scholar 

62.
Fitzpatrick, M. C. & Keller, S. R. Ecological genomics meets community-level modelling of biodiversity: mapping the genomic landscape of current and future environmental adaptation. Ecol. Lett. 18, 1–16 (2015).
PubMed  Article  Google Scholar 

63.
Wegrzyn, J. L. et al. Cyberinfrastructure and resources to enable an integrative approach to studying forest trees. Evol. Appl. 13, 228–241 (2020).
PubMed  Article  Google Scholar 

64.
Kurz, W. A. et al. Carbon in Canada’s boreal forest —a synthesis. Environ. Rev. 21, 260–292 (2013).
CAS  Article  Google Scholar 

65.
Alberto, F. J. et al. Potential for evolutionary responses to climate change—evidence from tree populations. Glob. Change Biol. 19, 1645–1661 (2013).
ADS  Article  Google Scholar 

66.
Splawinski, T. B., Cyr, D., Gauthier, S., Jetté, J.-P. & Bergeron, Y. Analyzing risk of regeneration failure in the managed boreal forest of northwestern Quebec. Can. J. Res. 49, 680–691 (2019).
Article  Google Scholar 

67.
Chaste, E., Girardin, M. P., Kaplan, J. O., Bergeron, Y. & Hély, C. Increases in heat-induced tree mortality could drive reductions of biomass resources in Canada’s managed boreal forest. Landsc. Ecol. https://doi.org/10.1007/s10980-019-00780-4 (2019).

68.
McLane, S. C., Daniels, L. D. & Aitken, S. N. Climate impacts on lodgepole pine (Pinus contorta) radial growth in a provenance experiment. Ecol. Manag. 262, 115–123 (2011).
Article  Google Scholar 

69.
Larsson, L. CooRecorder and Cdendro Programs of the CooRecorder/Cdendro Package (Version 7.6) (Cybis Elektronik, 2013).

70.
Holmes, R. L. Computer-assisted quality control in tree-ring dating and measurement. Tree Ring Bull. 43, 69–78 (1983).

71.
de Lafontaine, G., Prunier, J., Gérardi, S. & Bousquet, J. Tracking the progression of speciation: variable patterns of introgression across the genome provide insights on the species delimitation between progenitor-derivative spruces (Picea mariana × P. rubens). Mol. Ecol. 24, 5229–5247 (2015).
PubMed  Article  Google Scholar 

72.
Ehrich, M. et al. Quantitative high-throughput analysis of DNA methylation patterns by base-specific cleavage and mass spectrometry. Proc. Natl Acad. Sci. USA. 102, 15785–15790 (2005).
ADS  CAS  PubMed  Article  Google Scholar 

73.
Ung, C.-H., Jing Guo, X. & Fortin, M. Canadian national taper models. For. Chron. 89, 211–224 (2013).
Article  Google Scholar 

74.
Pritchard, J. K., Stephens, M. & Donnelly, P. Inference of population structure using multilocus genotype data. Genetics 155, 945–959 (2000).
CAS  PubMed  PubMed Central  Google Scholar 

75.
Wang, J. The computer program structure for assigning individuals to populations: easy to use but easier to misuse. Mol. Ecol. Resour. 17, 981–990 (2017).
CAS  PubMed  Article  Google Scholar 

76.
Evanno, G., Regnaut, S. & Goudet, J. Detecting the number of clusters of individuals using the software structure: a simulation study. Mol. Ecol. 14, 2611–2620 (2005).
CAS  PubMed  Article  Google Scholar 

77.
Excoffier, L. Evolution of human mitochondrial DNA: evidence for departure from a pure neutral model of populations at equilibrium. J. Mol. Evol. 30, 125–139 (1990).
ADS  CAS  PubMed  Article  Google Scholar 

78.
Meirmans, P. G. genodive version 3.0: easy‐to‐use software for the analysis of genetic data of diploids and polyploids. Mol. Ecol. Resour. 20, 1126–1131 (2020).
CAS  PubMed  PubMed Central  Article  Google Scholar 

79.
Regniere, J. & Bolstad, P. Statistical simulation of daily air temperature patterns Eastern North America to forecast seasonal events in insect pest management. Environ. Entomol. 23, 1368–1380 (1994).
Article  Google Scholar 

80.
Hogg, E. H., Barr, A. G. & Black, T. A. A simple soil moisture index for representing multi-year drought impacts on aspen productivity in the western Canadian interior. Agric. Meteorol. 178–179, 173–182 (2013).
Article  Google Scholar 

81.
Wood, S. N. Thin plate regression splines. J. R. Stat. Soc. Ser. B Stat. Methodol. 65, 95–114 (2003).
MathSciNet  MATH  Article  Google Scholar 

82.
Rossi, S., Morin, H. & Deslauriers, A. Causes and correlations in cambium phenology: towards an integrated framework of xylogenesis. J. Exp. Bot. 63, 2117–2126 (2012).
CAS  PubMed  Article  Google Scholar 

83.
Wood, S. Generalized Additive Models: An Introduction with R. 2nd edn Vol. 66 (Chapman and Hall/CRC, 2006).

84.
R Development Core Team. R: A Language and Environment for Statistical Computing (2013).

85.
Dray, S. & Dufour, A.-B. The ade4 package: implementing the duality diagram for ecologists. J. Stat. Softw. 22, https://doi.org/10.18637/jss.v022.i04 (2007).

86.
Fisher R. A. Statistical Methods for Research Workers 4th edn (Oliver and Boyd, London, 1932).

87.
Legendre, P. & Legendre, L. Numerical Ecology Vol. 24, 3rd edn (Elsevier Science BV, Amsterdam, 2012).

88.
Reiss, P. T. & Ogden, R. T. Smoothing parameter selection for a class of semiparametric linear models. J. R. Stat. Soc. Ser. B Stat. Methodol. 71, 505–523 (2009).
MathSciNet  MATH  Article  Google Scholar 

89.
Wood, S. N. Fast stable restricted maximum likelihood and marginal likelihood estimation of semiparametric generalized linear models. J. R. Stat. Soc. Ser. B Stat. Methodol. 73, 3–36 (2011).
MathSciNet  MATH  Article  Google Scholar 

90.
Beaudoin, A. et al. Mapping attributes of Canada’s forests at moderate resolution through k NN and MODIS imagery. Can. J. Res. 44, 521–532 (2014).
Article  Google Scholar  More