Ecology

The research in this study was performed in agreement with permission M45-14 issued by Malmö/Lund Ethical Committee for Animal Research, Sweden, which granted capture and blood sampling of wild birdsSamplesNine willow warblers, determined to be males (based on a wing length > 69 mm), were caught opportunistically with mist nets during the time of autumn migration in September 2016 at Krankesjön, 15 km East of Lund, Southern Sweden. While most of the individuals were phenotypically similar to willow warblers breeding in Southern Scandinavia, some were slightly larger and had a greyer plumage, which is more commonly seen in Northern Scandinavia12. The set of samples thus potentially contained willow warblers of each of the two major migratory phenotypes. Blood from each bird was collected through a puncture of the brachial vein and was stored in two vails containing SET buffer and 70% ethanol, respectively. An aliquot of the blood was used for DNA extraction with a phenol-chloroform protocol. From the extracted DNA, we genotyped the samples for two loci located on chromosomes 1 and 5, respectively (NBEA and FADS2)45,46, and for a bi-allelic marker within the divergent region on chromosome 3 (AFLP-ww1)47. Based on the genotyping results we selected two samples that were homozygous northern or homozygous southern for all three loci, respectively. We also included a sample from a chiffchaff Phylloscopus collybita (female) for de novo genome sequencing of a closely related outgroup species, as well as an additional willow warbler (DD81063, male) to confirm breakpoint differences with linked read sequencing. Both of these birds were opportunistically caught at the same site as above during autumn migration in 2019, and collection of blood followed the same approach as for the other birds.Optical mapsDNA from the northern and southern willow warbler was extracted from blood stored in ethanol using a Plug Lysis protocol (v.30026D; Bionano Genomics, CA, USA). The blood was first separated from the ethanol through gentle centrifugation and embedded in molten 2% agarose plugs (DNA plug kit; Bio-Rad, CA, USA). The solidified plugs were submerged in Lysis Buffer solution (Bionano Genomics) and 66.8 µl per ml Buffer Puregene Proteinase K (Qiagen,MD, USA) for 2 h at 50 °C. The plugs were subsequently washed in 1× Wash buffer (Bio-Rad DNA plug kit) followed by TE buffer. In the following step, the plugs were treated with RNase (Qiagen, 20 µl in 1 ml TE buffer) for 1 h at 37 °C, followed by another washing step using the same buffers as in the previous step. Next, the plugs were melted for 2 min at 70 °C and treated with GELase (Epicenter, WI, USA) for 45 min at 43 °C. The DNA was then purified from digested agarose using drop dialysis against TE buffer on a 0.1 µm dialysis membrane (MF-Millipore, Merck KGaA, Germany) for 2.5 h.Optical maps for each of the two samples were produced using Bionano Genomic’s commercial Irys system48. BspQ1 was determined to be the most suitable nicking enzyme after using the software LabelDensityCalculator v.1.3.0 and Knickers v.1.5.5 to analyze a previous short-read assembly13. Bionano Genomic’s IrysPrep Labeling-NLRS protocol (v.30024) was used for the NLRS reaction. For this step, DNA was treated with Nt.BspQ1 (NEB, MA, USA) to create single-stranded nicks in a molecule-specific pattern. These were then labeled with Bionano Genomic’s (CA, USA) labeling mix (NLRS kit), aided by Taq Polymerase (NEB), and repaired using Bionano Genomics’s repair mix (NLRS kit), in the presence of Thermopol Rxn buffer, NAD+, and Taq DNA Ligase (NEB). Finally, the DNA backbone was stained using DNA stain from Bionano Genomics’s NLRS kit. Each sample was then loaded on two IrysChips (Bionano Genomics) each, and the DNA with stained BspQ1 nicks was visualized using an Irys instrument, following Bionano Genomics’s Irys user guide (v.30047). This resulted in 200 and 182 Gb of data for the northern and southern sample, respectively.Genome maps were assembled de novo using Bionano Genomic’s in house software IrysView v.2.5.1, with noise parameter set to “autonoise” and using a human arguments xml file. The genome map was then further refined by re-assembling all data, but using the first assembly version as a reference. The final assemblies were both 1.3 Gb in total size, with an average coverage of 92.3 and 96.4×, and N50 of 0.93 and 0.95 Mb, for the northern and southern sample, respectively.Linked read sequencingFor the southern sample and sample DD81063, DNA for chromium sequencing (10× Genomics, CA, USA) was extracted from blood stored in SET buffer using a MagAttract HMW DNAkit (Qiagen) at Scilifelab, Stockholm, Sweden. For the northern sample the extraction for bionano optical maps was used. The libraries of the northern and southern sample were each sequenced on a separate lane of a HiSeqX (Illumina, CA, USA) and the DD81063 sample was sequenced on a NovaSeq6000 (Illumina). For all samples sequencing was performed using a 2 × 150 bp setup.Northern willow warbler de novo assemblyLibrary preparation for long read sequencing was done on DNA previously extracted for the optical map and followed Pacific Bioscience’s (CA, USA) standard protocol for 10–20 kb libraries. No shearing was performed prior to the library construction, but the library was size selected using the BluePippin pulse field size selection system (Sage Science, MA, USA), with a size cut-off >25 kb. The library was sequenced on eight SMRT cells on a Sequel platform (Pacific Biosciences). The sequencing yielded 63.66 Gbp of data comprised of 4,690,365 subreads with a mean length of 13,573 bp (range: 50–170,531 bp).The Pacbio reads were assembled de novo in HGAP449 in the SMRT Link package with default settings except for specifying an expected genome size of 1.2 Gbp and setting the polishing algorithm to “Arrow”. We ran Falcon unzip50 on the assembly to obtain partially phased primary contigs and fully phased haplotigs. Within the software, Arrow was used to polish the assembly using reads assigned to each haplotype. We evaluated two unzipped assemblies based on 30× or 40× coverage of seed reads in the preassembly step in HGAP4. A lower coverage threshold will lead to longer reads in the initial assembly step, which may increase the contiguity of the assembly, but will on the other hand, limit the number of reads that can be used in the phasing and polishing step. Although the unzipped assemblies were very similar, the 40× version was chosen for downstream analyses as it was slightly more contiguous and contained a higher number of single-copy bird orthologues as determined by BUSCO version 3.0.251.The assembly was further polished with Pilon 1.2252 with Illumina chromium reads from the same sample. The Illumina reads were mapped to the assembly using bwa version 0.7.17-r118853 and duplicated reads were marked using picardtools 2.10.3 (http://broadinstitute.github.io/picard). Pilon was run by only correcting indels and in total the software made 1,043,827 insertions and 275,457 deletions, respectively, of which the vast majority (94%) were single basepair changes. The Illumina polishing had a pronounced effect on the number of single-copy bird orthologues that could be detected in the primary contigs (Supplementary Table 1).For further assembly steps, we extracted the Illumina-polished primary Pacbio contigs (N = 2737, N50 of 2.1 Mb and a length of 1.29 Gb). These contigs showed an unexpectedly high level of duplicated single-copy orthologues (7.4%), which suggested partial or complete overlap between some contigs. As a first step to reduce the redundancy and increase the contiguity of the assembly, we hybridized the primary contigs to the optical map of the same sample using bionano solve version 3.2.2 (BioNano Genomics) with default settings except for specifying aggressive scaffolding parameters. The hybrid scaffolding resulted in 19 cuts to the bionano maps and 259 cuts to the Pacbio contigs and created 363 super-scaffolds. Most of the gaps between the contigs in the super-scaffolds were estimated to be negative (i.e., some overlap between sequences). However, in the hybrid assembly, sequences on either side of these gaps were not collapsed and thus formed false segmental duplications. To remedy this problem we extracted 304 sets of overlapping contigs (“supercontigs”) and used GAP5 in the staden package 2.0.0.b1154 to find potential joins between the contig ends. Using this approach, we merged contigs at 558 (87%) of the putative overlaps. The mean alignment length in the overlaps was 111 kb (range: 0.259–661 kb) with a mean sequence divergence of 3.28% (range: 0.31–15.55%). The highest divergence was caused by the presence of large indels. By trimming off one or both ends of the contigs at the gaps (mean 23 kb, range: 0.6–60 kb), we were able to close 23 further gaps. For the remainder of gaps, GAP5 failed to find potential joins between contigs or the ends supposed to be joined were considered to have too high divergence. The new assembly, including supercontigs consisted of 2401 contigs with an N50 of 6.5 Mb and had a considerably lower amount of duplicated single-copy genes (4.6% vs 7.4%).To further reduce the redundancy, we used the purge haplotig pipeline55 (downloaded 2019-02-15) to remove contigs that could be mapped over most of their length to larger contigs and that showed limited diploid coverage. We first estimated coverage by mapping the Pacbio subreads used for the de novo assembly with minimap2 version2.13-r86056 using default settings for Pacbio reads (-x map-pb). To minimize the loss of repetitive sequences that could be separated and scaffolded by the bionano optical map, we used the first bionano hybrid assembly (363 superscaffolds and 1500 cut and unscaffolded contigs) as a reference for mapping. From the mapped data we detected a clear haploid and diploid peak and set a threshold of diploid coverage above 34× and below 85×. Any scaffold where less than 80% of its positions had diploid coverage was considered a putative haplotig and was mapped to other scaffolds using minimap2 within the software. We removed 1209 scaffolds (mean size: 107,655 bp, range: 598–495,788 bp) with a coverage to the best hit of at least 70% (mean: 97.4%). Using this approach, we specifically excluded contigs that could not be incorporated in superscaffolds. However, we also removed three contigs that each entirely made up short superscaffolds that could be uniquely assigned to larger superscaffolds and that had a high degree of haploid coverage. At this stage, we also removed five additional contigs shorter than 1000 bp that were the result of cutting the assembly with the bionano optical map. This led to an assembly with 1187 contigs, a length of 1.1 Gbp and a N50 of 7.9 Mb. The filtered assembly showed a large reduction in single-copy orthologue bird genes (1.3 vs 4.6%).To provide an intermediate level of scaffolding to the optical map, we mapped the 10× chromium reads of the same sample to the assembly using bwa and used arcs version 1.0.557 and LINKS version 1.8.658 for scaffolding. Arcs was run with default settings except for enabling gap size estimation (–dist_est) and LINKS was run by setting the number of supporting links to at least 5 (-l = 5) and the maximum link ratio between the two best contig pairs to 0.3 (-a = 0.3). The scaffolding resulted in 739 scaffolds with a N50 of 16.4 Mb and a length 1.12 Gb.As a final scaffolding step, we hybridized the 10× chromium-Pacbio scaffolds to the bionano optical map using the same settings as before. The hybrid scaffolding made 23 cuts to the optical map, 122 cuts to the scaffolds and resulted in 497 scaffolds with an N50 of 16.8 Mb. Two contigs representing the divergent region on chromosome 1 had been scaffolded together by arcs but were separated and not re-scaffolded with other sequences in the bionano hybrid assembly. Since the mismatched end of the optical map was short, located at a large gap, and the gene order is the same as seen in other bird genomes, we decided to keep the scaffold generated by arcs.For this round of hybrid scaffolding, there were 52 gaps that were estimated to be negative. Using the same approach as when creating supercontigs, we were able to close 10 of these gaps. We additionally closed gaps using PBJelly59 from PBSuite 15.8.24 with default settings except for specifying –spanOnly –capturedOnly”. The software filled 97 gaps, extended one end of 12 gaps, extended both ends of 18 gaps and overfilled 28 gaps (extended both ends but detected no overlap despite the extension is larger than the predicted gap).We further checked for potential misjoins between scaffolds that originate from different chromosomes. To this end, we used SatsumaSynteny 2.060 to produce whole-genome alignments between the assembly and the genomes of chicken (version GRCg6a) and zebra finch (version taeGut3.2.4), both downloaded from Ensembl (www.ensembl.org). Using this approach, we detected a scaffold that showed good alignments to both chromosomes 10 and 23 in both of the other species. We considered this join unlikely and decided to split the scaffold.Next, we performed a second round of polishing with the 10× chromium Illumina data from the same sample. For this round, since we had fewer than 500 scaffolds, we used the longranger 2.1.14 align pipeline61 to map reads in a barcode-aware way. Pilon was then run with the same settings as before and resulted in the correction of 417,032 indels, of which 78.7% were single-basepair changes. The second round of polishing considerably increased the number of single-copy bird orthologues that could be identified in the assembly (Supplementary Table 1).The mitochondrial genome was not found in the original Pacbio genome assembly. We obtained this genome by adding the complete mitochondrial sequence from a previous short-read assembly13. We then used bwa to map the 10× chromium reads from the northern sample to the assembly and extracted alignments on the mitochondrial sequence. Next, freebayes was used with a haploid setting to detect differences present in the aligned reads. The raw variant file was filtered with vcftools for sites with a quality less than 30 and for two intervals with excessive read coverage (possibly reads from unassembled NUMTs). The filtered variant file contained 11 substitutions and three indels, and was used with bcftools version 1.1462 to create a new mitochondrial reference.For the extraction and removal of sequences in the different assembly steps we used kentUtils 370 (https://github.com/ucscGenomeBrowser/kent). Summary statistics for each assembly (e.g., N50) were calculated using the assemblathon_stats.pl script63.Southern willow warbler and chiffchaff de novo assembliesThe southern willow warbler and the chiffchaff were each sequenced on two lanes on a Sequel II (Pacific Biosciences) using a high-fidelity (HiFi) setup. Sequencing libraries for the southern willow warbler was prepared from a previous extraction used for optical maps (see above), whereas for the chiffchaff, DNA was extracted from blood using a Nanobind extraction kit (Circulomics, MD, USA). The southern willow sample yielded 2,576,876 HiFi reads with a mean length 19,303 bp and representing 49.7 Gbp. The chiffchaff sample yielded 2,612,165 HiFi reads with a mean length of 19,829 bp and representing 51.8 Gbp.The HiFi reads were assembled de novo using hifiasm version 0.15.5-r35064 with default settings and primary contigs were selected for downstream analyses. For the chiffchaff hifiasm assembly, we removed the first 6 Mb part of a contig overlapping with another contig and removed a short interval at the end of a contig containing adaptor sequences. For the southern willow warbler, the primary contigs (N = 540, Supplementary Table 1) were hybridized to the optical map of the same sample using the same pipeline as for the northern sample. Although we had access to chromium data from the same sample, we did not include it to perform an intermediate scaffolding step (as we did for the northern willow warbler assembly) because the long-read assembly was already highly contiguous. The hybridization step made 39 cuts to the contigs and 20 cuts to the optical maps, resulting in an assembly with 111 superscaffolds and 439 non-scaffolded contigs. We decided to ignore an optical map-supported fusion of contigs that mapped to separate chromosomes in other bird species, as this fusion was made in a large repetitive region. We further excluded a 45 bp sequence resulting from the hybrid assembly cutting and masked four short intervals containing adaptor sequences. The assembly of the mitochondrion in the southern assembly followed the same pipeline as used for the northern assembly (see above). In this case, 10 substitutions and two indels were added to the mitochondrial sequence from the previous short-read assembly based on alignments of linked reads from the southern sample.Repeat annotationWe used Repeatmodeler version 1.0.865 for de novo identification of repeats in the southern assembly. The repeats detected by repeatmodeler were combined with 1,023 bird-specific repeats into a custom library. Next, we used repeatmasker version 4.0.766 with the custom library and by using a more sensitive search (-s flag) to annotate repeats in the genome. Bedtools v2.29.267, together with the annotated repeats, was used to create a softmasked version of the southern assembly, which was used in the gene annotation step. The same repeat library was also used to annotate repeats in the de novo assembly of the northern sample. For the chiffchaff assembly we used the same annotation approach as for the southern willow warbler, but included a species-specific library generated with repeatmodeler, and also included a tandem-repeat associated sequence associated with the divergent regions on chromosomes 1 and 3 from the willow warbler library. Intervals with tandem repeats in divergent regions were also analyzed with tandem repeats finder version 4.0.968 using default settings except for specifying a maximum period size of 2000 bp.Duplicated intervals within divergent scaffolds were identified with Minimap2 and subsequently aligned with EMBOSS Stretcher 6.6.0 (https://www.ebi.ac.uk/Tools/psa/emboss_stretcher/).RNA sequencingWe used total RNA extracted from whole brain from six samples used in an earlier study quantifying differential expression in migratory and breeding willow warblers69 (Supplementary Table 3). The quality of the RNA was checked with a Bioanalyzer version 2100 (Agilent, CA, USA). All of the extractions had a RNA Integrity Number (RIN) of at least > 7.10. RNA libraries for sequencing were prepared using a TruSeq Stranded mRNA Sample prep kit with 96 dual indexes (Illumina) according to the instructions of the manufacturer with the exception of automating the protocols using an NGS workstation (Agilent) and using purification steps as described in Lundin et al70. and Borgström et al71. The raw RNA data was trimmed using cutadapt version 1.872 within Trim Galore version 0.4.0 (https://github.com/FelixKrueger/TrimGalore) with default settings.We used Stringtie version 1.3.373 to create transcripts from the RNAseq data. These transcripts were not used directly in the generation of gene models, but used in the manual curation step as potential alternative transcripts. For the software, we first mapped the reads with Hisat2 version 2.1.074 using default settings for stranded sequence libraries and downstream transcript analyses.Gene annotationWe used Augustus version 3.2.375 to create gene models using hints provided from RNAseq data and protein data from other bird species. For the RNAseq data, we mapped the trimmed reads to the assembly using STAR version 2.7.9a76. Accessory scripts in the Augustus package were used to filter the alignments for paired and uniquely mapped reads and for extracting intron hints. We additionally generated coverage wig files for each strand from the filtered alignment file using the software stranded-coverage (https://github.com/pmenzel/stranded-coverage) and used these as input for the august wig2hints.pl to generate exonpart hints.For homology evidence, we downloaded a set of bird proteins from NCBI (https://www.ncbi.nlm.nih.gov/). This data set included 49,673 proteins from chicken, 41,214 proteins from zebra finch and 38,619 proteins from great tit. We also downloaded an additional dataset from Uniprot (www.uniprot.org) that consisted of 3175 manually reviewed bird proteins and 204 and 12,263 bird proteins that were not manually reviewed but supported by protein or transcript data, respectively. The protein data was mapped to the genome using exonerate version 2.4.077. We used the script align2hints.pl from braker 2.1.678 to generate CDSpart, intron, start and stop hints from the data.Augustus was run with species-specific parameters (see training Augustus below) and with default settings except for specifying “softmasking=true”, “–alternatives-from-evidence=true”, “–UTR = on”, “–gff3=on” and “–allow_hinted_splicesites=atac”. In the extrinsic configuration file, we changed the malus for introns from 0.34 to 0.001, which increases the penalty for predicted introns that are not supported by the extrinsic data (RNAseq and protein hints). The prediction resulted in 28,491 genes and 35,389 transcripts.The Augustus-derived gene models were assigned names based on overlap with synteny-transferred zebra finch genes. For this purpose, we used SatsumaSynteny with default settings to obtain whole-genome alignments between our assembly and the zebra finch genome version bTaeGut1.4.pri79. Based on the alignment, we used kraken80 (downloaded 2020-04-14) to transfer the zebra finch genome annotations (NCBI Release 106) to the willow warbler assembly. We then extracted the CDS from the Augustus gene models and the kraken genes and used bedtools intersect to quantify the overlap. The gene models were also searched against the longest translation of each of the chicken, zebra finch and great tit Parus major genes used as evidence for the gene prediction step and against 86,131 swissprot vertebrate proteins using blastp 2.5.0+81 with an E value threshold of 1e−5. Gene models that were not annotated through synteny were assigned a gene name based on the blast results. Protein domains in the gene models were annotated with interproscan v 5.30–69.082. To reduce the number of false positive predictions we removed 5697 genes that were not supported by synteny to zebra finch genes, showed no significant similarity to vertebrate proteins or did not contain any annotated protein domains.We used Webapollo 2.6.583 to manually curate gene models in the previously identified divergent chromosome regions and in other regions where differences were present. In the curation step, we specifically validated the support for the coding sequence and the UTR and also removed genes that were likely to be pseudogenes based on a truncated coding sequence compared to homologous genes in other vertebrates, had no support from synteny in other bird species and/or that were located in repeat-rich regions.Training AugustusWe used a previous repeat-masked short-read assembly13 and the trimmed RNAseq data used in this study to obtain species-specific parameters for Augustus. The RNAseq data was assembled into transcripts using Trinity version 2.0.284 to create a de novo and a genome-guided assembly that together were comprised of 1,929,396 transcripts. The genome-guided transcript assembly was based on RNAseq mapped to the genome using GSNAP version 2016-07-1185 with default settings. We used PASA version 2.0.286 to create high-quality transcripts, which were imported into Webapollo. To assess the completeness of the transcripts, we compared them to synteny-transferred models from the chicken genome using Kraken. We selected 1249 transcripts that appeared complete, were not overlapping with other genes and showed less than 80% amino acid similarity to another gene in the training set. From this set, we excluded 21 genes that were giving initial training errors, which gave us a training set of 1228 genes. This gene set was randomly split into 1028 training genes and 200 genes used for testing. For training, we used the optimize_augustus.pl script with default settings except for the flag –UTR = on.Whole-genome resequencing and variant callingWe used the whole-genome resequencing data from nine samples of each migratory phenotype provided in Lundberg et al13. and sequenced an additional two high-coverage samples from each migratory phenotype (Supplementary Table 4). Sequencing libraries for the new samples were prepared with a TruSeq DNA PCR-Free kit (Illumina) with a targeted insert size of 670 bp or with a Truseq DNA nano (Illumina) with a targeted insert size of 350 bp. All of the new samples were sequenced on a HiSeqX (Illumina). The raw reads were trimmed with trimmomatic 0.3687 with the parameters “ILLUMINACLIP:TruSeq3-PE-2.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:30”.Quality-trimmed reads were mapped to the southern assembly using bwa mem with default settings except for specifying -M flag to ensure compatibility with the downstream duplicate removal steps and converted into binary alignment map (bam) files using samtools. For samples sequenced across multiple lanes, reads from each lane were mapped independently and the resulting bam files were merged with samtools. Read duplicates were removed with the markduplicates tool provided in picardtools.From the aligned whole-genome resequencing data set, we called variants with freebayes v1.1.0 using default settings and parallelizing the analyses of separate scaffolds using GNU parallel88. Vcflib version 2017-04-0489 was used to filter the raw set of variants for sites with quality score >30 and for alternate alleles that were supported by at least one read on each strand (SAF  > 0 & SAR  > 0) and had at least one read balanced to the right and the left (RPL  > 0 & RPR  > 0). Next, we used vcftools 0.1.1690 to filter genotypes with a coverage of at least 5x and removed sites a maximum of four genotypes missing in each of the populations. The variants were also filtered for collapsed repeats by removing sites with a mean coverage of more than twice the median mean coverage (30×). We next used vcflib to decompose haplotype calls and complex alleles into indels and SNPs and removed any variants that were overlapping with annotated repeats. This gave us a final of 51 million variants of which 45 million were bi-allelic SNPs. We used vcftools to calculate FST91 for each variant and for bi-allelic SNPs in non-overlapping windows of 10 kb. As many rare variants segregate in the willow warbler populations, which may downwardly bias differentiation estimates92, we focused on variants with a minor allele frequency of at least 0.1.Coverage for each resequenced sample was calculated in non-overlapping 1 kb windows using bedtools and only included properly paired reads with a mapping quality of at least 1. The raw coverage values for each sample were normalized by its median coverage across all windows.Structural variant callingWe used a combination of delly 0.9.193 and GraphTyper 2.7.494 to call structural variants in the resequenced samples. To identify a set of high confidence variants, we first mapped the long reads from the northern willow warbler to the southern assembly using minimap 2.22-r110156 with default settings for Pacbio reads and from the alignments called variants using delly. Next, GraphTyper was used to genotype the resequenced samples for the delly variants in the scaffolds containing the divergent chromosome regions. The raw set of variants were filtered to contain only sites with a “PASS” flag and, for each variant, the aggregated genotype, which is the genotype model out of breakpoint alignments and coverage that has the highest genotyping quality, was chosen for downstream analyses. Genetic differentiation (FST) was calculated in vcftools and variants with FST ≥ 0.7 between homozygotes in each divergent chromosome region were extracted and checked for overlap with genes and gene features using bedtools. To get more reliable differentiation estimates, we only included sites where at least 80% of the southern and northern homozygotes had genotypes.Inversion genotypes for resequenced samplesThe resequenced samples were assigned a genotype of southern and northern haplotypes for each of the divergent regions based on a multidimensional scaling (MDS)-based clustering in invclust95 of SNP array genotypes in Lundberg et al.13. To obtain genotypes of the SNPs included on the array in the resequenced samples, we mapped the SNP array probe sequences to the northern assembly using gmap and from the alignments extracted the positions of the focal SNPs. Next, we used freebayes to genotype the resequenced samples for these positions and plink version 1.996 to combine the genotypes with the genotypes from the SNP array. In the genotyping step, we also included mapped 10× chromium libraries for the northern and southern reference samples and the additional willow warbler sample. From the combined dataset, we extracted genotypes for SNPs located in each of the divergent regions and used invclust to assign each sample a genotype of inverted and non-inverted haplotypes. The inverted and non-inverted haplotypes were recoded as southern or northern haplotypes based on their frequency in each subspecies.Breakpoint analysesWe used MUMmer 4.0.0rc197 to align the genomes of the southern and northern willow warblers, and the southern willow warbler genome to the genomes of the chiffchaff, zebra finch (3.2.4) and collared flycatcher FicAlb (1.5)98.To provide further evidence of breakpoints, we mapped the 10× chromium reads of each sample to both the northern and the southern assembly and called structural variants using the longranger wgs pipeline. For the southern genome, we selected the 499 largest scaffolds and concatenated the rest into a single scaffold to make it compatible with the software. We also checked for differences in linked read molecule coverage between the samples. For this purpose, the raw reads of each sample were first processed with longranger basic for quality trimming and barcode processing. The trimmed reads were mapped to the assemblies using bwa mem using a -C flag to extract the barcode information of each read and alignments converted into bam files using samtools. To estimate coverage of barcodes, we first used the tigmint-molecule script from tigmint 1.1.299 to obtain positional information of barcodes (molecules) in each divergent region. The software was run with default settings except for only using reads with a mapping quality of at least 1 and only to report molecules that were estimated to be at least 10 kb. We next used bedtools to count the number of overlapping molecules in 1 kb windows.We explored differences between optical maps by using the runSV.py script in bionano solve with the southern optical map as a query and the northern assembly as target and the reciprocal analysis with the northern optical map as a query and the southern assembly as a target. We also used the bionano solve hybrid assembly pipeline to visualize differences between the optical maps and the genome assemblies at breakpoint regions.Functional annotation of differencesWe used bedtools to quantify the distance between breakpoint intervals and annotated genes. To provide a functional annotation of the SNPs and short indels, we selected variants that showed a FST ≥ 0.7 between southern and northern homozygotes for each of the region and used these as input to Snpeff 5.0.0e100 together with the annotation and reference genome. We used Snpsift 5.0.0e101 to select variants that were predicted to have a moderate to high effect on genes. Gene ontology terms for the genes were extracted from orthologous genes in other bird genomes in ensembl (www.ensembl.org) or through domain searches of the proteins with interproscan.Age estimation and demographic analyses of divergent regionsIn order to estimate the timing of the inversion events, we used high-coverage resequencing data from two southern samples, two northern samples and, as an outgroup, one dusky warbler Phylloscopus fuscatus (Supplementary Table 4). The willow warbler samples were chosen so that they were either homozygous southern or northern for all of three divergent regions. The dusky warbler library was prepared using a TruSeq Nano DNA library prep kit for Neoprep (Illumina) according to the instructions of the manufacturer and sequenced on a HiSeq X (Illumina). Quality-trimming of the raw reads and mapping of the trimmed reads to the northern reference genome followed the same approach as used for the willow warbler resequencing samples (see above).Variants were called using freebayes and the raw set of variants were filtered using gIMble’s preprocess module (v0.6.0). Sample-specific callable sites were identified using gIMble preprocess and were defined as those with a minimum coverage of 8× and a maximum of 0.75 standard deviations above the mean coverage. Genic and repetitive regions of the genome were removed from the callable sites in order to limit downstream analyses to intergenic regions.Summary statistics of genetic variation (π and dxy) within the divergent regions were calculated using gIMble. Following this, net divergence (da) between northern and southern samples was calculated as dnorth–south − (πnorth + πsouth)/2. To convert the net divergence into years we used the germline mutation rate (4.6 × 10−9) estimated in the collared flycatcher21. Relative node depth (RND) using the dusky warbler (DW) as an outgroup was calculated as dnorth–south/(dDW-north + dDW-south)/2. For each divergent region, a blockwise site frequency spectrum (bSFS) was generated with gIMble using blocks of 64 bp in length. This length refers to the number of callable sites within a block, while the physical length of blocks was allowed to vary due to missing data but was limited to 128 bp. Downstream analyses that relied on a bSFS used a kmax of 2, meaning that only marginal probabilities were calculated for mutation counts >2. The composite likelihood (CL) of a model, given the bSFS of one of the divergent regions, was optimized using the Nelder-Mead algorithm with the maximum number of iterations set to 1000. Within the software we evaluated three different population models. The first model was a strict isolation model (SI), with parameters ancestral effective population size, effective population sizes for southern and northern willow warblers and divergence time. The second model was an isolation with migration model (IM1) that also included a migration rate from northern to southern samples, and the third model (IM2) instead had a migration rate from southern to northern willow warblers.Simulations were carried out by msprime 0.7.4102 through gIMble. The recombination rates used for these simulations were chromosome-specific estimates from a high-density recombination map of the collared flycatcher98 and were 2.04, 1.95, and 2.63 cM/Mb for chromosomes 1, 3, and 5, respectively. A total of 100 replicates were simulated for the optimized SI parameters of each region. These simulated bSFSs were then optimized under both an SI model as well as the best fitting IM model for that region. The improvement in CL between these models was used as a null distribution for testing whether improvements in CL observed for the real data were greater than expected given a history of no migration. For each parameter, we calculated 95% CI as Maximum Composite Likelihood (MCL) estimate ± 1.96 * standard deviation of simulations (Supplementary Table 7). As a result, our estimates of uncertainty are affected by the recombination rates that we assumed for simulations. We also used the results of simulations to quantify the potential bias in MCL estimates due to intra-block recombination (Supplementary Table 7). However, we did not attempt to correct for this bias as it is relatively small (e.g., the MCL divergence times are estimated to be biased upwards by 7, 24, and 10%) and our estimation of the bias itself is largely dependent on the recombination rates we assumed.MSMC224 was used to explore genome-wide changes in Ne through time. As input to the software, we used the callable intergenic bed file and filtered vcf file mentioned above, with the addition of further filtering the bed file to only include autosomal scaffolds ≥500 kb and excluding the divergent regions. The input files for MSMC2, i.e., an unphased set of heterozygous sites for each sample, were generated using the generate_multihetsep.py script from msmc-tools. MSMC2 was run with a starting ρ/μ of 1 for 30 expectation-maximum iterations. For both the demographic modeling and MSMC2, we used the collared flycatcher germline mutation rate21 and a generation time of 1.7 years11 to convert divergence times into years.To infer the effects of demographic events and selection, we also calculated several genetic summary statistics. To this end, we first imputed missing genotypes and inferred haplotypes for the filtered set of variants using beagle version 5.4103. From the full set of samples, we selected 10 and seven samples that were homozygous southern or northern for the three divergent regions, respectively, as determined from the MDS analysis (see above), and extracted bi-allelic SNPs. To identify ancestral and derived alleles, we extracted genotypes for the focal SNP positions from the aligned chiffchaff and dusky warblers reads using bcftools 1.1462 with the mpileup command. As a conservative approach, we considered any site with the presence of both the reference and alternate allele as heterozygous (regardless of their frequencies) and only included sites where the coverage was at least one-third of the mean coverage among all sites for each outgroup species. We next used a customized script to extract the sites from the original vcf files, and, if necessary, switch the reference and alternate allele and swap the genotypes accordingly. With the polarized genotype data, we used PopGenome 2.7.5104 to calculate Fay and Wu’s H and vcftools to get counts for the derived allele. We further used selscan 1.3.0105 to calculate XP-nsl106 between the southern and northern samples, Sweepfinder2107 to calculate a composite likelihood ratio (CLR) between a model where a selective sweep has had an effect on the allele frequency and a model based on the genome-wide allele frequency spectrum and used vcftools to calculate nucleotide diversity, Tajima’s D and linkage disequilibrium (D’).The use of the southern assembly as a reference could potentially lead to a mapping bias for reads from southern samples, particularly in regions of higher divergence between the subspecies. This, in turn, could have an effect on genetic summary statistics and demographic modeling estimates. To explore the effect of reference bias, we therefore also mapped the resequencing data to the northern assembly, performed variant calling and calculated nucleotide diversity and Tajima’s D in 10 kb windows. For the northern assembly, we also used the same demographic modeling as used for the southern assembly. Contrasting average genetic summary statistics and demographic parameter estimates, we found negligible differences between the two genome assemblies (Supplementary Table 10).Reporting summaryFurther information on research design is available in the Nature Portfolio Reporting Summary linked to this article. More

Mohn WW, Mertens B, Neufeld JD, Verstraete W, de Lorenzo V. Distribution and phylogeny of hexachlorocyclohexane-degrading bacteria in soils from Spain. Environ Microbiol. 2006;8:60–8.CAS 

Google Scholar 
Sharma P, Raina V, Kumari R, Malhotra S, Dogra C, Kumari H, et al. Haloalkane Dehalogenase LinB is responsible for β- and δ-hexachlorocyclohexane transformation in Sphingobium indicum B90A. Appl Environ Microbiol. 2006;72:5720–7.CAS 

Google Scholar 
Lal R, Dogra C, Malhotra S, Sharma P, Pal R. Diversity, distribution and divergence of lin genes in hexachlorocyclohexane-degrading sphingomonads. Trends Biotechnol. 2006;24:121–30.CAS 

Google Scholar 
Singh A, Lal R. Sphingobium ummariense sp. nov., a hexachlorocyclohexane (HCH)-degrading bacterium, isolated from HCH-contaminated soil. Int J Syst Evol Microbiol. 2009;59:162–6.CAS 

Google Scholar 
Singh BK, Kuhad RC. Biodegradation of lindane (γ-hexachlorocyclohexane) by the white-rot fungus Trametes hirsutus. Lett Appl Microbiol. 1999;28:238–41.CAS 

Google Scholar 
Kumar D, Pannu R. Perspectives of lindane (γ-hexachlorocyclohexane) biodegradation from the environment: a review. Bioresour Bioprocess. 2018;5:29.CAS 

Google Scholar 
Lal R, Pandey G, Sharma P, Kumari K, Malhotra S, Pandey R, et al. Biochemistry of microbial degradation of hexachlorocyclohexane and prospects for bioremediation. Microbiol Mol Biol Rev. 2010;74:58–80.CAS 

Google Scholar 
Nayyar N, Lal R. Hexachlorocyclohexane contamination and solutions: Brief history and beyond. J Bioremed Biodeg. 2016;07:338.Bumpus JA, Tien M, Wright D, Aust SD. Oxidation of persistent environmental pollutants by a white rot fungus. Science. 1985;228:1434–6.CAS 

Google Scholar 
Phillips TM, Seech AG, Lee H, Trevors JT. Biodegradation of hexachlorocyclohexane (HCH) by microorganisms. Biodegr. 2005;16:363–92.CAS 

Google Scholar 
Dogra C, Raina V, Pal R, Suar M, Lal S, Gartemann KH, et al. Organization of lin genes and IS6100 among different strains of hexachlorocyclohexane-degrading Sphingomonas paucimobilis: Evidence for horizontal gene transfer. J Bacteriol. 2004;186:2225–35.CAS 

Google Scholar 
Pal R, Bala S, Dadhwal M, Kumar M, Dhingra G, Prakash O, et al. Hexachlorocyclohexane-degrading bacterial strains Sphingomonas paucimobilis B90A, UT26 and Sp+, having similar lin genes, represent three distinct species, Sphingobium indicum sp. nov., Sphingobium japonicum sp. nov. and Sphingobium francense sp. nov., and reclassification of [Sphingomonas] chungbukensis as Sphingobium chungbukense comb. nov. Int J Syst Evol Microbiol. 2005;55:1965–72.CAS 

Google Scholar 
Rijnaarts HHM, Bachmann A, Jumelet JC, Zehnder AJB. Effect of desorption and intraparticle mass transfer on the aerobic biomineralization of. alpha-hexachlorocyclohexane in a contaminated calcareous soil. Environ Sci Technol. 1990;24:1349–54.CAS 

Google Scholar 
Harms H, Schlosser D, Wick LY. Untapped potential: exploiting fungi in bioremediation of hazardous chemicals. Nat Rev Microbiol. 2011;9:177–92.CAS 

Google Scholar 
Wick LY Bioavailability as a microbial system property: Lessons learned from biodegradation in the mycosphere. In: Ortega-Calvo JJ, Parsons JR, (eds). Bioavailability of organic chemicals in soil and sediment. Cham: Springer International Publishing; 2020 p. 267–89.Jennings DH. Translocation of Solutes in fungi. Biol Rev. 1987;62:215–43.CAS 

Google Scholar 
Nazir R, Warmink JA, Boersma H, van Elsas JD. Mechanisms that promote bacterial fitness in fungal-affected soil microhabitats. FEMS Microbiol Ecol. 2010;71:169–85.CAS 

Google Scholar 
Worrich A, Stryhanyuk H, Musat N, König S, Banitz T, Centler F, et al. Mycelium-mediated transfer of water and nutrients stimulates bacterial activity in dry and oligotrophic environments. Nat Commun. 2017;8:15472.CAS 

Google Scholar 
Deveau A, Bonito G, Uehling J, Paoletti M, Becker M, Bindschedler S, et al. Bacterial–fungal interactions: ecology, mechanisms and challenges. FEMS Microbiol Rev. 2018;42:335–52.CAS 

Google Scholar 
Kohlmeier S, Smits THM, Ford RM, Keel C, Harms H, Wick LY. Taking the fungal highway: Mobilization of pollutant-degrading bacteria by fungi. Environ Sci Technol. 2005;39:4640–6.CAS 

Google Scholar 
Wick LY, Furuno S, Harms H Fungi as transport vectors for contaminants and contaminant-degrading bacteria. In: Timmis KN, (eds). Handbook of hydrocarbon and lipid microbiology. Berlin, Heidelberg: Springer Berlin Heidelberg; 2010. p. 1555–61.Espinosa-Ortiz EJ, Rene ER, Gerlach R. Potential use of fungal-bacterial co-cultures for the removal of organic pollutants. Crit Rev Biotechnol. 2021;0:1–23.
Google Scholar 
Haq IU, Zhang M, Yang P, van Elsas JD The interactions of bacteria with fungi in soil. In: Adv Appl Microbiol. Elsevier; 2014. p. 185–215.Furuno S, Foss S, Wild E, Jones KC, Semple KT, Harms H, et al. Mycelia promote active transport and spatial dispersion of polycyclic aromatic hydrocarbons. Environ Sci Technol. 2012;46:5463–70.CAS 

Google Scholar 
Schamfuß S, Neu TR, van der Meer JR, Tecon R, Harms H, Wick LY. Impact of mycelia on the accessibility of fluorene to PAH-degrading bacteria. Environ Sci Technol. 2013;47:6908–15.
Google Scholar 
Schamfuß S, Neu TR, Harms H, Wick LY. A Whole Cell Bioreporter approach to assess transport and bioavailability of organic contaminants in water unsaturated systems. J Vis Exp. 2014;24:52334.
Google Scholar 
You X, Kallies R, Kühn I, Schmidt M, Harms H, Chatzinotas A, et al. Phage co-transport with hyphal-riding bacteria fuels bacterial invasion in a water-unsaturated microbial model system. ISME J. 2022;16:1275–83.CAS 

Google Scholar 
Wick LY, Remer R, Würz B, Reichenbach J, Braun S, Schäfer F, et al. Effect of fungal hyphae on the access of bacteria to phenanthrene in soil. Environ Sci Technol. 2007;41:500–5.CAS 

Google Scholar 
Boer W, de, Folman LB, Summerbell RC, Boddy L. Living in a fungal world: impact of fungi on soil bacterial niche development. FEMS Microbiol Rev. 2005;29:795–811.
Google Scholar 
Frey‐Klett P, Garbaye J, Tarkka M. The mycorrhiza helper bacteria revisited. N. Phytol. 2007;176:22–36.
Google Scholar 
Hazen TC Cometabolic Bioremediation. In: Timmis KN, editor. Handbook of hydrocarbon and lipid microbiology. Berlin, Heidelberg: Springer; 2010. p. 2505–14.Jehmlich N, Vogt C, Lünsmann V, Richnow HH, von Bergen M. Protein-SIP in environmental studies. Curr Opin Biotechnol. 2016;41:26–33.CAS 

Google Scholar 
Sachsenberg T, Herbst FA, Taubert M, Kermer R, Jehmlich N, von Bergen M, et al. MetaProSIP: automated inference of stable isotope incorporation rates in proteins for functional metaproteomics. J Proteome Res. 2015;14:619–27.CAS 

Google Scholar 
Taubert M, Jehmlich N, Vogt C, Richnow HH, Schmidt F, von Bergen M, et al. Time resolved protein-based stable isotope probing (Protein-SIP) analysis allows quantification of induced proteins in substrate shift experiments. Proteomics 2011;11:2265–74.CAS 

Google Scholar 
Lünsmann V, Kappelmeyer U, Benndorf R, Martinez-Lavanchy PM, Taubert A, Adrian L, et al. In situ protein-SIP highlights Burkholderiaceae as key players degrading toluene by para ring hydroxylation in a constructed wetland model: Protein-SIP in a toluene-degrading wetland. Environ Microbiol. 2016;18:1176–86.
Google Scholar 
Khan N, Toscan RB, Lunayo A, Wamalwa B, Muge E, Mulaa FJ, et al. Draft genome sequence of Fusarium equiseti K3, a fungal species isolated from hexachlorocyclohexane-contaminated soil. Microbiol Resour Announc. 2021;10:e00885–21.CAS 

Google Scholar 
Khan N, Brizola Toscan R, Lunayo A, Wamalwa B, Muge E, Mulaa FJ, et al. Draft genome sequences of two Sphingobium species associated with hexachlorocyclohexane (HCH) degradation isolated from an HCH-Contaminated Soil. Microbiol Resour Announc. 2022;11:e00886–21.
Google Scholar 
Yang CH, Menge JA, Cooksey DA. Mutations Affecting Hyphal Colonization and Pyoverdine Production in Pseudomonads Antagonistic toward Phytophthora parasitica. Appl Environ Microbiol. 1994;9:473–81.Senoo K, Wada H. Isolation and identification of an aerobic γ-HCH-decomposing bacterium from soil. Soil Sci Plant Nutr. 1989;35:79–87.CAS 

Google Scholar 
Jehmlich N, Schmidt F, Taubert M, Seifert J, Bastida F, von Bergen M, et al. Protein-based stable isotope probing. Nat Protoc. 2010;5:1957–66.CAS 

Google Scholar 
Graves S, Dorai-Raj HPP and LS with help from S. multcompView: Visualizations of paired comparisons. 2019. Available from: https://CRAN.R-project.org/package=multcompViewWickham H, François R, Henry L, Müller K, RStudio. dplyr: A Grammar of data manipulation. 2022. Available from: https://CRAN.R-project.org/package=dplyrggplot2: Create elegant data visualisations using the grammar of graphics—ggplot2-package. Available from: https://ggplot2.tidyverse.org/reference/ggplot2-package.htmlSchramm FD, Heinrich K, Thüring M, Bernhardt J, Jonas K. An essential regulatory function of the DnaK chaperone dictates the decision between proliferation and maintenance in Caulobacter crescentus. PLoS Genet. 2017;13:e1007148.
Google Scholar 
Li SX, Wu HT, Liu YT, Jiang YY, Zhang YS, Liu WD, et al. The F1Fo-ATP synthase β subunit is required for Candida albicans pathogenicity due to its role in carbon flexibility. Front Microbiol. 2018;9:1025.
Google Scholar 
Godlewska R, Wiśniewska K, Pietras Z, Jagusztyn-Krynicka EK. Peptidoglycan-associated lipoprotein (Pal) of Gram-negative bacteria: function, structure, role in pathogenesis and potential application in immunoprophylaxis. FEMS Microbiol Lett. 2009;298:1–11.CAS 

Google Scholar 
Taubert M, Vogt C, Wubet T, Kleinsteuber S, Tarkka MT, Harms H, et al. Protein-SIP enables time-resolved analysis of the carbon flux in a sulfate-reducing, benzene-degrading microbial consortium. ISME J. 2012;6:2291–301.CAS 

Google Scholar 
Little AEF, Robinson CJ, Peterson SB, Raffa KF, Handelsman J. Rules of engagement: Interspecies interactions that regulate microbial communities. Annu Rev Microbiol. 2008;62:375–401.CAS 

Google Scholar 
Morris BEL, Henneberger R, Huber H, Moissl-Eichinger C. Microbial syntrophy: interaction for the common good. FEMS Microbiol Rev. 2013;37:384–406.CAS 

Google Scholar 
Van Hees PAW, Rosling A, Essén S, Godbold DL, Jones DL, Finlay RD. Oxalate and ferricrocin exudation by the extramatrical mycelium of an ectomycorrhizal fungus in symbiosis with Pinus sylvestris. N Phytol. 2006;169:367–78.
Google Scholar 
Sun YP, Unestam T, Lucas SD, Johanson KJ, Kenne L, Finlay R. Exudation-reabsorption in a mycorrhizal fungus, the dynamic interface for interaction with soil and soil microorganisms. Mycorrhiza 1999;9:137–44.CAS 

Google Scholar 
Leveau JHJ, Preston GM. Bacterial mycophagy: definition and diagnosis of a unique bacterial-fungal interaction. N Phytol. 2008;177:859–76.
Google Scholar 
Bassler BL, Losick R. Bacterially speaking. Cell 2006;125:237–46.CAS 

Google Scholar 
Kim HJ, Boedicker JQ, Choi JW, Ismagilov RF. Defined spatial structure stabilizes a synthetic multispecies bacterial community. Proc Natl Acad Sci. 2008;105:18188–93.CAS 

Google Scholar 
Nayyar N, Sangwan N, Kohli P, Verma H, Kumar R, Negi V, et al. Hexachlorocyclohexane: persistence, toxicity and decontamination. Rev Environ Health. 2014;29:49–52.CAS 

Google Scholar 
Willett KL, Ulrich EM, Hites RA. Differential toxicity and environmental fates of hexachlorocyclohexane isomers. Environ Sci Technol. 1998;32:2197–207.CAS 

Google Scholar 
Ellegaard-Jensen L, Knudsen BE, Johansen A, Albers CN, Aamand J, Rosendahl S. Fungal-bacterial consortia increase diuron degradation in water-unsaturated systems. Sci Total Environ. 2014;466–467:699–705.
Google Scholar 
Knudsen BE, Ellegaard-Jensen L, Albers CN, Rosendahl S, Aamand J. Fungal hyphae stimulate bacterial degradation of 2,6-dichlorobenzamide (BAM). Environ Pollut. 2013;181:122–7.CAS 

Google Scholar 
Saez JM, Alvarez A, Fuentes MS, Amoroso MJ, Benimeli CS An Overview on Microbial Degradation of Lindane. In: Singh SN, (eds). Microbe-induced degradation of pesticides. Cham: Springer International Publishing; 2017. p. 191–212.Nzila A. Update on the cometabolism of organic pollutants by bacteria. Environ Pollut. 2013;178:474–82.CAS 

Google Scholar 
Benimeli CS, González AJ, Chaile AP, Amoroso MJ. Temperature and pH effect on lindane removal by Streptomyces sp. M7 in soil extract. J Basic Microbiol. 2007;47:468–73.CAS 

Google Scholar 
Álvarez A, Yañez ML, Benimeli CS, Amoroso MJ. Maize plants (Zea mays) root exudates enhance lindane removal by native Streptomyces strains. Int Biodeterior Biodegrad. 2012;66:14–8.
Google Scholar 
Boltner D, Moreno-Morillas S, Ramos JL. 16S rDNA phylogeny and distribution of lin genes in novel hexachlorocyclohexane-degrading Sphingomonas strains. Environ Microbiol. 2005;7:1329–38.CAS 

Google Scholar  More

Elrick, M. et al. Major Early-Middle Devonian oceanic oxygenation linked to early land plant evolution detected using high-resolution U isotopes of marine limestones. Earth Planet. Sci. Lett. 581, 117410 (2022).CAS 

Google Scholar 
Algeo, T. J. & Scheckler, S. E. Terrestrial-marine teleconnections in the Devonian: Links between the evolution of land plants, weathering processes, and marine anoxic events. Philos. Trans. R. Soc. Lond. Ser. B Biol. Sci. 353(1365), 113–130 (1998).
Google Scholar 
Capel, E. et al. The Silurian-Devonian terrestrial revolution: Diversity patterns and sampling bias of the vascular plant macrofossil record. Earth Sci. Rev. 231, 104085 (2022).
Google Scholar 
Racki, G., Joachimski, M. M. & Morrow, J. R. A major perturbation of the global carbon budget in the Early-Middle Frasnian transition (Late Devonian). Palaeogeogr. Palaeoclimatol. Palaeoecol. 269(3–4), 127–129 (2008).
Google Scholar 
Stein, W. E., Berry, C. M., Hernick, L. V. & Mannolini, F. Surprisingly complex community discovered in the mid-Devonian fossil forest at Gilboa. Nature 483(7387), 78–81 (2012).ADS 
CAS 

Google Scholar 
Retallack, G. J. & Huang, C. Ecology and evolution of Devonian trees in New York, USA. Palaeogeogr. Palaeoclimatol. Palaeoecol. 299(1–2), 110–128 (2011).
Google Scholar 
Qie, W., Algeo, T. J., Luo, G. & Herrmann, A. Global events of the late Paleozoic (early Devonian to Middle Permian): A review. Palaeogeogr. Palaeoclimatol. Palaeoecol. 531, 109259 (2019).
Google Scholar 
Smart, M. S., Filippelli, G., Gilhooly III, W. P., Marshall, J. E. & Whiteside, J. H. Enhanced terrestrial nutrient release during the Devonian emergence and expansion of forests: Evidence from lacustrine phosphorus and geochemical records. GSA Bulletin. Nov. 9 (2022).Śliwiński, M. G., Whalen, M. T. & Day, J. Trace element variations in the Middle Frasnian punctata zone (Late Devonian) in the western Canada sedimentary basin— changes in oceanic bioproductivity and paleoredox spurred by a pulse of terrestrial afforestation?. Geol. Belg. 4, 459–482 (2010).
Google Scholar 
Filippelli, G. M. & Souch, C. Effects of climate and landscape development on the terrestrial phosphorus cycle. Geology 27(2), 171–174 (1999).ADS 
CAS 

Google Scholar 
Filippelli, G. M., Souch, C., Horn, S. P. & Newkirk, D. The pre-Colombian footprint on terrestrial nutrient cycling in Costa Rica: Insights from phosphorus in a lake sediment record. J. Paleolimnol. 43(4), 843–856 (2010).ADS 

Google Scholar 
Pisarzowska, A. & Racki, G. Comparative carbon isotope chemostratigraphy of major Late Devonian biotic crises. In Stratigraphy & Timescales. 387–466, vol. 5. (Academic Press, 2020).Mortlock, R. A. & Froelich, P. N. A simple method for the rapid determination of biogenic opal in pelagic marine sediments. Deep Sea Res. Part A Oceanogr. Res. Pap. 36(9), 1415–1426 (1989).ADS 
CAS 

Google Scholar 
Schieber, J., Krinsley, D. & Riciputi, L. Diagenetic origin of quartz silt in mudstones and implications for silica cycling. Nature 406(6799), 981–985 (2000).ADS 
CAS 

Google Scholar 
Buckman, J., Mahoney, C., März, C., Wagner, T. & Blanco, V. Identifying biogenic silica: Mudrock micro-fabric explored through charge contrast imaging. Am. Miner. 102(4), 833–844 (2017).ADS 

Google Scholar 
Gao, P., He, Z., Lash, G. G., Zhou, Q. & Xiao, X. Controls on silica enrichment of Lower Cambrian organic-rich shale deposits. Mar. Pet. Geol. 130, 105126 (2021).CAS 

Google Scholar 
Schieber, J. Early diagenetic silica deposition in algal cysts and spores; a source of sand in black shales?. J. Sediment. Res. 66(1), 175–183 (1996).
Google Scholar 
Śliwiński, M. G., Whalen, M. T., Newberry, R. J., Payne, J. H. & Day, J. E. Stable isotope (δ13Ccarb and org, δ15Norg) and trace element anomalies during the Late Devonian ‘punctata Event’in the Western Canada Sedimentary Basin. Palaeogeogr. Palaeoclimatol. Palaeoecol. 307(1–4), 245–271 (2011).
Google Scholar 
Papazis, P. K. & Milliken, K. Cathodoluminescent textures and the origin of quartz in the Mississippian Barnett Shale, Fort Worth Basin, Texas. In AAPG Annual Meeting, Volume Abstracts: Calgary, Alberta, American Association of Petroleum Geologists A, 105 (2005).Ross, D. J. & Bustin, R. M. Investigating the use of sedimentary geochemical proxies for paleoenvironment interpretation of thermally mature organic-rich strata: Examples from the Devonian-Mississippian shales, Western Canadian Sedimentary Basin. Chem. Geol. 260(1–2), 1–19 (2009).ADS 
CAS 

Google Scholar 
Götze, J., Plötze, M. & Habermann, D. Origin, spectral characteristics and practical applications of the cathodoluminescence (CL) of quartz–a review. Mineral. Petrol. 71(3), 225–250 (2001).ADS 

Google Scholar 
Milliken, K. L., Ergene, S. M. & Ozkan, A. Quartz types, authigenic and detrital, in the Upper Cretaceous Eagle Ford Formation, south Texas, USA. Sed. Geol. 339, 273–288 (2016).CAS 

Google Scholar 
Blatt, H. Perspectives; Oxygen isotopes and the origin of quartz. J. Sediment. Res. 57(2), 373–377 (1987).ADS 
CAS 

Google Scholar 
Rowe, H. D., Loucks, R. G., Ruppel, S. C. & Rimmer, S. M. Mississippian Barnett Formation, Fort Worth Basin, Texas: Bulk geochemical inferences and Mo–TOC constraints on the severity of hydrographic restriction. Chem. Geol. 257(1–2), 16–25 (2008).ADS 
CAS 

Google Scholar 
Wright, A. M., Ratcliffe, K. T., Zaitlin, B. A. & Wray, D. S. The application of chemostratigraphic techniques to distinguish compound incised valleys in low-accommodation incised-valley systems in a foreland-basin setting: An example from the Lower Cretaceous Mannville Group and Basal Colorado Sandstone (Colorado Group), Western Canadian Sedimentary Basin, in K.T. Ratcliffe, and B.A. Zaitlin (eds.), Application of Modern Stratigraphic Techniques: Theory and Case Histories: SEPM SP PUB no. 94 (2010).Murata, K. J. & Norman, M. B. An index of crystallinity for quartz. Am. J. Sci. 276(9), 1120–1130 (1976).ADS 
CAS 

Google Scholar 
Tréguer, P. J. et al. Reviews and syntheses: The biogeochemical cycle of silicon in the modern ocean. Biogeosciences 18(4), 1269–1289 (2021).ADS 

Google Scholar 
Rivard, B., Harris, N. B., Feng, J. & Dong, T. Inferring total organic carbon and major element geochemical and mineralogical characteristics of shale core from hyperspectral imagery. AAPG Bull. 102(10), 2101–2121 (2018).
Google Scholar 
Lippincott, E. R., Van Valkenburg, A., Weir, C. E. & Bunting, E. N. Infrared studies on polymorphs of silicon dioxide and germanium dioxide. J. Res. Natl. Bur. Stand 61(1), 61–70 (1958).CAS 

Google Scholar 
Salisbury, J. W., D’Aria, D. M. & Jarosewich, E. Midinfrared (2.5–13.5 μm) reflectance spectra of powdered stony meteorites. Icarus 92(2), 280–297 (1991).ADS 
CAS 

Google Scholar 
Wong, P. K., Weissenberger, J. A. W., Gilhooly, M. G., Playton, T. E. & Kerans, C. Revised regional Frasnian sequence stratigraphic framework, Alberta outcrop and subsurface. New Adv. Devonian Carbonates: Outcrop Analogs, Reservoirs, and Chronostratigr. 49(1), 37–85 (2016).
Google Scholar 
Wendte, J. C. Cooking Lake platform evolution and its control on Late Devonian Leduc reef inception and localization, Redwater, Alberta. Bull. Can. Pet. Geol. 42(4), 499–528 (1994).
Google Scholar 
Wendte, J., Stoakes, F. A. & Campbell, C. V. Cyclicity of Devonian strata in the Western Canada Sedimentary Basin. In: Devonian-Early Mississippian Carbonates of the Western Canada Sedimentary Basin: A sequence stratigraphic framework. J. Wendte (ed.). Society of Economic Paleontologists and Mineralogists, Short Course no. 28, p. 25–40 (1995).Stoakes, F. A. Nature and control of shale basin fill and its effect on reef growth and termination: Upper Devonian Duvernay and Ireton Formations of Alberta, Canada. Bull. Can. Pet. Geol. 28(3), 345–410 (1980).
Google Scholar 
Alberta Energy Regulator Duvernay Reserves and Resources Report: A Comprehensive Analysis of Alberta’s Foremost Liquids-Rich Shale Resource, December 2016.Knapp, L. J., McMillan, J. M. & Harris, N. B. A depositional model for organic-rich Duvernay Formation mudstones. Sed. Geol. 347, 160–182 (2017).CAS 

Google Scholar 
Andrichuk, J. M. Stratigraphic evidence for tectonic and current control of Upper Devonian reef sedimentation, Duhamel area, Alberta, Canada. AAPG Bull. 45(5), 612–632 (1961).
Google Scholar 
Harris, N. B., McMillan, J. M., Knapp, L. J. & Mastalerz, M. Organic matter accumulation in the Upper Devonian Duvernay Formation, Western Canada Sedimentary Basin, from sequence stratigraphic analysis and geochemical proxies. Sed. Geol. 376, 185–203 (2018).CAS 

Google Scholar 
Hildred, G. V., Ratcliffe, K. T., Wright, A. M., Zaitlin, B. A. & Wray, D. S. Chemostratigraphic applications to low-accommodation fluvial incised-valley settings: An example from the Lower Mannville Formation of Alberta, Canada. J. Sedim. Res. 80(11), 1032–1045 (2010).
Google Scholar 
Wedepohl, K. H. Environmental influences on the chemical composition of shales and clays. Phys. Chem. Earth 8, 307–333 (1971).ADS 

Google Scholar 
Pearce, T. J., Martin, J. H., Cooper, D. & Wray, D. S. Chemostratigraphy of upper carboniferous (Pennsylvanian) sequences from the Southern North Sea (United Kingdom). Application of Modern Stratigraphic Techniques: Theory and Case Histories. SEPM Spec. Publ. 94, 109–127 (2010).
Google Scholar 
Adachi, M., Yamamoto, K. & Sugisaki, R. Hydrothermal chert and associated siliceous rocks from the northern Pacific their geological significance as indication of ocean ridge activity. Sed. Geol. 47(1–2), 125–148 (1986).CAS 

Google Scholar 
Abercrombie, H. J., Hutcheon, I. E., Bloch, J. D. & Caritat, P. D. Silica activity and the smectite-illite reaction. Geology 22(6), 539–542 (1994).ADS 
CAS 

Google Scholar 
McLennan, S. M. Weathering and global denudation. J. Geol. 101(2), 295–303 (1993).ADS 

Google Scholar 
Nesbitt, H. W. & Young, G. M. Formation and diagenesis of weathering profiles. J. Geol. 97(2), 129–147 (1989).ADS 
CAS 

Google Scholar 
Fedo, C. M., Wayne Nesbitt, H. & Young, G. M. Unraveling the effects of potassium metasomatism in sedimentary rocks and paleosols, with implications for paleoweathering conditions and provenance. Geology 23(10), 921–924 (1995).ADS 
CAS 

Google Scholar 
Nesbitt, H. W. & Young, G. M. Early Proterozoic climates and plate motions inferred from major element chemistry of lutites. Nature 299, 715–717 (1982).ADS 
CAS 

Google Scholar 
von Eynatten, H., Barceló-Vidal, C. & Pawlowsky-Glahn, V. Modelling compositional change: The example of chemical weathering of granitoid rocks. Math. Geol. 35(3), 231–251 (2003).
Google Scholar 
Clark, R. N. & Rencz, A. N. Spectroscopy of rocks and minerals, and principles of spectroscopy. Manual Remote Sens. 3(11), 3–58 (1999).
Google Scholar 
Kump, L. R. & Arthur, M. A. Interpreting carbon-isotope excursions: Carbonates and organic matter. Chem. Geol. 161(1–3), 181–198 (1999).ADS 
CAS 

Google Scholar 
Holmden, C. et al. Carbon isotope chemostratigraphy of Frasnian sequences in Western Canada. Saskatchewan Geol. Surv. Summary Investig. 1, 1–6 (2006).
Google Scholar 
Pisarzowska, A. & Racki, G. Isotopic chemostratigraphy across the Early-Middle Frasnian transition (Late Devonian) on the South Polish carbonate shelf: A reference for the global punctata Event. Chem. Geol. 334, 199–220 (2012).ADS 
CAS 

Google Scholar 
Racki, G. & Bultynck, P. Conodont biostratigraphy of the Middle to Upper Devonian boundary Beds in the Kielce area of the Holy Cross Mts. Acta Geol. Pol. 44, 1–25 (1993).
Google Scholar 
Ziegler and Sandberg. The Late Devonian standard conodont zonation CFS, Cour. Forschungsinst. Senckenberg, 121 (1990).Klapper, G., The Montagne Noire Frasnian (Upper Devonian) conodont succession. In McMillan, N.J., et al., eds., Devonian of the world, Volume III: Canadian Society of Petroleum Geologists Memoir 14, p. 449–468 (1988).Jiao, X. et al. Mixed biogenic and hydrothermal quartz in Permian lacustrine shale of Santanghu Basin, NW China: Implications for penecontemporaneous transformation of silica minerals. Int. J. Earth Sci. 107(6), 1989–2009 (2018).CAS 

Google Scholar 
Peltonen, C., Marcussen, Ø., Bjørlykke, K. & Jahren, J. Clay mineral diagenesis and quartz cementation in mudstones: The effects of smectite to illite reaction on rock properties. Mar. Pet. Geol. 26(6), 887–898 (2009).CAS 

Google Scholar 
Pearce, T. J., Besly, B. M., Wray, D. S. & Wright, D. K. Chemostratigraphy: A method to improve interwell correlation in barren sequences—a case study using onshore Duckmantian/Stephanian sequences (West Midlands, UK). Sed. Geol. 124(1–4), 197–220 (1999).CAS 

Google Scholar 
Calvert, S. E. & Pedersen, T. F. Chapter fourteen elemental proxies for palaeoclimatic and palaeoceanographic variability in marine sediments: Interpretation and application. Dev. Mar. Geol. 1, 567–644 (2007).
Google Scholar 
Perri, F., Cirrincione, R., Critelli, S., Mazzoleni, P. & Pappalardo, A. Clay mineral assemblages and sandstone compositions of the Mesozoic Longobucco Group, northeastern Calabria: Implications for burial history and diagenetic evolution. Int. Geol. Rev. 50(12), 1116–1131 (2008).
Google Scholar 
Johnson, J. G., Klapper, G. & Sandberg, C. A. Devonian eustatic fluctuations in Euramerica. Geol. Soc. Am. Bull. 96(5), 567–587 (1985).ADS 

Google Scholar 
Warme, J. E. & Sandberg, C. A. Alamo megabreccia: Record of a Late Devonian impact in southern Nevada. GSA Today 6(1), 1–7 (1996).
Google Scholar 
Ernst, R. E., Rodygin, S. A. & Grinev, O. M. Age correlation of Large Igneous Provinces with Devonian biotic crises. Glob. Planet. Change 185, 103097 (2020).
Google Scholar 
Schiffbauer, J. D. et al. Decoupling biogeochemical records, extinction, and environmental change during the Cambrian SPICE event. Sci. Adv. 3(3), e1602158 (2017).ADS 

Google Scholar 
Duller, R. A., Armitage, J. J., Manners, H. R., Grimes, S. & Jones, T. D. Delayed sedimentary response to abrupt climate change at the Paleocene-Eocene boundary, northern Spain. Geology 47(2), 159–162 (2019).ADS 
CAS 

Google Scholar  More

We used a custom designed biologging tag package with onboard video to describe a 3D high-resolution pursuit between a solitary sailfish and an individual small tuna in open water, representing the first time such an interaction has been documented. The sailfish was tagged at 09:53 on 18 October 2019, and the tag package remained attached to the sailfish for 67 h. However, analyses here are limited to the 24 h period in which the predation event took place (19 October–20 October; ~ 14 h after tagging and ~ 9 h after post-release recovery18) because this coincides with the time period the video camera was recording during daylight hours (on at 0600, off at 1800, sunrise and sunset, respectively) enabling us to ground-truth acceleration signals. Biologging data and accompanying video show the sailfish performing oscillatory dives between the surface and depths of 40–50 m during daylight hours. At night, fewer dives were performed and the sailfish generally remained within the top 10–20 m of the water column (Fig. 1a), leading to a greater range of temperatures experienced during the day (day 20.9–27.9 °C; night 26.5–28.2 °C). Due to the temperature dependence of the estimated active metabolic rate (AMRE), the cooler temperatures at depth led to a reduced AMRE during daylight hours (212.9 ± 89.1 mgO2 kg−1 h−1) compared to night (224.7 ± 44.4 mgO2 kg−1 h−1). Additionally, AMRE initially increases with depth due to increased swim speeds during diving (Fig. 1b), until the thermocline is reached in the 30–40 m depth bin, at which point AMRE decreases with further increased depth (Fig. 1b, c). However, due to thermal inertia of large-bodied fishes19,20,21, it is possible that the sailfish’s body retained heat during the short (14.7 ± 1.7 min) excursions below the thermocline and did not drop to ambient temperature. As such, the metabolic rate calculated at depth may be underestimated with the temperature correction performed here. For example, during the dive in which the predation event occurred (Fig. 1; Table 1), if body temperature was assumed equivalent to surface temperature throughout the dive, estimated metabolic rates would increase by 18% compared to if the metabolic rates were temperature corrected according to the tag’s external temperature reading (Table S2). Yet, because the majority ( > 90%) of time over the 24 h was spent above the thermocline, the temperature correction has little impact on the daily calculated AMRE and subsequent energy expenditure ( More

Eckert, K. L. & Eckert, A. E. An Atlas of Sea Turtle Nesting Habitat for the Wider Caribbean Region Revised edn. WIDECAST Technical Report no. 19 (2019).Guzmán-Hernández, V. Informe técnico 2017 del Programa de conservación de tortugas marinas en Laguna de Términos, Campeche, México. Contiene información de: 1. CPCTM Isla Aguada y 2. Reseña estatal (APFFLT/RPCyGM/CONANP/SEMARNAT, 2018).Guzmán-Hernández, V. et al. Recovery of green turtle populations and their interactions with coastal dune as a baseline for an integral ecological restoration. Acta Bot. Mex. 129, 1. https://doi.org/10.21829/abm129.2022.1954 (2022).Article 

Google Scholar 
Garduño-Andrade, M., Guzmán, V., Miranda, E., Briseño-Dueñas, R. & Abreu-Grobois, F. A. Increases in hawksbill turtle (Eretmochelys imbricata) nestings in the Yucatán Peninsula, Mexico, 1977–1996: Data in support of successful conservation?. Chelonian Conserv. Biol. 3, 286–295 (1999).
Google Scholar 
Guzmán-Hernández, V., Escanero-Figueroa, G. & Márquez, R. Programa tortuguero en el Centro Regional de Investigación Pesquera de Ciudad del Carmen, Campeche: Retrospectiva, avances y perspectivas. In Tortugas Marinas (eds Márquez, R. & Garduño-Dionate, M.) (Instituto Nacional de Pesca, 2014).
Google Scholar 
Guzmán-Hernández, V. Informe técnico 2021 del Programa de conservación de tortugas marinas en Laguna de Términos, Campeche, México. Contiene información de: 1. CPCTM Isla Aguada y 2. Reseña estatal (APFFLT/RPCyGM/CONANP/SEMARNAT, 2022).Troëng, S. & Rankin, E. Long-term conservation efforts contribute to positive green turtle Chelonia mydas nesting trend at Tortuguero, Costa Rica. Biol. Conserv. 121, 111–116 (2005).Article 

Google Scholar 
Lira-Reyes, D. et al. Informe final de la temporada de anidación 2021 del programa para la conservación de la tortuga marina en El Cuyo, Yucatán e Isla Holbox, Quintana Roo (Dirección General de Vida Silvestre de SEMARNAT/Pronatura Península de Yucatán, 2021).Piacenza, S. E., Balazs, G. H., Hargrove, S. K., Richards, P. M. & Heppell, S. S. Trends and variability in demographic indicators of a recovering population of green sea turtles Chelonia mydas. Endanger. Species Res. 31, 103–117 (2016).Article 

Google Scholar 
Iles, T. C. & Beverton, R. J. H. Stock, recruitment and moderating processes in flatfish. J. Sea Res. 39, 41–55 (1998).Article 
ADS 

Google Scholar 
Heppell, S. S., Crowder, L. B., Crouse, D. T., Epperly, S. P. & Frazer, N. B. Population models for Atlantic Loggerheads: Past, present, and future. In Loggerhead Sea Turtles (eds Bolten, A. & Witherington, B.) 255–273 (Smithsonian Institution Press, 2003).
Google Scholar 
Beverton, R. J. H. & Holt, S. J. On the dynamics of exploited fish populations (Fishery Investigation Series II, 1957).Ricker, W. E. Stock and recruitment. J. Fish. Res. Board Can. 11, 559–623 (1954).Article 

Google Scholar 
Cushing, D. H. The dependence of recruitment on parent stock in different groups of fishes. ICES J. Mar. Sci. 33, 340–362. https://doi.org/10.1093/icesjms/33.3.340 (1971).Article 

Google Scholar 
Iles, T. C. A review of stock-recruitment relationships with reference to flatfish populations. Neth. J. Sea Res. 32, 399–420 (1994).Article 

Google Scholar 
Hilborn, R. & Walters, C. Quantitative Fisheries Stock Assessment: Choice, Dynamics and Uncertainty (Chapman & Hall, 1992).Book 

Google Scholar 
Subbey, S., Devine, J. A., Schaarscmidt, U. & Nash, R. D. M. Modeling and forecasting stock recruitment: Current and future perspectives. ICES J. Mar. Sci. 71, 2307–2322 (2014).Article 

Google Scholar 
del Monte-Luna, P., Villalobos-Ortíz, H. & Arreguín-Sánchez, F. Variability of sea surface temperature in the southwestern Gulf of Mexico. Cont. Shelf Res. 102, 73–79 (2015).Article 
ADS 

Google Scholar 
Van Houtan, K. S. & Halley, J. M. Long-term climate forcing in loggerhead sea turtle nesting. PLoS ONE 6, e19043. https://doi.org/10.1371/journal.pone.0019043 (2011).Article 
ADS 
CAS 

Google Scholar 
SWOT Scientific Advisory Board. The State of the World’s Sea Turtles (SWOT) Minimum Data Standards for Nesting Beach Monitoring, version 1.0 (2011).Cuevas, D., Garrido-Chávez, E. & Raymundo-Sánchez, A. Monitoreo del género Chelonia en playas con alta densidad de anidación: reporte final (PROCER/DGOR/15/2013) para la Comisión Nacional de Áreas Naturales Protegidas (Pronatura Península de Yucatán-Comisión Nacional de Áreas Naturales Protegidas, 2013).Chim-Vera, Y. A. Evaluación del esfuerzo de monitoreo del éxito de emergencia en nidos de tortuga carey y blanca en la Península de Yucatán (Instituto Tecnológico de Conkal, 2009).Guzmán-Hernández, V., Cuevas-Flores, E., García-Alvarado, P. & González-Ruíz, T. Biological monitoring of sea turtles on nesting beaches: Datasets and basic evaluations. In Successful Conservation Strategies for Sea Turtles. Achievements and Challenges (eds Lara-Uc, M. M. et al.) 41–78 (Nova Science Publishers, 2015).
Google Scholar 
Méndez-Matos, V. C., Guzmán-Hernández, V. & Rivas-Hernández, G. Dinámica poblacional de hembras de tortuga blanca (Chelonia mydas) en el estado de Campeche, México. In El uso del conocimiento de las tortugas marinas como herramienta para la restauración de sus poblaciones y hábitats asociados (eds Cuevas, E. A. et al.) 171–187 (Universidad Autónoma del Carmen, 2019).
Google Scholar 
Xavier, R., Cortez, L. P., Cuevas, E., Barata, A. & Queiroz, N. Hawksbill turtle (Eretmochelys imbricata Linnaeus 1766) and green turtle (Chelonia mydas Linnaeus 1754) nesting activity (2002–2004) at El Cuyo beach, Mexico. Amphibia-Reptilia 27, 539–547 (2006).Article 

Google Scholar 
Whiting, A. U., Chaloupka, M. & Limpus, C. J. Sampling nesting sea turtles: Impact of survey error on trend detection. Mar. Ecol. Prog. Ser. 634, 213–223 (2020).Article 
ADS 

Google Scholar 
Harrell, F. E. Regression Modeling Strategies: With Applications to Linear Models, Logistic and Ordinal Regression, and Survival Analysis. Springer Series in Statistics (Springer, 2015).Book 
MATH 

Google Scholar 
R Core Team. R: A Language and Environment for Statistical Computing (R Foundation for Statistical Computing, Vienna, 2021).Florida Fish and Wildlife Conservation Commission. Index Nesting Beach Survey Totals (1989–2020). https://myfwc.com/research/wildlife/sea-turtles/nesting/beach-survey-totals. Accessed 2022-05-04 (2022).Marjomäki, T. J. Analysis of the spawning stock-recruitment relationship of vendace (Coregonus albula (L.)) with evaluation of alternative models, additional variables, biases and errors. Ecol. Freshw. Fish 13, 46–60 (2004).Article 

Google Scholar 
Arendt, M. D., Schwenter, J. A., Witherington, B. E., Meylan, A. B. & Saba, V. S. Historical versus contemporary climate forcing on the annual nesting variability of loggerhead sea turtles in the Northwest Atlantic Ocean. PLoS ONE 8, e81097. https://doi.org/10.1371/journal.pone.0081097 (2013).Article 
ADS 
CAS 

Google Scholar 
del Monte-Luna, P., Guzmán-Hernández, V., Cuevas, E. A., Arreguín-Sánchez, F. & Lluch-Belda, D. Effect of North Atlantic climate variability on hawksbill turtles in the Southern Gulf of Mexico. J. Exp. Mar. Biol. Ecol. 412, 103–109 (2012).Article 

Google Scholar 
Howard, R., Bell, I. & Pike, D. A. Thermal tolerances of sea turtle embryos: Current understanding and future directions. Endanger. Species Res. 26, 75–86. https://doi.org/10.3354/esr00636 (2014).Article 

Google Scholar 
Broderick, A. C., Godley, B. J. & Hays, G. C. Trophic status drives interannual variability in nesting numbers of marine turtles. Proc. R. Soc. Lond. B. Biol. 268, 1481–1487 (2001).Article 
CAS 

Google Scholar 
Manzano-Sarabia, M. M. & Salinas-Zavala, C. A. Variabilidad estacional e interanual de la concentración de clorofila y temperatura superficial del mar en la región occidental del Golfo de México: 1996–2007. Interciencia 33, 628–634 (2008).
Google Scholar 
Hays, G. C. The implications of variable remigration intervals for the assessment of population size in marine turtles. J. Theor. Biol. 206, 221–227. https://doi.org/10.1006/jtbi.2000.2116 (2000).Article 
ADS 
CAS 

Google Scholar 
Saragoça-Bruno, R., Restrepo, J. A. & Valverde, R. A. Effects of El Niño Southern Oscillation and local ocean temperature on the reproductive output of green turtles (Chelonia mydas) nesting at Tortuguero, Costa Rica. Mar. Biol. 167, 1 (2020).Article 

Google Scholar 
Patrício, A. R., Hawkes, L. A., Monsinjon, J. R., Godley, B. J. & Fuentes, M. M. P. B. Climate change and marine turtles: Recent advances and future directions. Endanger. Species Res. 44, 363–395. https://doi.org/10.3354/esr01110 (2021).Article 

Google Scholar 
Valverde-Cantillo, V., Robinson, N. J. & Santidrián Tomillo, P. Influence of oceanographic conditions on nesting abundance, phenology and internesting periods of east Pacific green turtles. Mar. Biol. 166, 93. https://doi.org/10.1007/s00227-019-3541-1 (2019).Article 

Google Scholar 
Chaloupka, M. Encouraging outlook for recovery of a once severely exploited marine megaherbivore. Glob. Ecol. Biogeogr. 17, 297–304 (2008).Article 

Google Scholar 
Ariza-Gallego, M., Herrera-Carmona, J., Payán, L. & Giraldo, A. Relationship between sea surface temperature and the nesting of the Olive Ridley sea turtle Lepidochelys olivacea (Testudines: Cheloniidae) in Gorgona Island, Colombian Pacific. Rev. Biol. Trop. 68, 528–540. https://doi.org/10.15517/RBT.V68I2.38642 (2020).Article 

Google Scholar 
Ceriani, S. A., Casale, P., Brost, M., Leone, E. H. & Witherington, B. E. Conservation implications of sea turtle nesting trends: Elusive recovery of a globally important loggerhead population. Ecosphere 10, e02936. https://doi.org/10.1002/ecs2.2936 (2019).Article 

Google Scholar 
Arendt, M. D., Schwenter, J. A., Owens, D. W. & Valverde, R. A. Theoretical modeling and neritic monitoring of loggerhead Caretta caretta [Linnaeus, 1758] sea turtle sex ratio in the southeast United States do not substantiate fears of a male-limited population. Glob. Change Biol. 27, 4849–4859. https://doi.org/10.1111/gcb.15808 (2021).Article 
CAS 

Google Scholar 
Doi, T., Márquez, R., Kimoto, H. & Azeno, N. Diagnosis and Conservation of Hawksbill Turtle Population in the Cuban Archipelago. Technical Report 40 (Japan Bekko Association, 1992).Broderick, A. C. Are green turtles globally endangered?. Glob. Ecol. Biogeogr. 14, 21–26 (2006).Article 

Google Scholar 
Mazaris, A. D., Schofield, G., Gkazinou, C., Almpanidou, V. & Hays, G. C. Global sea turtle conservation successes. Sci. Adv. 3, e1600730. https://doi.org/10.1126/sciadv.1600730 (2017).Article 
ADS 

Google Scholar 
Early-Capistrán, M. M. et al. Integrating local ecological knowledge, ecological monitoring, and computer simulation to evaluate conservation outcomes. Conserv. Lett.https://doi.org/10.1111/conl.12921 (2022).Article 

Google Scholar 
Taylor, M. A., Stephenson, T. S., Chen, A. A. & Stephenson, K. A. Climate change and the Caribbean: Review and response. Caribb. Stud. 40, 169–200 (2012).Article 

Google Scholar  More

Charlesworth, B. & Charlesworth, D. Population genetics from 1966 to 2016. Heredity 118, 2–9 (2017).CAS 

Google Scholar 
Orsini, L., Vanoverbeke, J., Swillen, I., Mergeay, J. & Meester, L. Drivers of population genetic differentiation in the wild: Isolation by dispersal limitation, isolation by adaptation and isolation by colonization. Mol. Ecol. 22, 5983–5999 (2013).
Google Scholar 
Vendrami, D. L. J. et al. RAD sequencing resolves fine-scale population structure in a benthic invertebrate: Implications for understanding phenotypic plasticity. R. Soc. Open Sci. 4, 160548 (2017).ADS 

Google Scholar 
Dufresnes, C., Rodrigues, N. & Savary, R. Slow and steady wins the race: Contrasted phylogeographic signatures in two Alpine amphibians. Integr. Zool. 17, 181–190 (2021).
Google Scholar 
Frankham, R. Genetics and extinction. Biol. Conserv. 126, 131–140 (2005).
Google Scholar 
Schwartz, M. K., Luikart, G. & Waples, R. S. Genetic monitoring as a promising tool for conservation and management. Trends in Ecol. Evol. 22, 25–33 (2007).
Google Scholar 
Ottewell, K. M., Bickerton, D. C., Byrne, M. & Lowe, A. J. Bridging the gap: A genetic assessment framework for population-level threatened plant conservation prioritization and decision-making. Divers. Distrib. 22, 174–188 (2016).
Google Scholar 
Frankham, R., Bradshaw, C. J. A. & Brook, B. W. Genetics in conservation management: Revised recommendations for the 50/500 rules, red list criteria and population viability analyses. Biol. Conserv. 170, 56–63 (2014).
Google Scholar 
Hohenlohe, P. A., Funk, C. W. & Rajora, O. P. Population genomics for wildlife conservation and management. Mol. Ecol. 30, 62–82 (2020).
Google Scholar 
Angelone, S. & Holderegger, R. Population genetics suggests effectiveness of habitat connectivity measures for the European tree frog in Switzerland. J. Appl. Ecol. 46, 879–887 (2009).
Google Scholar 
Griffiths, S. M., Taylor-Cox, E. D., Behringer, D. C., Butler, M. J. IV. & Preziosi, R. F. Using genetics to inform restoration and predict resilience in declining populations of a keystone marine sponge. Biodivers. Conserv. 29, 1383–1410 (2020).
Google Scholar 
Moritz, C. Conservation units and translocations: Strategies for conserving evolutionary processes. Hereditas 130, 217–228 (1999).
Google Scholar 
Bohonak, A. J. Dispersal, gene flow, and population structure. Q. Rev. Biol. 74, 21–45 (1999).CAS 

Google Scholar 
Arguedas, N. & Parker, P. G. Seasonal migration and genetic population structure in house wrens. Condor 102, 517–528 (2000).
Google Scholar 
Quillfeldt, P. et al. Does genetic structure reflect differences in non-breeding movements? A case study in small, highly mobile seabirds. BMC Evol. Biol. 17, 160 (2017).
Google Scholar 
Charlesworth, B., Sniegowski, P. & Stephan, W. The evolutionary dynamics of repetitive DNA in eukaryotes. Nature 371, 215–220 (1994).ADS 
CAS 

Google Scholar 
Schlötterer, C. Evolutionary dynamics of microsatellite DNA. Chromosoma 109, 365–371 (2000).
Google Scholar 
Ellegren, H. Microsatellites: Simple sequences with complex evolution. Nat. Rev. Genet. 5, 435–445 (2004).CAS 

Google Scholar 
Hodel, R. G. J. et al. The report of my death was an exaggeration: A review for researchers using microsatellites in the 21st century. Appl. Plant Sci. 4, 1600025 (2016).
Google Scholar 
Dufresnes, C. & Litvinchuk, S. N. Diversity, distribution and molecular species delimitation in frogs and toads from the Eastern Palearctic. Zool. J. Linn. Soc. 195, 695–760 (2022).
Google Scholar 
Galtier, N., Nabholz, B., Glémin, S. & Hurst, G. D. D. Mitochondrial DNA as a marker of molecular diversity: A reappraisal. Mol. Evol. 18, 4541–4550 (2009).CAS 

Google Scholar 
Zink, R. M. & Barrowclough, G. Mitochondrial DNA under siege in avian phylogeography. Mol. Ecol. 17, 2107–2121 (2008).CAS 

Google Scholar 
Toews, D. P. L. & Brelsford, A. The biogeography of mitochondrial and nuclear discordance in animals. Mol. Ecol. 21, 3907–3930 (2012).CAS 

Google Scholar 
Bonnet, T., Leblois, R., Rousset, F. & Crochet, P.-A. A reassessment of explanations for discordant introgressions of mitochondrial and nuclear genomes. Evolution 71, 2140–2218 (2017).
Google Scholar 
Davey, J. W. & Blaxter, M. L. RADSeq: Next-generation population genetics. Brief Funct. Genomics 9, 416–423 (2010).CAS 

Google Scholar 
Lexer, C. et al. ‘Next generation’ biogeography: Towards understanding the drivers of species diversification and persistence. J. Biogeogr. 40, 1013–1022 (2013).
Google Scholar 
Baird, N. A. et al. Rapid SNP discovery and genetic mapping using sequenced RAD markers. PLoS ONE 3, e3376 (2008).ADS 

Google Scholar 
Peterson, B. K., Weber, J. N., Kay, E. H., Fisher, H. S. & Hoekstra, H. E. Double Digest RADseq: An inexpensive method for de novo SNP discovery and genotyping in model and non-model species. PLoS ONE 7, e37135 (2012).ADS 
CAS 

Google Scholar 
Dufresnes, C. et al. Phylogeography of a cryptic speciation continuum in Eurasian spadefoot toads (Pelobates). Mol. Ecol. 28, 3257–3270 (2019).
Google Scholar 
Sunde, J., Yıldırım, Y., Tibblin, P. & Forsman, A. Comparing the performance of microsatellites and RADseq in population genetic studies: Analysis of data for pike (Esox Lucius) and a synthesis of previous studies. Front. Genet. 11, 218 (2020).
Google Scholar 
Moussy, C. et al. Migration and dispersal patterns of bats and their influence on genetic structure. Mammal Rev. 43, 183–195 (2013).
Google Scholar 
Berthier, P., Excoffier, L. & Ruedi, M. Recurrent replacement of mtDNA and cryptic hybridization between two sibling bat species Myotis myotis and Myotis blythii. Proc. R. Soc. B: Biol. Sci. 273, 3101–3109 (2007).
Google Scholar 
Wright, P. G. R. et al. Hydrogen isotopes reveal evidence of migration of Miniopterus schreibersii in Europe. BMC Ecol. 20, 52 (2020).CAS 

Google Scholar 
Schnetter, W. Beringungsergebnisse an der Langflügelfledermaus (Miniopterus schreibersi Kühl) im Kaiserstuhl. Bonn. Zool. Beitr. 11, 150–165 (1960).
Google Scholar 
Rodrigues, L. Miniopterus schreibersii. In The Atlas of European Mammals (eds Mitchell-Jones, A. J. et al.) 154–155 (Academic Press, 1999).
Google Scholar 
Rodrigues, L., Ramos Pereira, M. J., Rainho, A. & Palmeirim, J. M. Behavioral determinants of gene flow in the bat Miniopterus schreibersii. Behav. Ecol. Sociobiol. 64, 835–843 (2010).
Google Scholar 
Rodrigues, L. & Palmeirim, J. M. Migratory behaviour of Miniopterus schreibersii (Chiroptera): When, where, and why do cave bats migrate in a Mediterranean region?. J. Zool. 274, 116–125 (2008).
Google Scholar 
Ramos Pereira, M. J., Salgueiro, P., Rodrigues, L., Coelho, M. M. & Palmeirim, J. M. Population structure of a cave-dwelling bat, Miniopterus schreibersii: Does it reflect history and social organization?. J. Hered. 100, 533–544 (2009).
Google Scholar 
Bilgin, R. et al. Circum-Mediterranean phylogeography of a bat coupled with past environmental niche modeling: A new paradigm for the recolonization of Europe?. Mol. Phylogenet. Evol. 99, 323–336 (2016).
Google Scholar 
Gürün, K. et al. A continent-scale study of the social structure and phylogeography of the bent-wing bat, Miniopterus schreibersii (Mammalia: Chiroptera), using new microsatellite data. J. Mammal. 100, 1865–1878 (2019).
Google Scholar 
Gazaryan, S., Bücs, S., Çoraman, E. Miniopterus schreibersii (errata version published in 2021). The IUCN Red List of Threatened Species 2020: e.T81633057A195856522 (2020).Miller-Butterworth, C. M., Jacobs, D. S. & Harley, E. H. Isolation and characterization of highly polymorphic microsatellite loci in Schreibers’ long-fingered bat, Miniopterus schreibersii (Chiroptera: Vespertilionidae). Mol. Ecol. Notes 2, 139–141 (2002).CAS 

Google Scholar 
Wood, R., Weyeneth, N. & Appleton, B. Development and characterisation of 20 microsatellite loci isolated from the large bent-wing bat, Miniopterus schreibersii (Chiroptera: Miniopteridae) and their cross-taxa utility in the family Miniopteridae. Mol. Ecol. Resour. 11, 675–685 (2011).
Google Scholar 
Witsenburg, F. et al. How a haemosporidian parasite of bats gets around: The genetic structure of a parasite, vector and host compared. Mol. Ecol. 24, 926–940 (2015).CAS 

Google Scholar 
Van Oosterhout, C., Hutchinson, W. F., Wills, D. P. M. & Shipley, P. micro-checker: Software for identifying and correcting genotyping errors in microsatellite data. Mol. Ecol. Notes 4, 535–538 (2004).
Google Scholar 
Parchman, T. L. et al. Genome wide association mapping of an adaptive trait in lodgepole pine. Mol. Ecol. 21, 2991–3005 (2012).CAS 

Google Scholar 
Catchen, J. M., Amores, A., Hohenlohe, P., Cresko, W. & Postlethwait, J. H. Stacks: Building and genotyping loci de novo from short-read sequences. G3 1, 171–182 (2011).CAS 

Google Scholar 
Weir, B. S. & Cockerham, C. C. Estimating F-statistics for the analyses of population structure. Evolution 38, 1358–1370 (1984).CAS 

Google Scholar 
Goudet, J. hierfstat, a package for r to compute and test hierarchical F-statistics. Mol. Ecol. Notes 5, 184–186 (2005).
Google Scholar 
Frankham, R., Ballou, J. D. & Briscoe, D. A. A Primer of Conservation Genetics (Cambridge University Press, 2004).
Google Scholar 
Weir, B. S. & Goudet, J. A unified characterization of population structure and relatedness. Genetics 206, 2085–2103 (2017).
Google Scholar 
Mantel, N. A. The detection of disease clustering and a generalized regression approach. Cancer Res. 27, 209–220 (1967).CAS 

Google Scholar 
Wright, S. Isolation by distance. Genetics 28, 114–138 (1943).CAS 

Google Scholar 
Cavalli-Sforza, L. L. & Edwards, A. W. F. Phylogenetic analysis: Model and estimation procedures. Am. J. Hum. Genet. 19, 233–257 (1967).CAS 

Google Scholar 
Paradis, E. & Schliep, K. ape 5.0: An environment for modern phylogenetics and evolutionary analyses in R. Bioinformatics 35, 526–528 (2019).CAS 

Google Scholar 
Goudet, J., Perrin, N. & Waser, P. Tests for sex-biased dispersal using bi-parentally inherited genetic markers. Mol. Ecol. 11, 1103–1114 (2002).CAS 

Google Scholar 
Frichot, E. & François, O. lea: An r package for landscape and ecological association studies. Methods Ecol. Evol. 6, 925–929 (2015).
Google Scholar 
Yannic, G. et al. High connectivity in a long-lived High-Arctic seabird, the ivory gull Pagophila eburnea. Polar Biol. 39, 221–236 (2016).
Google Scholar 
Cumer, T. et al. Landscape and climatic variations of the Quaternary shaped multiple secondary contacts among barn owls (Tyto alba) of the Western Palearctic. Mol. Biol. Evol. 39, msab343 (2022).CAS 

Google Scholar 
Boston, E. S. M., Montgomery, W. I., Hynes, R. & Prodöhl, P. A. New insights on postglacial colonization in western Europe: The phylogeography of the Leisler’s bat (Nyctalus leisleri). Proc. R. Soc. B: Biol. Sci. 282, 20142605 (2015).
Google Scholar 
Razgour, O. et al. The shaping of genetic variation in edge-of-range populations under past and future climate change. Ecol. Lett. 16, 1258–1266 (2013).
Google Scholar 
Petit, E., Balloux, F. & Goudet, J. Sex-biased dispersal in a migratory bat: A characterization using sex-specific demographic parameters. Evolution 55, 635–640 (2001).CAS 

Google Scholar 
Moussy, C. et al. Population genetic structure of serotine bats (Eptesicus serotinus) across Europe and implications for the potential spread of bat rabies (European bat lyssavirus EBLV-1). Heredity 115, 83–92 (2015).CAS 

Google Scholar 
Rossiter, S. J., Benda, P., Dietz, C., Zhang, S. & Jones, G. Rangewide phylogeography in the greater horseshoe bat inferred from microsatellites: Implications for population history, taxonomy and conservation. Mol. Ecol. 16, 4699–4714 (2007).CAS 

Google Scholar 
Dool, S. E. et al. Phylogeography and postglacial recolonization of Europe by Rhinolophus hipposideros: Evidence from multiple genetic markers. Mol. Ecol. 22, 4055–4070 (2013).CAS 

Google Scholar 
Kerth, G. et al. Communally breeding Bechstein’s bats have a stable social system that is independent from the postglacial history and location of the populations. Mol. Ecol. 17, 2368–2381 (2008).CAS 

Google Scholar 
Garrick, R. C., Banusiewicz, J. D., Burgess, S., Hyseni, C. & Symula, R. E. Extending phylogeography to account for lineage fusion. J. Biogeogr. 46, 268–278 (2019).
Google Scholar 
Burri, R. et al. The genetic basis of color-related local adaptation in a ring-like colonization around the Mediterranean. Evolution 70, 140–153 (2016).
Google Scholar 
Taberlet, P., Fumagalli, L., Wust-Saucy, A.-G. & Cosson, J.-F. Comparative phylogeography and postglacial colonization routes in Europe. Mol. Ecol. 7, 453–464 (1998).CAS 

Google Scholar 
Hewitt, G. M. Post-glacial re-colonization of European biota. Biol. J. Linn. Soc. 68, 87–112 (1999).
Google Scholar 
Gómez, A. & Lunt, D. H. Refugia within refugia: Patterns of phylogeographic concordance in the Iberian Peninsula. In Phylogeography of Southern European Refugia (eds Weiss, S. & Ferrand, N.) 155–188 (Springer, 2007).
Google Scholar 
Vonhof, M. J., Russell, A. L. & Miller-Butterworth, M. Range-wide genetic analysis of little brown bat (Myotis lucifugus) populations: Estimating the risk of spread of white-nose syndrome. PLoS ONE 10, e0128713 (2015).
Google Scholar 
Auteri, G. G. & Knowles, L. L. Decimated little brown bats show potential for adaptive change. Sci. Rep. 10, 3023 (2020).ADS 
CAS 

Google Scholar 
Gignoux-Wolfsohn, S. A. et al. Genomic signatures of selection in bats surviving white-nose syndrome. Mol. Ecol. 30, 5643–5657 (2021).
Google Scholar 
Rivers, N. M., Butlin, R. K. & Altringham, J. D. Autumn swarming behaviour of Natterer’s bats in the UK: Population size, catchment area and dispersal. Biol. Conserv. 127, 215–226 (2006).
Google Scholar 
Reis, N. R., Fregonezi, M. N., Peracchi, A. L. & Rossaneis, B. K. Metapopulation in bats of Southern Brazil. Braz. J. Biol. 72, 605–609 (2012).CAS 

Google Scholar 
Humphrey, S. R. & Oli, M. K. Population dynamics and site fidelity of the cave bat, Myotis velifer, Oklahoma. J. Mammal. 96, 946–956 (2015).
Google Scholar 
Jeffries, D. L. et al. Comparing RADseq and microsatellites to infer complex phylogeographic patterns, an empirical perspective in the Crucian carp, Carassius carassius. L. Mol. Ecol. 25, 2997–3018 (2016).
Google Scholar 
Hodel, R. G. J. et al. Adding loci improves phylogeographic resolution in red mangroves despite increased missing data: Comparing microsatellites and RAD-Seq and investigating loci filtering. Sci. Rep. 7, 17598 (2017).ADS 

Google Scholar 
Lemopoulos, A. et al. Comparing RADseq and microsatellites for estimating genetic diversity and relatedness—Implications for brown trout conservation. Ecol. Evol. 9, 2106–2120 (2019).
Google Scholar 
Zimmerman, S. J., Aldridge, C. L. & Oyler-McCance, S. J. An empirical comparison of population genetic analyses using microsatellite and SNP data for a species of conservation concern. BMC Genom. 21, 382 (2020).CAS 

Google Scholar 
Hale, M. L., Burg, T. M. & Steeves, T. E. Sampling for microsatellite-based population genetic studies: 25 to 30 individuals per population is enough to accurately estimate allele frequencies. PLoS ONE 7, e45170 (2012).ADS 
CAS 

Google Scholar 
Quetglas, J., Gonzalez, F. & Paz, O. Estudian la extraña mortandad de miles de murcielago de cuevas. Quercus 203, 50 (2003).
Google Scholar 
Negredo, A. et al. Discovery of an ebolavirus-like filovirus in Europe. PLoS Pathog. 7, e1002304 (2011).CAS 

Google Scholar 
Reed, D. H. & Frankham, R. Correlation between fitness and genetic diversity. Conserv. Biol. 17, 230–237 (2003).
Google Scholar 
Alcalde, J. T., Artácoz, A. & Meijide, F. Recuperación de la colonia de Miniopterus schreibersii de la cueva de Cueva de Ágreda (Soria). Barbastella 5, 32–35 (2012).
Google Scholar 
Kemenesi, G. et al. Re-emergence of Lloviu virus in Miniopterus schreibersii bats, Hungary, 2016. Emerg. Microbes Infect. 7, 66 (2018).
Google Scholar 
Kemenesi, et al. Isolation of infectious Lloviu virus from Schreiber’s bats in Hungary. Nat. Commun. 13, 1706 (2022).ADS 
CAS 

Google Scholar 
Stoffel, C. et al. Genetic consequences of population expansions and contractions in the common hippopotamus (Hippopotamus amphibius) since the late Pleistocene. Mol. Ecol. 24, 2507–2520 (2015).
Google Scholar  More

Quality controls for the validation of MICODE protocolTo validate the MICODE experimental approach in the version of fecal batch of the human proximal colon, we chose to monitor and check some parameters as quality controls (QC) related to metabolites and microbes at the end of fermentations, and in comparison, to the baseline. QCs adopted were; (i) the Firmicutes/Bacteroidetes ratio (F/B), which is related to health and disease11, was maintained at a low level, confirming the capacity to simulate a healthy in vivo condition for 24 h. (ii) The presence of Archea (e.g., Methanobrevibacter smithii and Methanosphaera stadtmanae), which are pretty sensible to oxygen content12, was retained from the baseline to the end point in each vessel and repetition, indicating that the environmental conditions were strictly maintained. (iii) Good’s rarity index of alpha biodiversity remained similar during time of fermentation (p  > 0.05), indicating enough support to the growth of rare species. (iv) Observed OTUs richness index scored approximately 400 OTUs at the end point. (v) The paradigm of prebiotics was confirmed when the positive control (FOS) was applied on MICODE; high probiotic and SCFAs increases and limitation of enteropathogens. (vi) Each GC/MS analysis had quantified some stool-related compounds (urea, 1-propanol, and butylated hydroxy toluene), that ranged across the complete chromatogram and were adsorbed at the same retention times.Changes in bacterial alpha and beta diversitiesThe microbiota diversity indices were analyzed to study the impact of HPBA on microbial population, to assess population’s stability during fermentation, and to compare its microbiota to that of other bioreactors (Figure S1). The baseline of value was compared to the endpoints of fermentation of different treatments. It is undisputable that a part of the effect of reduction in richness (Observed OTUs) was derived by the passage from in vivo to in vitro condition, but the focus must be set on the different trend that other alpha diversity indices had. For example, abundance (Chao 1) for HBPA was significantly higher at the end of fermentation (p  0.05) and HPBA (p  0.05), while oppositely, FOS decreased in evenness (p  > 0.05) and raised in dominance (p  0.05). Among these, 31 variables were significant and their Log2 fold changes in respect to the baseline were compared by post-hoc test (Table 1). The 41 OTUs selected were those that recorded shifts after fermentation and that from literature are susceptible to the effect of prebiotic or fiber substrates. We have included even three OTUs of Archea relative to QC of the experiments (previously discussed).Table 1 Abundances (% ± S.D.) and changes in phylum taxa (Log2 F/C) after 24 h in vitro fecal batch culture fermentations from healthy donors and administrated with HBPA, HB, and FOS as the substrates, and also including a blank control.Full size tableThe first group of OTUs included beneficial or commensal bacteria that usually respond to prebiotics. In this group, three Bifidobacterium were picked showing increases on the substrates and reduction on the blank control. HB and HBPA fostered Bif. bifidum, but just the latter did it significantly, making this taxon grew up to the 3.30% of relative abundance (p  More

Data collectionPeer-reviewed articles published up to Dec 31, 2021 were searched through the Web of Science (http://webofknowledge.com), Google scholar (http://scholar.google.com), and the China National Knowledge Infrastructure (CNKI, http://www.cnki.net). Here we employed a combination of the following search terms: “(tree biomass OR aboveground biomass OR plant biomass OR plant productivity) and (allometric biomass equation OR allometric model OR productivity model OR biomass equation OR biomass model)”. To avoid potential selection bias and duplicates, we conducted a cross-check between the references of relevant articles, which resulted in the selection of 729 relevant articles from the thousands of the appearing articles initially. Subsequently, eligible articles were selected using the following criteria: (1) Allometric models built for specific species with confirmed locations without disturbances were selected, generalized species, large-scales (e.g., province or nation), as well as recently disturbed tree models were excluded. (2) The method employed to develop the model was destructive harvesting and weighing, with at least twenty sample trees, were selected; articles were excluded that did not include measurements and used less than twenty sample trees. (3) The model forms were W = a*Db and LnW = a + b*Ln(D) or W = a*(D2H)b and LnW = a + b*Ln(D2H), where W is the aboveground biomass, and D is the diameter at breast height, H is the tree height, were selected. Consequently, we excluded articles with other variables and other forms of models. Finally, 426 articles remained from the original 729 (Supplementary Fig. S1).We then distilled data from the articles for the following variables: (1) Allometric models, in the form of W = a*Db and LnW = a + b*Ln(D), W = a*(D2H)b and LnW = a + b*Ln(D2H) including the parameters a, b in the D range and H range. (2) Tree species corresponding to the models, including families, genera, and species. (3) Location data, including longitude, latitude, and study sites. (4) Climate data, including mean annual temperature (MAT, °C) and mean annual precipitation (MAP, mm) of the tree species location. (5) Terrain data, including slope and aspect. (6) Soil data, including soil organic carbon (SOC), clay, and soil type.Since not all articles provided the location, climate, soil, and terrain data of the studies, we estimated the missing data as follows, (1) we supplemented the longitude and latitude with the study location using Google Earth. (2) We extracted the missing climate data by using geographic coordinates from WorldClim version 2.0 (http://worldclim.org/current)16. (3) We obtained the shuttle radar topographic mission DEM data with 30 m resolution from NASA, and used SAGA-GIS software to derive various terrain data from the DEM such as altitude, slope, and aspect17, 18. (4) The missing soil data was derived from the Regridded Harmonized World Soil Database v1.219. In particular, we established the soil type according to Soil Taxonomy to increase the accuracy of the analysis and prediction. Furthermore, if the experiments were performed at multiple sites in one study, they were treated as independent observations. In light of above criteria, 817 allometric models in the form of W = a*Db or LnW = a + b*Ln(D) and 612 allometric models in the form of W = a*(D2H)b or LnW = a + b*Ln(D2H) were collected from the 426 articles.Allometric modelThe relationship between the diameter and aboveground biomass was in the form of the power function20:$$begin{array}{c}Wi=atimes D{i}^{b},end{array}$$
(1)
where Wi is the dry mass of the ith tree (kg), Di is diameter at breast height (cm), and a and b are the parameters of the model.$$Wi=atimes (D{i}^{2}Hi{)}^{b},$$
(2)
where Wi is the dry mass of the ith tree (kg), Di is diameter at breast height (cm), Hi is the tree height (cm), and a and b are the parameters of the model.However, a heteroscedasticity exists when directly fitting the tree biomass. The logarithmic transformation of Eq. (1) or Eq. (2), is convenient to facilitate model fitting and deal with heterocedasticity21. The logarithmic transformation allometric model:$$begin{array}{c}Lnleft(Wiright)=a+btimes Lnleft(Diright),end{array}$$
(3)

was used in this function, where a (Eq. 3) represents Ln(a) (Eq. 1), and b (Eq. 3) is the same as b (Eq. 1), respectively.$$begin{array}{c}Lnleft(Wiright)=a+btimes Lnleft(D{i}^{2}Hright),end{array}$$
(4)
was used in this function, where a (Eq. 4) represents Ln(a) (Eq. 2), and b (Eq. 4) is the same as b (Eq. 2), respectively. To unify the models, we transformed the collected Eqs. (1) to (3) and Eqs. (2) to (4).Data analysisTo establish the relationship between variables with parameters a and b for making a parameter prediction on a global scale, Random Forest (RF) (an example of a machine learning model) was employed, which consists of an ensemble of randomized classification and regression trees (CART)21. In short, the RF will generate a number of trees and aggregate these to provide a single prediction. In regression problems the prediction is the average of the individual tree outputs, whereas in classification the trees vote by majority on the correct classification22, 23. Generated trees called ntree are based on a bootstrapped 2/3 sample of the original data to decrease correlations by choosing different training sets in the RF modeling process15. In addition to this normal bagging function, the best split at each node of the tree was searched only among a randomly selected subset (mtry) of predictors24. The tree growing procedure is performed recursively until the size of the node reaches a minimum, k, which is parameterized by the user. For the rest of the original data, RF provides a believable error estimation using the data called Out-Of-Bag (OOB), which is employed to obtain a running unbiased estimate of the classification error as trees are added to the forest15.Predictive variable selectionThe variables included stand factors such as density, family, and diameters, as well as non-stand factors such as MAT, MAP, and SOC. Considering that the prediction was on a global scale, the first step was to exclude the factors that it was not possible to completely extract. Next, we selected variables through the following22: (1) the RF classifier was initially applied using all of the predictor variables, and variable importance was used to rank them based on the mean decrease in accuracy. (2) Removing the least important variables by the variable importance ranking, (3) the training data were then partitioned five-fold for cross-validation and the error rates for each of the five cross-validation partitions were aggregated into a mean error rate, and 20 replicates of the five-fold CV were performed25.By means of the above, eleven variables, including family, genus, species, MAT, MAP, altitude, aspect, SOC, slope, clay, and soil type, were remained to predict parameters. Since the combinations of variables were different, five combinations were performed to make predictions from the eleven variables above. Among the five combinations, each were used by RF to predict and select via the model evaluation index VaR explained and the mean of squared residual (Supplementary Table S1).Optimization of Random Forest parametersRF depends primarily on three parameters that are set by users. (1) ntree, the number of trees in the forest. (2) nodesize, the minimum number of data points in each terminal node. (3) mtry, the number of features tried at each node. To obtain the optimization of RF parameters, we set ntree = 1000, 2000, 3000 and the selection criterion was that ntree was small enough to maximize computational efficiency as well as produced stable OOB error25. As for nodesize, we used 3, 5, 7, and 5 as the default for regression RF, given that the mtry value always is always one third of the number of variables. Here we also set the mtry values (ranging from 2 to 4), which were tested, and we accessed the OOB error rates from 50 replicates for each mtry value25. The primary tuning parameter above were optimized, as well as each combination of the three RF parameters through a grided search, which were used to predict and set RF parameters according to the predictive effect of each combination (Supplementary Table S2).All above data analysis were conducted in R 4.0.326. And the output is the spatial pattern of allometric model parameters at 0.5° resolution.Predicted parameter validationFurther to assess the accuracy of the predicted parameters, we applied them to estimate the AGB at six sites. And the actual AGB of the sites had been obtained via destructive sampling from 209 plots, which were located in Hubei, Liaoning, Gansu, Hebei and Heilongjiang provinces, and Inner Mongolia autonomous region from 2009 to 201327 (Table 1). First, we selected the sample trees according the dominant, average and inferior tree outside the plot. Then the sample trees were felled as carefully as possible and tree height (H), tree diameter in the breast (DBH) and live crown length were recorded. To divide trees into several sub-samples, including branches, leaves, stem wood and stem bark, all of the branches were removed and leaves were picked. Besides, stem was divided into 1 m sections and bark of the stem was removed. Finally, all sub-samples of aboveground part of trees were oven-dried at 80 °C until a constant weight was reached and the sum of all the sun-samples weight was the actual AGB. Through the above process, 249 actual AGB data were obtained. Meanwhile, the predicted parameters of the models together with the DBH and H estimated the predicted AGB. The actual AGB data of 249 sample trees were compared with the predicted AGB by making fitting curves between them in R to show the availability of predicted parameters according root mean square error (RMSE) and R2.Table 1 The basic features of the sampling sites.Full size tableThe experimental research and field studies on plants in this study, including the collection of plant material, complied with the relevant institutional, national, and international guidelines and legislation. And we ensured that we have permission for the plant sampling, all of the steps were allowed in our study for the plant research. In addition, plant identification in this study was conducted by X.Z according to World Plants (https://www.worldplants.de) in the herbarium of School of Forestry & Landscape of Architecture, Anhui Agricultural University, and the voucher specimen of all plant material has been deposited in a publicly available herbarium. More