Ecology

Friedlingstein, P. et al. Earth Syst. Sci. Data 12, 3269–3340 (2020).Article 

Google Scholar 
LeBauer, D. S. & Treseder, K. K. Ecology 89, 371–379 (2008).Article 
PubMed 

Google Scholar 
Vitousek, P. M. & Howarth, R. W. Biogeochemistry 13, 87–115 (1991).Article 

Google Scholar 
Barker, W. et al. Nature https://doi.org/10.1038/s41586-022-05502-6 (2022).Article 

Google Scholar 
Sprent, J. I. Legume Nodulation: A Global Perspective (Wiley-Blackwell, 2009).
Google Scholar 
Gei, M. et al. Nature Ecol. Evol. 2, 1104–1111 (2018).Article 
PubMed 

Google Scholar 
Peng, J. et al. Glob. Biogeochem. Cycles 34, e2019GB006296 (2020).Article 

Google Scholar 
Batterman, S. A. et al. Nature 502, 224–227 (2013).Article 
PubMed 

Google Scholar 
Taylor, B. N. & Menge, D. N. L. Nature Plants 4, 655–661 (2018).
Google Scholar 
McCulloch, L. A. & Porder, S. New Phytol. 231, 1734–1745 (2021).Article 
PubMed 

Google Scholar 
Sheffer, E., Batterman, S. A., Levin, S. A. & Hedin, L. O. Nature Plants 1, 15182 (2015).Article 

Google Scholar 
Barron, A. R., Purves, D. W. & Hedin, L. O. Oecologia 165, 511–520 (2011).Article 
PubMed 

Google Scholar 
Adams, M. A., Turnbull, T. L., Sprent, J. I. & Buchmann, N. Proc. Natl Acad. Sci. USA 113, 4098–4103 (2016).Article 
PubMed 

Google Scholar 
Taylor, B. N. & Ostrowsky, L. R. J. Trop. Ecol. 35, 270–279 (2019).Article 

Google Scholar 
Sprent, J. I. New Phytol. 174, 11–25 (2007).Article 
PubMed 

Google Scholar  More

SamplingSediment samples were obtained from the Kap København Formation in North Greenland (82° 24′ 00″ N 22° 12′ 00″ W) in the summers of 2006, 2012 and 2016 (see Supplementary Table 3.1.1). Sampled material consisted of organic-rich permafrost and dry permafrost. Prior to sampling, profiles were cleaned to expose fresh material. Samples were hereafter collected vertically from the slope of the hills either using a 10 cm diameter diamond headed drill bit or cutting out ~40 × 40 × 40 cm blocks. Sediments were kept frozen in the field and during transportation to the lab facility in Copenhagen. Disposable gloves and scalpels were used and changed between each sample to avoid cross-contamination. In a controlled laboratory environment, the cores and blocks were further sub-sampled for material taking only the inner part of sediment cores, leaving 1.5–2 cm between the inner core and the surface that provided a subsample of approximately 6–10 g. Subsequently, all samples were stored at temperatures below −22 °C.We sampled organic-rich sediment by taking samples and biological replicates across the three stratigraphic units B1, B2 and B3, spanning 5 different sites, site: 50 (B3), 69 (B2), 74a (B1), 74b (B1) and 119 (B3). Each biological replicate from each unit at each site was further sampled in different sublayers (numbered L0–L4, Source Data 1, sheet 1).Absolute age datingIn 2014, Be and Al oxide targets from 8× 1 kg quartz-rich sand samples collected at modern depths ranging from 3 to 21 m below stream cut terraces were analysed by accelerator mass spectrometry and the cosmogenic isotope concentrations interpreted as maximum ages using a simple burial dating approach1 (26Al:10Be versus normalized 10Be). The 26Al and 10Be isotopes were produced by cosmic ray interactions with exposed quartz in regolith and bedrock surfaces in the mountains above Kap København prior to deposition. We assume that the 26Al:10Be was uniform and steady for long time periods in the upper few metres of these gradually eroding palaeo-surfaces. Once eroded by streams and hillslope processes, the quartz sand was deposited in sandy braided stream sediment, deltaic distributary systems, or the near-shore environment and remained effectively shielded from cosmic ray nucleons buried (many tens of metres) under sediment, intermittent ice shelf or ice sheet cover, and—at least during interglacials—the marine water column until final emergence. The simple burial dating approach assumes that the sand grains experienced only one burial event. If multiple burial events separated by periods of re-exposure occurred, then the starting 26Al:10Be before the last burial event would be less than the initial production ratio (6.75 to 7.42, see discussion below) owing to the relatively faster decay of 26Al during burial, and therefore the calculated burial age would be a maximum limiting age. Multiple burial events can be caused by shielding by thick glacier ice in the source area, or by sediment storage in the catchment prior to final deposition. These shielding events mean that the 26Al:10Be is lower, and therefore a calculated burial age assuming the initial production ratio would overestimate the final burial duration. We also consider that once buried, the sand grains may have been exposed to secondary cosmogenic muons (their depth would be too great for submarine nucleonic production). As sedimentation rates in these glaciated near-shore environments are relatively rapid, we show that even the muonic production would be negligible (see Supplemental Information). However, once the marine sediments emerged above sea level, in-situ production by both nucleogenic and muogenic production could alter the 26Al:10Be. The 26Al versus 10Be isochron plot reveals this complex burial history (Supplementary Information, section 3) and the concentration versus depth composite profiles for both 26Al and 10Be reveal that the shallowest samples may have been exposed during a period of time (~15,000 years ago) that is consistent with deglaciation in the area (Supplemental Information). While we interpret the individual simple burial age of all samples as a maximum limiting age of deposition of the Kap København Formation Member B, we recommend using the three most deeply shielded samples in a single depth profile to minimize the effect of post-depositional production. We then calculate a convolved probability distribution age for these three samples (KK06A, B and C). However, this calculation depends on the 26Al:10Be production ratio we use (that is, between 6.75 and 7.42) and on whether we adjust for erosion in the catchment. So, we repeat the convolved probability distribution function age for the lowest and highest production ratio and zero to maximum possible erosion rate, to obtain the minimum and maximum limiting age range at 1σ confidence (Supplementary Information, section 3). Taking the midpoint between the negative and positive 3σ confidence limits, we obtain a maximum burial age of 2.70 ± 0.46 Myr. This age is also supported by the position of those three samples on the isochron plot, which suggests the true age may not be significantly different that this maximum limiting age.Thermal ageThe extent of thermal degradation of the Kap København DNA was compared to the DNA from the Krestovka Mammoth molar. Published kinetic parameters for DNA degradation64 were used to calculate the relative rate difference over a given interval of the long-term temperature record and to quantify the offset from the reference temperature of 10 °C, thus estimating the thermal age in years at 10 °C for each sample (Supplementary Information, section 4). The mean annual air temperature (MAT) for the the Kap København sediment was taken from Funder et al. (2001)6 and for the Krestovka Mammoth the MAT was calculated using temperature data from the Cerskij Weather Station (WMO no. 251230) 68.80° N 161.28° E, 32 m from the International Research Institute Data Library (https://iri.columbia.edu/) (Supplementary Table 4.4.1).We did not correct for seasonal fluctuation for the thermal age calculation of the Kap København sediments or from the Krestovka Mammoth. We do provide theoretical average fragment length for four different thermal scenarios for the DNA in the Kap København sediments (Supplementary Table 4.4.2). A correction in the thermal age calculation was applied for altitude using the environmental lapse rate (6.49 °C km−1). We scaled the long-term temperature model of Hansen et al. (2013)65 to local estimates of current MATs by a scaling factor sufficient to account for the estimates of the local temperature decline at the last glacial maximum and then estimated the integrated rate using an activation energy (Ea) of 127 kJ mol−1 (ref. 64).Mineralogic compositionThe minerals in each of the Kap København sediment samples were identified using X-ray diffraction and their proportions were quantified using Rietveld refinement. The samples were homogenized by grinding ~1 g of sediment with ethanol for 10 min in a McCrone Mill. The samples were dried at 60 °C and added corundum (CR-1, Baikowski) as the internal standard to a final concentration of 20.0 wt%. Diffractograms were collected using a Bruker D8 Advance (Θ–Θ geometry) and the LynxEye detector (opening 2.71°), with Cu Kα1,2 radiation (1.54 Å; 40 kV, 40 mA) using a Ni-filter with thickness of 0.2 mm on the diffracted beam and a beam knife set at 3 mm. We scanned from 5–90° 2θ with a step size of 0.1° and a step time of 4 s while the sample was spun at 20 rpm. The opening of the divergence slit was 0.3° and of the antiscatter slit 3°. Primary and secondary Soller slits had an opening of 2.5° and the opening of the detector window was 2.71°. For the Rietveld analysis, we used the Profex interface for the BGMN software66,67. The instrumental parameters and peak broadening were determined by the fundamental parameters ray-tracing procedure68. A detailed description of identification of clay minerals can be found in the supporting information.AdsorptionWe used pure or purified minerals for adsorption studies. The minerals used and treatments for purifying them are listed in Supplementary Table 4.2.6. The purity of minerals was checked using X-ray diffraction with the same instrumental parameters and procedures as listed in the above section i.e., mineralogical composition. Notes on the origin, purification and impurities can be found in the Supplementary Information section 4. We used artificial seawater69 and salmon sperm DNA (low molecular weight, lyophilized powder, Sigma Aldrich) as a model for eDNA adsorption. A known amount of mineral powder was mixed with seawater and sonicated in an ultrasonic bath for 15 min. The DNA stock was then added to the suspension to reach a final concentration between 20–800 μg ml−1. The suspensions were equilibrated on a rotary shaker for 4 h. The samples were then centrifuged and the DNA concentration in the supernatant determined with UV spectrometry (Biophotometer, Eppendorf), with both positive and negative controls. All measurements were done in triplicates, and we made five to eight DNA concentrations per mineral. We used Langmuir and Freundlich equations to fit the model to the experimental isotherm and to obtain adsorption capacity of a mineral at a given equilibrium concentration.PollenThe pollen samples were extracted using the modified Grischuk protocol adopted in the Geological Institute of the Russian Academy of Science which utilizes sodium pyrophosphate and hydrofluoric acid70. Slides prepared from 6 samples were scanned at 400× magnification with a Motic BA 400 compound microscope and photographed using a Moticam 2300 camera. Pollen percentages were calculated as a proportion of the total palynomorphs including the unidentified grains. Only 4 of the 6 samples yielded terrestrial pollen counts ≥50. In these, the total palynomorphs identified ranged from 225 to 71 (mean = 170.25; median = 192.5). Identifications were made using several published keys71,72. The pollen diagram was initially compiled using Tilia version 1.5.1273 but replotted for this study using Psimpoll 4.1074.DNA recoveryFor recovery calculation, we saturated mineral surfaces with DNA. For this, we used the same protocol as for the determination of adsorption isotherms with an added step to remove DNA not adsorbed but only trapped in the interstitial pores of wet paste. This step was important because interstitial DNA would increase the amount of apparently adsorbed DNA and overestimate the recovery. To remove trapped DNA after adsorption, we redispersed the minerals in seawater. The process of redispersing the wet paste in seawater, ultracentrifugation and removal of supernatant lasted less than 2.5 min. After the second centrifugation, the wet pastes were kept frozen until extraction. We used the same extraction protocol as for the Kap København sediments. After the extraction, the DNA concentration was again determined using UV spectrometry.MetagenomesA total of 41 samples were extracted for DNA75 and converted to 65 dual-indexed Illumina sequencing libraries (including 13 negative extraction- and library controls)30. 34 libraries were thereafter subjected to ddPCR using a QX200 AutoDG Droplet Digital PCR System (Bio-Rad) following manufacturer’s protocol. Assays for ddPCR include a P7 index primer (5′-AGCAGAAGACGGCATAC-3′) (900nM), gene-targeting primer (900 nM), and a gene-targeting probe (250nM). We screened for Viridiplantae psbD (primer: 5′-TCATAATTGGACGTTGAACC-3′, probe: 5′-(FAM)ACTCCCATCATATGAAA(BHQ1)-3′) and Poaceae psbA (primer: 5′-CTCACAACTTCCCTCTAGAC-3′, probe 5′-(HEX)AGCTGCTGTTGAAGTTC(BHQ1)-3′). Additionally, 34 of the 65 libraries were enriched using targeted capture enrichment, for mammalian mitochondrial DNA using the PaleoChip Arctic1.0 bait-set31 and all libraries were hereafter sequenced on an Illumina HiSeq 4000 80 bp PE or a NovaSeq 6000 100 bp PE. We sequenced a total of 16,882,114,068 reads which, after low complexity filtering (Dust = 1), quality trimming (q ≥ 25), duplicate removal and filtering for reads longer than 29 bp (only paired read mates for NovaSeq data) resulted in 2,873,998,429 reads that were parsed for further downstream analysis. We next estimated kmer similarity between all samples using simka32 (setting heuristic count for max number of reads (-max-reads 0) and a kmer size of 31 (-kmer-size 31)), and performed a principal component analysis (PCA) on the obtained distance matrix (see Supplementary Information, ‘DNA’). We hereafter parsed all QC reads through HOLI33 for taxonomic assignment. To increase resolution and sensitivity of our taxonomic assignment, we supplemented the RefSeq (92 excluding bacteria) and the nucleotide database (NCBI) with a recently published Arctic-boreal plant database (PhyloNorway) and Arctic animal database34 as well as searched the NCBI SRA for 139 genomes of boreal animal taxa (March 2020) of which 16 partial-full genomes were found and added (Source Data 1, sheet 4) and used the GTDB microbial database version 95 as decoy. All alignments were hereafter merged using samtools and sorted using gz-sort (v. 1). Cytosine deamination frequencies were then estimated using the newly developed metaDMG, by first finding the lowest common ancestor across all possible alignments for each read and then calculating damage patterns for each taxonomic level36 (Supplementary Information, section 6). In parallel, we computed the mean read length as well as number of reads per taxonomic node (Supplementary Information, section 6). Our analysis of the DNA damage across all taxonomic levels pointed to a minimum filter for all samples at all taxonomic levels with a D-max ≥ 25% and a likelihood ratio (λ-LR) ≥ 1.5. This ensured that only taxa showing ancient DNA characteristics were parsed for downstream profiling and analysis and resulted in no taxa within any controls being found (Supplementary Information, section 6).Marine eukaryotic metagenomeWe sought to identify marine eukaryotes by first taxonomically labelling all quality-controlled reads as Eukaryota, Archaea, Bacteria or Virus using Kraken 276 with the parameters ‘–confidence 0.5 –minimum-hit-groups 3’ combined with an extra filtering step that only kept those reads with root-to-leaf score >0.25. For the initial Kraken 2 search, we used a coarse database created by the taxdb-integration workflow (https://github.com/aMG-tk/taxdb-integration) covering all domains of life and including a genomic database of marine planktonic eukaryotes63 that contain 683 metagenome-assembled genomes (MAGs) and 30 single-cell genomes (SAGs) from Tara Oceans77, following the naming convention in Delmont et al.63, we will refer to them as SMAGs. Reads labelled as root, unclassified, archaea, bacteria and virus were refined through a second Kraken 2 labelling step using a high-resolution database containing archaea, bacteria and virus created by the taxdb-integration workflow. We used the same Kraken 2 parameters and filtering thresholds as the initial search. Both Kraken 2 databases were built with parameters optimized for the study read length (–kmer-len 25 –minimizer-len 23 –minimizer-spaces 4).Reads labelled as eukaryota, root and unclassified were hereafter mapped with Bowtie278 against the SMAGs. We used MarkDuplicates from Picard (https://github.com/broadinstitute/picard) to remove duplicates and then we calculated the mapping statistics for each SMAG in the BAM files with the filterBAM program (https://github.com/aMG-tk/bam-filter). We furthermore estimated the postmortem damage of the filtered BAM files with the Bayesian methods in metaDMG and selected those SMAGs with a D-max ≥ 0.25 and a fit quality (λ-LR) higher than 1.5. The SMAGs with fewer than 500 reads mapped, a mean read average nucleotide identity (ANI) of less than than 93% and a breadth of coverage ratio and coverage evenness of less than 0.75 were removed. We followed a data-driven approach to select the mean read ANI threshold, where we explored the variation of mapped reads as a function of the mean read ANI values from 90% to 100% and identified the elbow point in the curve (Supplementary Fig. 6.11.1). We used anvi’o79 in manual mode to plot the mapping and damage results using the SMAGs phylogenomic tree inferred by Delmont et al. as reference. We used the oceanic signal of Delmont et al. as a proxy to the contemporary distribution of the SMAGs in each ocean and sea (Fig. 5 and Supplementary Information, section 6).Comparison of DNA, macrofossil and pollenTo allow comparison between records in DNA, macrofossil and pollen, the taxonomy was harmonized following the Pan Arctic Flora checklist43 and NCBI. For example, since Bennike (1990)18, Potamogeton has been split into Potamogeton and Stuckenia, Polygonym has been split to Polygonum and Bistorta, and Saxifraga was split to Saxifraga and Micranthes, whereas others have been merged, such as Melandrium with Silene40. Plant families have changed names—for instance, Gramineae is now called Poaceae and Scrophulariaceae has been re-circumscribed to exclude Plantaginaceae and Orobancheae80. We then classified the taxa into the following: category 1 all identical genus recorded by DNA and macrofossils or pollen, category 2 genera recorded by DNA also found by macrofossils or pollen including genus contained within family level classifications, category 3 taxa only recorded by DNA, category 4 taxa only recorded by macrofossils or pollen (Source Data 1).Phylogenetic placementWe sought to phylogenetically place the set of ancient taxa with the most abundant number of reads assigned, and with a sufficient number of reference sequences to build a phylogeny. These taxa include reads mapped to the chloroplast genomes of the flora genera Salix, Populus and Betula, and to the mitochondrial genomes of the fauna families Elephantidae, Cricetidae, Leporidae, as well as the subfamilies Capreolinae and Anserinae. Although the evolution of the chloroplast genome is somewhat less stable than that of the plant mitochondrial genome, it has a faster rate of evolution, and is non-recombining, and hence is more likely to contain more informative sites for our analysis than the plant mitochondria81. Like the mitochondrial genome, the chloroplast genome also has a high copy number, so that we would expect a high number of sedimentary reads mapping to it.For each of these taxa, we downloaded a representative set of either whole chloroplast or whole mitochondrial genome fasta sequences from NCBI Genbank82, including a single representative sequence from a recently diverged outgroup. For the Betula genus, we also included three chloroplast genomes from the PhyloNorway database34,83. We changed all ambiguous bases in the fasta files to N. We used MAFFT84 to align each of these sets of reference sequences, and inspected multiple sequence alignments in NCBI MSAViewer to confirm quality85. We trimmed mitochondrial alignments with insufficient quality due to highly variable control regions for Leporidae, Cricetidae and Anserinae by removing the d-loop in MegaX86.The BEAST suite49 was used with default parameters to create ultrametric phylogenetic trees for each of the five sets of taxa from the multiple sequence alignments (MSAs) of reference sequences, which were converted from Nexus to Newick format in Figtree (https://github.com/rambaut/figtree). We then passed the multiple sequence alignments to the Python module AlignIO from BioPython87 to create a reference consensus fasta sequence for each set of taxa. Furthermore, we used SNPSites88 to create a vcf file from each of the MSAs. Since SNPSites outputs a slightly different format for missing data than needed for downstream analysis, we used a custom R script to modify the vcf format appropriately. We also filtered out non-biallelic SNPs.From the damage filtered ngsLCA output, we extracted all readIDs uniquely classified to reference sequences within these respective taxa or assigned to any common ancestor inside the taxonomic group and converted these back to fastq files using seqtk (https://github.com/lh3/seqtk). We merged reads from all sites and layers to create a single read set for each respective taxon. Next, since these extracted reads were mapped against a reference database including multiple sequences from each taxon, the output files were not on the same coordinate system. To circumvent this issue and avoid mapping bias, we re-mapped each read set to the consensus sequence generated above for that taxon using bwa89 with ancient DNA parameters (bwa aln -n 0.001). We converted these reads to bam files, removed unmapped reads, and filtered for mapping quality  > 25 using samtools90. This produced 103,042, 39,306, 91,272, 182 and 129 reads for Salix, Populus, Betula, Elephantidae and Capreolinae, respectively.We next used pathPhynder62, a phylogenetic placement algorithm that identifies informative markers on a phylogeny from a reference panel, evaluates SNPs in the ancient sample overlapping these markers, and traverses the tree to place the ancient sample according to its derived and ancestral SNPs on each branch. We used the transversions-only filter to avoid errors due to deamination, except for Betula, Salix and Populus in which we used no filter due to sufficiently high coverage. Last, we investigated the pathPhynder output in each taxon set to determine the phylogenetic placement of our ancient samples (see Supplementary Information for discussion on phylogenetic placement).Based on the analysis described above we further investigated the phylogenetic placement within the genus Mammut, or mastodons. To avoid mapping reference biases in the downstream results, we first built a consensus sequence from all comparative mitochondrial genomes used in said analysis and mapped the reads identified in ngsLCA as Elephantidae to the consensus sequence. Consensus sequences were constructed by first aligning all sequences of interest using MAFFT84 and taking a majority rule consensus base in Geneious v2020.0.5 (https://www.geneious.com). We performed three analyses for phylogenetic placement of our sequence: (1) Comparison against a single representative from each Elephantidae species including the sea cow (Dugong dugon) as outgroup, (2) Comparison against a single representative from each Elephantidae species, and (3) Comparison against all published mastodon mitochondrial genomes including the Asian elephant as outgroup.For each of these analyses we first built a new reference tree using BEAST v1.10.4 (ref. 47) and repeated the previously described pathPhynder steps, with the exception that the pathPhynder tree path analysis for the Mammut SNPs was based on transitions and transversions, not restricting to only transversions due to low coverage.
Mammut americanum
We confirmed the phylogenetic placement of our sequence using a selection of Elephantidae mitochondrial reference sequences, GTR+G, strict clock, a birth-death substitution model, and ran the MCMC chain for 20,000,000 runs, sampling every 20,000 steps. Convergence was assessed using Tracer91 v1.7.2 and an effective sample size (ESS)  > 200. To determine the approximate age of our recovered mastodon mitogenome we performed a molecular dating analysis with BEAST47 v1.10.4. We used two separate approaches when dating our mastodon mitogenome, as demonstrated in a recent publication92. First, we determined the age of our sequence by comparing against a dataset of radiocarbon-dated specimens (n = 13) only. Secondly, we estimated the age of our sequence including both molecularly (n = 22) and radiocarbon-dated (n = 13) specimens using the molecular dates previously determined92. We utilized the same BEAST parameters as Karpinski et al.92 and set the age of our sample with a gamma distribution (5% quantile: 8.72 × 104, Median: 1.178 × 106; 95% quantile: 5.093 × 106; initial value: 74,900; shape: 1; scale: 1,700,000). In short, we specified a substitution model of GTR+G4, a strict clock, constant population size, and ran the Markov Chain Monte Carlo chain for 50,000,000 runs, sampling every 50,000 steps. Convergence of the run was again determined using Tracer.Molecular dating methodsIn this section, we describe molecular dating of the ancient birch (Betula) chloroplast genome using BEAST v1.10.4 (ref. 47). In principle, the genera Betula, Populus and Salix had both sufficiently high chloroplast genome coverage (with mean depth 24.16×, 57.06× and 27.04×, respectively, although this coverage is highly uneven across the chloroplast genome) and enough reference sequences to attempt molecular dating on these samples. Notably, this is one of the reasons we included a recently diverged outgroup with a divergence time estimate in each of these phylogenetic trees. However, our Populus sample clearly contained a mixture of different species, as seen from its inconsistent placement in the pathPhynder output. In particular, there were multiple supporting SNPs to both Populus balsamifera and Populus trichocarpa, and both supporting and conflicting SNPs on branches above. Furthermore, upon inspection, our Salix sample contained a surprisingly high number of private SNPs which is inconsistent with any ancient or even modern age, especially considering the number of SNPs assigned to the edges of the phylogenetic tree leading to other Salix sequences. We are unsure what causes this inconsistency but hypothesize that our Salix sample is also a mixed sample, containing multiple Salix species that diverged from the same placement branch on the phylogenetic tree at different time periods. This is supported by looking at all the reads that cover these private SNP sites, which generally appear to be from a mixed sample, with reads containing both alternate and reference alleles present at a high proportion in many cases. Alternatively, or potentially jointly in parallel, this could be a consequence of the high number of nuclear plastid DNA sequences (NUPTs) in Salix93. Because of this, we continued with only Betula.First, we downloaded 27 complete reference Betula chloroplast genome sequences and a single Alnus chloroplast genome sequence to use as an outgroup from the NCBI Genbank repository, and supplemented this with three Betula chloroplast sequences from the PhyloNorway database generated in a recent study29, for a total of 31 reference sequences. Since chloroplast sequences are circular, downloaded sequences may not always be in the same orientation or at the same starting point as is necessary for alignment, so we used custom code (https://github.com/miwipe/KapCopenhagen) that uses an anchor string to rotate the reference sequences to the same orientation and start them all from the same point. We created a MSA of these transformed reference sequences with Mafft84 and checked the quality of our alignment by eye in Seqotron94 and NCBI MsaViewer. Next, we called a consensus sequence from this MSA using the BioAlign consensus function87 in Python, which is a majority rule consensus caller. We will use this consensus sequence to map the ancient Betula reads to, both to avoid reference bias and to get the ancient Betula sample on the same coordinates as the reference MSA.From the last common ancestor output in metaDMG36, we extracted read sets for all units, sites and levels that were uniquely classified to the taxonomic level of Betula or lower, with at a minimum sequence similarity of 90% or higher to any Betula sequence, using Seqtk95. We mapped these read sets against the consensus Betula chloroplast genome using BWA89 with ancient DNA parameters (-o 2 -n 0.001 -t 20), then removed unmapped reads, quality filtered for read quality ≥25, and sorted the resulting bam files using samtools89. For the purpose of molecular dating, it is appropriate to consider these read sets as a single sample, and so we merged the resulting bam files into one sample using samtools. We used bcftools89 to make an mpileup and call a vcf file, using options for haploidy and disabling the default calling algorithm, which can slightly biases the calls towards the reference sequence, in favour of a majority call on bases that passed the default base quality cut-off of 13. We included the default option using base alignment qualities96, which we found greatly reduced the read depths of some bases and removed spurious SNPs around indel regions. Lastly, we filtered the vcf file to include only single nucleotide variants, because we do not believe other variants such as insertions or deletions in an ancient environmental sample of this type to be of sufficiently high confidence to include in molecular dating.We downloaded the gff3 annotation file for the longest Betula reference sequence, MG386368.1, from NCBI. Using custom R code97, we parsed this file and the associated fasta to label individual sites as protein-coding regions (in which we labelled the base with its position in the codon according to the phase and strand noted in the gff3 file), RNA, or neither coding nor RNA. We extracted the coding regions and checked in Seqotron94 and R that they translated to a protein alignment well (for example, no premature stop codons), both in the reference sequence and the associated positions in the ancient sequence. Though the modern reference sequence’s coding regions translated to a high-quality protein alignment, translating the associated positions in the ancient sequence with no depth cut-off leads to premature stop codons and an overall poor quality protein alignment. On the other hand, when using a depth cut-off of 20 and replacing sites in the ancient sequence which did not meet this filter with N, we see a high-quality protein alignment (except for the N sites). We also interrogated any positions in the ancient sequence which differed from the consensus, and found that any suspicious regions (for example, with multiple SNPs clustered closely together spatially in the genome) were removed with a depth cut-off of 20. Because of this, we moved forward only with sites in both the ancient and modern samples which met a depth cut-off of at least 20 in the ancient sample, which consisted of about 30% of the total sites.Next, we parsed this annotation through the multiple sequence alignment to create partitions for BEAST47. After checking how many polymorphic and total sites were in each, we decided to use four partitions: (1) sites belonging to protein-coding positions 1 and 2, (2) coding position 3, (3) RNA, or (4) non-coding and non-RNA. To ensure that these were high confidence sites, each partition also only included those positions which had at least depth 20 in the ancient sequence and had less than 3 total gaps in the multiple sequence alignment. This gave partitions which had 11,668, 5,828, 2,690 and 29,538 sites, respectively. We used these four partitions to run BEAST47 v1.10.4, with unlinked substitution models for each partition and a strict clock, with a different relative rate for each partition. (There was insufficient information in these data to infer between-lineage rate variation from a single calibration). We assigned an age of 0 to all of the reference sequences, and used a normal distribution prior with mean 61.1 Myr and standard deviation 1.633 Myr for the root height48; standard deviation was obtained by conservatively converting the 95% HPD to z-scores. For the overall tree prior, we selected the coalescent model. The age of the ancient sequence was estimated following the overall procedures of Shapiro et al. (2011)98. To assess sensitivity to prior choice for this unknown date, we used two different priors, namely a gamma distribution metric towards a younger age (shape = 1, scale = 1.7); and a uniform prior on the range (0, 10 Myr). We also compared two different models of rate variation among sites and substitution types within each partition, namely a GTR+G with four rate categories, and base frequencies estimated from the data, and the much simpler Jukes Cantor model, which assumed no variation between substitution types nor sites within each partition. All other priors were set at their defaults. Neither rate model nor prior choice had a qualitative effect on results (Extended Data Fig. 10). We also ran the coding regions alone, since they translated correctly and are therefore highly reliable sites and found that they gave the same median and a much larger confidence interval, as expected when using fewer sites (Extended Data Fig. 10). We ran each Markov chain Monte Carlo for a total of 100 million iterations. After removing a burn-in of the first 10%, we verified convergence in Tracer91 v1.7.2 (apparent stationarity of traces, and all parameters having an Effective Sample Size  > 100). We also verified that the resulting MCC tree from TreeAnnotator47 had placed the ancient sequence phylogenetically identically to pathPhynder62 placement, which is shown in Extended Data Fig. 9. For our major results, we report the uniform ancient age prior, and the GTR+G4 model applied to each of the four partitions. The associated XML is given in Source Data 3. The 95% HPD was (2.0172,0.6786) for the age of the ancient Betula chloroplast sequence, with a median estimate of 1.323 Myr, as shown in Fig. 2.Reporting summaryFurther information on research design is available in the Nature Portfolio Reporting Summary linked to this article. More

Identification, phylogeny, and distribution of five phylaTo advance our understanding of marine sediment microbial diversity, we obtained over 30 billion paired DNA sequences from 42 marine sediment samples (coastal and deep sea) (Supplementary Data 1). From this, we reconstructed over 8000 ( >50% complete, 95%) to genes from coastal waters (Venezuela), a hypersaline pond in Carpinteria (US), sediments in Garolim Bay (Korea), and others (Supplementary Data 6 and 7). The worldwide distribution of these five phyla suggests that they have potentially overlooked ecological roles across many environments.Detection of novel protein familiesTo explore novel metabolic capabilities of these bacteria, we employed a recently described approach to identify and characterize unknown genes exclusive to uncultivated taxa17. Using this computational method, we identified 1,934 novel protein families (NPFs) and 6,893 novel singletons (NSs) in the 55 MAGs. The former can be define as families that do not show any homology in broadly used databases (including eggNOG, pfamA, pfamB, and RefSeq, see “Methods”) while the latter (NSs) are NPFs that are detected only once in each given genome or group of genomes. To determine if this novelty was specific to the five phyla or distributed across other uncultivated prokaryotic taxa, we mapped these NPFs and NSs against a comprehensive dataset of 169,642 bacterial and archaeal genomes covered in Rodriguez del Río et al.17. Using an in-house pipeline (Supplementary Fig. 4), we found that 44.6% of these NPFs and NSs are present in other uncultured taxa, highlighting the novel and undescribed metabolic repertoire that these five phyla share with other uncultured prokaryotic lineages17. Specifically, we found that these proteins are also present in Marinisomatota, Bacteroidota, and WOR-3 from publicly available genomes obtained from both marine and terrestrial environments17. When comparing the total number of NPFs per genome in the novel bacterial phyla against the genomic dataset (approximately 170,000 genomes), we found that the novel taxa described in this study have a higher than average percentage of novel proteins per genome (5.68 ± 4.89%) (p  0.7) and widespread (coverage > 0.7) within each phylum are shown in dark purple bars. The number of novel protein families with conserved neighboring genes are shown in light gray bars. c, d, Selected examples of phylogenetic trees and novel protein family genomic context marked in gray with a black outline) in Blakebacterota and Arandabacterota. The protein families are similar between these two phyla and have conserved neighboring genes, including translation initiation factor IF-3 gene (infC), large subunit ribosomal protein L20 gene (rplT), phenylalanyl-tRNA synthetase genes (pheST), cell division protein gene (zapA), phosphodiesterase gene (ymdB), methenyltetrahydrofolate cyclohydrolase gene (folD), and exodeoxyribonuclease genes (xseAB). e Phylogenetic tree and genomic context of a novel protein family uniquely distributed in Joyebacterota. The novel protein family has conserved genomic neighbors related to energy conservation (Rnf complex genes, rnfABCDEG). The phylogeny was generated using FastTree2 and numbers on the top and bottom of the branch represent the bootstrap and branch length, respectively. Source data are provided as a Source Data file.Full size imageMetabolic pathways are often encoded by ‘genome neighborhoods’ (gene clusters and/or operons)18. Therefore, we calculated the genomic context conservation of the NPFs containing three or more sequences (3773 NPFs in total) and examined the annotation of genes found in genomic proximity of the NPFs to determine their potential function. Of the inspected families, 513 (14%) had a conservation score ≥ 0.9 (see “Methods”) indicating a high degree of conserved neighboring proteins. Manual annotation of these neighboring proteins indicated they are potentially involved in sulfur reduction, energy conservation, as well as the degradation of organics such as starch, fatty acids, and amino acids (highlighted in red in Supplementary Fig. 5). For example, a NPF predominantly found in Blakebacterota is neighbored by putative menaquinone reductases (QrcABCD), a conserved complex related to energy conservation in sulfate reducing bacteria19,20,21,22. However, metabolic annotations of Blakebacterota genomes that encode QrcABCD indicate that they largely lack the key enzymes for sulfate reduction, dissimilatory sulfite reductases (DsrABC), suggesting this QrcABCD complex may be involved in other bioenergetic contexts such as linking periplasmic hydrogen and formate oxidation to the menaquinone pool22.In some instances, we found NPFs coded near genes predicted to produce key proteins in nitrogen cycling. Two of the Joyebacterota MAGs code NPF neighboring proteins with homology to hydroxylamine dehydrogenases (HAO). HAO is a key enzyme in marine nitrogen cycling that has traditionally been thought to catalyze the oxidation of hydroxylamine (NH2OH) to nitrite (NO2−) in ammonia oxidizing bacteria. Recently, it has been suggested that HAO may also convert hydroxylamine to nitric oxide (NO) as an intermediate, which is then further oxidized to nitrite by an unknown mechanism. Hydroxylamine is also known to be an intermediate in the nitrogen cycle. It is a potential precursor of nitrous oxide (N2O), a potent greenhouse gas that is a byproduct of denitrification, nitrification23,24, and anaerobic ammonium oxidation25. The presence of HAO within the genomic context of these NPFs suggests they may be involved in mediating hydroxylamine metabolism, and thus may play an important role in nitrogen cycling.A number of NPFs are colocalized with genes predicted to be involved in the utilization of organic carbon. For example, one NPF found in Blakebacterota genomes is adjacent to a peptidase (PepQ; K01271) for dipeptide degradation. Another NPF, only detected in Blakebacterota, is neighbored by long-chain acyl-CoA synthetase (FadD; K01897), a key enzyme in fatty acid degradation (Supplementary Fig. 6). In Joyebacterota, as well as in publicly available Bacteroidetes and Latescibacteria we identified an NPF that is colocalized with amylo-alpha-1,6-glucosidase (Glycoside Hydrolase Family 57), suggesting a potential role in starch degradation.We also identified NPFs that are specific and very conserved in AABM5, Blakebacterota, Orphanbacterota, Arandabacterota, and Joyebacterota (2, 39, 3, 16, and 26 respectively). These NPFs were found in at least 70% of the MAGs belonging to each phylum, and rarely present in other genomes across the tree of life. Due to their unique nature, the 86 unique NPFs could be used as marker genes for future characterizations of the novel bacteria described in this study. When examining the genomic context of the phyla-specific NPFs, we found that more than half of the NPFs (49 of 86) shared the same gene order and are next to genes predicted to be involved in various catabolic and anabolic processes. For example, an NPF in Joyebacterota MAGs is adjacent to an Rnf complex26, which is important for energy conservation in numerous organisms21 (Fig. 2e). Also, two different NPFs in Blakebacterota and Arandabacterota MAGs were located next to tRNA synthesis genes (Fig. 2c, d). Additional phyla-specific NPFs were colocalized with genes predicted to be involved in other important processes, including peptidoglycan biosynthesis (Supplementary Fig. 6a), F-type ATPase (Supplementary Fig. 6b), acyl-CoA dehydrogenase, elements for transportation, sulfur assimilation (Supplementary Fig. 6c), and others (Supplementary Fig. 6d).Metabolic potential of the novel bacterial phylaIn addition to NPF-based analyses, we compared the predicted proteins in the novel lineages to a variety of databases and gene phylogenies to understand their metabolism (see “Methods”). The distribution of key metabolic proteins based on presence/absence of protein families (using MEBS: see methods) in the 61 MAGs is largely consistent with their phylogeny (Fig. 1a). Below, we detail the predicted metabolism of each novel bacterial phyla based on these analyses (Supplementary Fig. 5 and Supplementary Data 8 and 9, see details in Supplementary Information).JoyebacterotaJoyebacterota is composed of 20 MAGs predominantly reconstructed from hydrothermal vent sediments (blue, lower right side in the phylogeny shown in Fig. 1a). Metabolic inference suggests that these bacteria are obligate anaerobes encoding extracellular carbohydrate-active enzymes (CAZymes) with the potential to degrade pectate or pectin, photosynthetically fixed carbon in marine diatoms, macrophytes27, and terrestrial plants28. Furthermore, Joyebacterota seems to be involved in the sulfur cycle. Seven Joyebacterota MAGs encode sulfide:quinone oxidoreductases (SQR). Phylogenetic analysis indicate these SQR belong to the membrane-bound type I and III29. Interestingly, these SQR type I sequences are closely related to those sequences mostly found in terrestrial environments, e.g., freshwater, soil, and hot spring, while SQR-III  have been previously suggested to play a key role maintaining the sulfide homeostasis or bioenergetics in deep-sea sediments30. The presence of these pathways highlight the potential adaptation of Joyebacterota to several environments, contributing to recycling of carbon and sulfur.BlakebacterotaThe Blakebacterota phylum is composed of 11 MAGs predominantly reconstructed from the surface layer of GB sediments (0–6 cm). In this environment, temperatures range from 25 to 29 °C, CH4 measures 0.4–0.8 mM, CO2 reaches up to 10 mM, and SO42− concentrations are high (up to 28 mM)30. Metabolic inference using MEBS31 suggests Blakebacterota play an important role in N and S cycles. These findings were supported by the presence of key enzymes in these cycles. For example, we identified a nitrous oxide reductase in Blakebacterota, the only known enzyme to catalyze the reduction of nitrous oxide to nitrogen gas. This reaction acts as a sink for nitrous oxide, and thus is an important removal mechanism for this potent greenhouse gas. In addition to nitrogen cycling, we identified key genes involved in sulfur cycling in Blakebacterota. Six of the MAGs possess genes that code for SQR with sulfate or nitrous oxide as the final electron accepter. In addition, seven of the MAGs contain genes for thiosulfate dehydrogenase (doxD), which may convert thiosulfate to tetrathionate. Finally, one MAG is predicted to produce dimethyl sulfide (DMS) under oxic conditions via methanethiol S-methyltransferase (MddA) from methylate L-methionine or methanethiol (MeSH). Thus, these bacteria may play important roles in a variety of intermediate steps in nitrogen and sulfur cycling.ArandabacterotaLike Joyebacterota, Arandabacterota were largely recovered from shallow (2–14 cm) GB and deep (26–38 cm) BS sediments. This phylum contains 11 MAGs that are predicted to be anaerobic polysulfide and elemental sulfur reducers. They may mediate sulfur reduction via sulfhydrogenases (HydGB), which results in the production of sulfide32,33. Thus, Arandabacterota may contribute to sulfur cycling in marine sediments. Arandabacterota also code distinct hydrogenases, [NiFe] 3c and 4g types, (Fig. 3) for H2 oxidation. In addition, Arandabacterota may reduce nitrite via periplasmic dissimilatory nitrite reductases (NrfAH) present in Meg22_24_Bin_129, BHB10-38_Bin_9, and SY70-4-3_Bin_59. This mechanism for energy conservation is more efficient than polysulfide and elemental sulfur reduction. Therefore, they are likely to use sulfur species as electron donors in the absence of nitrite.Fig. 3: Maximum likelihood phylogenetic tree of NiFe hydrogenases from the novel phyla.The majority of NiFe hydrogenases identified from the five phyla in this study are highlighted in the gray background. Most hydrogenases are types 4g and 3c. Starred branches denote the minor NiFe hydrogenases identified in this study. Bootstrap values ≥ 80 are shown in circles. Source data are provided as a Source Data file.Full size imageOrphanbacterotaOrphanbacterota is composed of seven MAGs that were mostly obtained from the BS, and appear to be metabolically versatile, facultative aerobes. The BS has an average water depth of 18 m and is strongly influenced by anthropogenic activities in China, mainly the terrestrial input of nutrients and organic matter34. Orphanbacterota code a diversity of CAZymes for the degradation of complex carbohydrates. We identified genes coding for extracellular glycoside hydrolase family 16 (GH16), which may be involved in the degradation of laminarin, releasing glucose and oligosaccharides35. Six Orphanbacterota genomes also contain genes predicted to produce extracellular peptidases belonging to family M28 and S8, which are nonspecific peptidases (Supplementary Fig. 7 and Supplementary Data 10–14). The released amino acids could be taken up via ABC transporters coded by these bacteria.Consistent with their recovery from shallow sediment habitats (Supplementary Data 1), Orphanbacterota have a diverse repertoire of terminal cytochrome oxidase genes (Supplementary Data 9) suggesting they are capable of surviving in a range of oxygen concentrations. Based on the presence of isocitrate lyase and malate synthase, they may use the glyoxylate cycle for carbohydrate synthesis when sugar is not available, or use simple two-carbon compounds for energy conservation36,37. They also appear capable of reducing nitrate to nitrite via periplasmic nitrate reductases (NapAB)38. Moreover, they could reduce nitrate via the membrane-bound nitrate reductase for energy conservation and reducing nitrous oxide.One Orphanbacterota genome (M3-44_Bin_119) has genes predicted to mediate sulfate/sulfite reduction, including DsrABC, QmoABC, and membrane bound Rnf complexes (Supplementary Fig. 8a, b and Supplementary Data 8 and 9). Another Orphanbacterota (LQ108M_Bin_12) is predicted to contain diverse metabolic pathways, including MmdA for DMS production, SQR for sulfide oxidation, the Rnf complex for energy conservation21 or detoxification (Supplementary Fig. 8c), and sulfhydrogenases (HydABDG) for H2 oxidation. In addition to energy conservation and detoxification, sulfide oxidation is important for preventing the loss of sulfur through H2S volatilization. This is predicted to be an important process in sulfur-rich sediments, where large quantities of the self-produced H2S are produced during heterotrophic growth29.AABM5AABM5 (12 genomes, 7 obtained in this study) is an understudied bacterial group that has largely been recovered from shallow (4–12 cm) sediments in GB and deep (44–62 cm) sediments in BS. Despite the distinct environments where they have been found, genomes within this phylum have several shared metabolic abilities. In contrast to the strict anaerobic lifestyle that was previously reported in a subgroup within AABM5 (candidate division LCP–89)12, we predict they are facultative anaerobes. In support of this, we identified cytochrome c oxidase (CtaDCEF) and cytochrome bd ubiquinol oxidase (CydAB) for aerobic respiration39. In addition, we identified DsrABC in nine genomes (Supplementary Fig. 8 and Supplementary Data 15), indicating these organisms can potentially reduce sulfate/sulfite for energy conservation. Several AABM5 genomes are predicted to use H2 as an electron donor due to the presence of type 3c [NiFe] hydrogenase (MvhADG) (Fig. 3, Supplementary Fig. 9, and Supplementary Data 8 and 9). The metabolic versatility in this phylum better explains their global distribution.Ecological significance of the new phylaThese previously overlooked bacterial phyla appear to be involved in key biogeochemical processes in marine sediments, namely sulfur and nitrogen cycling, and the degradation of organic carbon. However, we did not find any evidence for complete autotrophic metabolisms (Wood-Ljungdahl pathway, Calvin–Benson–Bassham, reductive tricarboxylic acid, 3-hydroxypropionate bicycle, 3-hydroxypropionate-4-hydroxybutyrate, and dicarboxylate-4-hydroxybutyrate cycles) in any of these bacteria. Instead, they have a variety of pathways for the utilization of organic compounds as detailed above. These novel bacteria phyla (all except Blakebacterota) have the potential to degrade the algal glycan laminarin, one of the most important complex carbon compounds in the ocean40. These novel phyla encode extracellular laminarinases that specifically cleave the laminarin into more readily degradable sugars, e.g., glucose and oligosaccharide (Supplementary Fig. 7 and Supplementary Data 10–12). Laminarin glycan is produced in the surface ocean by microalgae that sequester CO2 as an important carbon sink in the oceans41. This is a key process of the global carbon cycle, and most studies have focused on understanding aerobic laminarin-degrading bacteria in the surface oceans41,42. Recently, it has been shown that laminarin plays a prominent role in oceanic carbon export and energy flow to higher trophic levels and the deep ocean40, yet the organisms responsible for laminarin degradation under anoxic conditions are unknown. The discovery of  these novel bacterial phyla opens new doors for future studies exploring laminarin degradation in the deep sea. In addition, most of them contain genes predicted to code for sulfatases. Blakebacterota, Orphanbacterota, Arandabacterota, and Joyebacterota code for arylsulfatase, mainly arylsulfatase A, for desulfation of galactosyl moiety of sulfatide. They also code choline sulfatase, iduronate 2-sulfatase and some uncharacterized sulfatases for different types of substrates43. This suggests they are capable of cleaving organic sulfate ester bonds as a source of sulfur and organic carbon on the ocean floor.Many metabolic processes identified here, including pathways for polysaccharide degradation, sulfur, and nitrogen metabolism are often incomplete (Fig. 4). This may be due to the incompleteness of these genomes, or it suggests that these processes occur via metabolic handoffs within the community. Some of the phyla are capable of mediating a variety of sulfur and nitrogen redox reactions (Fig. 4a, b). For example, four phyla code DsrABC, suggesting they play an overlooked role in inorganic matter degradation in marine sediments through sulfate reduction. The resultant sulfide may be reoxidized to sulfur intermediates and organic sulfur compounds by these newly described bacteria. Four phyla (Blakebacterota, Orphanbacterota, Arandabacterota, and Joyebacterota) code an SQR for producing elemental sulfur from sulfide. Methanethiol S-methyltransferase (MddA) is predicted to be produced by individual MAGs Blakebacterota (M3-38_Bin_215) and Orphanbacterota (LQ108M_Bin_12) for the production of DMS from methionine44. DMS is important in climate regulation and sulfur cycling in marine environments45,46, though little is known about the fate or production of DMS in anoxic environments like marine sediments. As detailed above, Blakebacterota contains genes for the conversion of thiosulfate to tetrathionate. Four phyla (AABM5, Orphanbacterota, Arandabacterota, and Joyebacterota) are predicted to disproportionate thiosulfate to sulfite via thiosulfate/3-mercaptopyruvate sulfurtransferase. Thus, we suspect these bacteria may be capable of mediating intermediate sulfur species in anoxic environments. These results provide a predictive framework for future physiological studiesto confirm our genomic-based predictions.Fig. 4: Genomic-based predictions of the potential metabolic role of the novel bacterial phyla.Key steps in the (a) sulfur and (b) nitrogen cycles predicted in the five bacterial phyla. Compounds (in gray triangle frames) were arranged according to the standard Gibbs free energy of formation of each sulfur or nitrogen compound (values next to the compound taken from Caspi et al.93). Star, square, triangle, pentagon, and diamond shapes correspond to AABM5, Blakebacterota, Orphanbacterota, Arandabacterota and Joyebacterota, respectively. Colored shapes represent the presence of genes in a given pathway. Fully colored shapes indicate the presence of genes in over 50% of the phyla. Half colored shapes signify that less than 50% of the phyla code for those genes. Uncolored shapes indicate presence of genes in only one MAG. Note that only pathways encoded in at least one MAG are shown. The red dotted line indicates the assimilatory process. The blue soild line indicates the confirmed pathway with phylogeny of key genes. c Phylogenetic tree and genomic context of a novel protein family (NPF) next to putative menaquinone reductase complex genes (qrcABCD) found in Blakebacterota and Orphanbacterota. d Phylogenetic tree and genomic context of a NPF next to hydroxylamine oxidoreductase genes (hao) in Joyebacterota.Full size imageIn addition to potential roles in sulfur cycling, the phyla described here may play key roles in nitrogen processes, for example several MAGs contain genes that code predicted hydroxylamine dehydrogenase proteins (HAO, confirmed by different databases)47,48. HAO is a precursor of nitrous oxide (N2O), a potent greenhouse gas and ozone destructing agent in the atmosphere. Marine N2O stems from nitrification and denitrification processes which depend on organic matter cycling and dissolved oxygen. Since hydroxylamine is a precursor of N2O, deciphering the organisms that can mediate the formation of N2O has important implications for Earth’s climate49. In addition, three phyla (AABM5, Blakebacterota, and Orphanbacterota) code for periplasmic and/or transmembrane nitrate reductase, and two phyla (AABM5 and Arandabacterota) are predicted to reduce nitrite via dissimilatory nitrite reductase.In recent years, there have been large advances in the exploration of novel microbial diversity. Genomic data has provided crucial insights into the ecological roles and biology of these new microbes. The recovery of bacterial genomes belonging to five overlooked, globally distributed phyla with considerably novel protein composition reminds us there is much to be learned about the microbial world. The identification of NPFs provides targets for future studies to elucidate the ecophysiology of these organisms. The presence of genes for organic carbon degradation and sulfur and nitrogen cycling in these new bacteria suggests they contribute to a variety of key processes in marine sediments. Thus, the addition of these bacterial genomes to ecosystem models will likely transform our understanding of how microbial communities drive carbon degradation and nutrient cycling in the oceans. More

Institute for Health Metrics and Evaluation, University of Washington, Seattle, WA, USARobert C. Reiner Jr., Catherine A. Welgan, Christopher E. Troeger, Mathew M. Baumann, Aniruddha Deshpande, Brigette F. Blacker, Molly K. Miller-Petrie, Lucas Earl, Daniel C. Casey, Aubrey J. Cook, Farah Daoud, Nicole Davis Weaver, Samath Dhamminda Dharmaratne, Laura Dwyer-Lindgren, Valery L. Feigin, Joseph Jon Frostad, Kimberly B. Johnson, Alice Lazzar-Atwood, Kate E. LeGrand, Stephen S. Lim, Paulina A. Lindstedt, Laurie B. Marczak, Benjamin K. Mayala, Ali H. Mokdad, Jonathan F. Mosser, Chrisopher J. L. Murray, QuynhAnh P. Nguyen, David M. Pigott, Puja C. Rao, David L. Smith, Emma Elizabeth Spurlock & Simon I. HayDepartment of Health Metrics Sciences, School of Medicine, University of Washington, Seattle, WA, USARobert C. Reiner Jr., Samath Dhamminda Dharmaratne, Laura Dwyer-Lindgren, Stephen S. Lim, Ali H. Mokdad, Chrisopher J. L. Murray, David M. Pigott, Benn Sartorius, David L. Smith & Simon I. HayMalaria Atlas Project, University of Oxford, Oxford, UKDaniel J. Weiss & Susan Fred RumishaImperial College London, London, UKSamir BhattDepartment of Laboratory Medicine, Karolinska University Hospital, Huddinge, SwedenHassan AbolhassaniResearch Center for Immunodeficiencies, Tehran University of Medical Sciences, Tehran, IranHassan Abolhassani & Nima RezaeiDepartment of Public Health, Debre Berhan University, Debre Berhan, EthiopiaAkine Eshete AbosetugnDepartment of Clinical Sciences, University of Sharjah, Sharjah, United Arab EmiratesEman Abu-GharbiehPopulation Health Sciences, King’s College London, London, EnglandVictor AdekanmbiCentre of Excellence for Epidemiological Modelling and Analysis, Stellenbosch University, Stellenbosch, South AfricaOlatunji O. AdetokunbohDepartment of Global Health, Stellenbosch University, Cape Town, South AfricaOlatunji O. AdetokunbohDepartment of Epidemiology and Biostatistics, Qom University of Medical Sciences, Qom, IranMohammad AghaaliFaculty of Medicine and Public Health, Jenderal Soedirman University, Purwokerto, IndonesiaBudi AjiMayo Evidence-based Practice Center, Mayo Clinic Foundation for Medical Education and Research, Rochester, MN, USAFares AlahdabJohn T. Milliken Department of Internal Medicine, Washington University in St. Louis, St. Louis, MO, USAZiyad Al-AlyClinical Epidemiology Center, Department of Veterans Affairs, St Louis, MO, USAZiyad Al-AlyInstitute of Health Research, University of Health and Allied Sciences, Ho, GhanaRobert Kaba AlhassanDepartment of Information Systems, College of Economics and Political Science, Sultan Qaboos University, Muscat, OmanSaqib AliInfectious and Tropical Disease Research Center, Hormozgan University of Medical Sciences, Bandar Abbas, IranHesam AlizadeDepartment of Health Policy and Management, Kuwait University, Safat, KuwaitSyed Mohamed AljunidInternational Centre for Casemix and Clinical Coding, National University of Malaysia, Bandar Tun Razak, MalaysiaSyed Mohamed AljunidDepartment of Epidemiology, Arak University of Medical Sciences, Arak, IranAmir Almasi-Hashiani, Rahmatollah Moradzadeh & Maryam ZamanianMedical Research Center, Jazan University, Jazan, Saudi ArabiaHesham M. Al-MekhlafiDepartment of Parasitology, Sana’a University, Sana’a, YemenHesham M. Al-MekhlafiPediatric Intensive Care Unit, King Saud University, Riyadh, Saudi ArabiaKhalid A. Altirkawi & Mohamad-Hani TemsahResearch Group in Health Economics, University of Cartagena, Cartagena, ColombiaNelson Alvis-GuzmanResearch Group in Hospital Management and Health Policies, ALZAK Foundation, Cartagena, ColombiaNelson Alvis-GuzmanSchool of Medicine, University of Adelaide, Adelaide, SA, AustraliaAzmeraw T. AmareCollege of Medicine and Health Science, Bahir Dar University, Bahir Dar, EthiopiaAzmeraw T. AmareHealth Services Management Department, Arak University of Medical Sciences, Arak, IranSaeed AminiMaternal and Child Wellbeing, African Population and Health Research Center, Nairobi, KenyaDickson A. AmugsiPharmacy Department, Carol Davila University of Medicine and Pharmacy, Bucharest, RomaniaRobert AncuceanuCardiology Department, Carol Davila University of Medicine and Pharmacy, Bucharest, RomaniaCatalina Liliana AndreiResearch Center for Evidence Based Medicine, Tabriz University of Medical Sciences, Tabriz, IranFereshteh AnsariRazi Vaccine and Serum Research Institute, Agricultural Research, Education, and Extension Organization (AREEO), Tehran, IranFereshteh AnsariDepartment of Parasitology, Mazandaran University of Medical Sciences, Sari, IranDavood AnvariDepartment of Parasitology, Iranshahr University of Medical Sciences, Iranshahr, IranDavood AnvariDepartment of Sociology and Social Work, Kwame Nkrumah University of Science and Technology, Kumasi, GhanaSeth Christopher Yaw AppiahCenter for International Health, Ludwig Maximilians University, Munich, GermanySeth Christopher Yaw AppiahHealth Management and Economics Research Center, Iran University of Medical Sciences, Tehran, IranJalal Arabloo & Ahmad GhashghaeeDepartment of Public Health, Birmingham City University, Birmingham, UKOlatunde AremuFaculty of Nursing, Philadelphia University, Amman, JordanMaha Moh’d Wahbi AtoutSchool of Business, University of Leicester, Leicester, UKMarcel AusloosDepartment of Statistics and Econometrics, Bucharest University of Economic Studies, Bucharest, RomaniaMarcel Ausloos, Claudiu Herteliu & Adrian PanaGastro-enterology Department, University of Liège, Liège, BelgiumFloriane AusloosDepartment of Health Policy Planning and Management, University of Health and Allied Sciences, Ho, GhanaMartin Amogre AyanoreDepartment of Nursing, Debre Berhan University, Debre Berhan, EthiopiaYared Asmare AynalemDepartment of Reproductive Health, University of Gondar, Gondar, EthiopiaZelalem Nigussie AzenePublic Health Risk Sciences Division, Public Health Agency of Canada, Toronto, ON, CanadaAlaa BadawiDepartment of Nutritional Sciences, University of Toronto, Toronto, ON, CanadaAlaa BadawiUnit of Biochemistry, Sultan Zainal Abidin University (Universiti Sultan Zainal Abidin), Kuala Terengganu, MalaysiaAtif Amin BaigDepartment of Hypertension, Medical University of Lodz, Lodz, PolandMaciej BanachPolish Mothers’ Memorial Hospital Research Institute, Lodz, PolandMaciej BanachDepartment of Community Medicine, Gandhi Medical College Bhopal, Bhopal, IndiaNeeraj BediJazan University, Jazan, Saudi ArabiaNeeraj BediDepartment of Social and Clinical Pharmacy, Charles University, Hradec Kralova, Czech RepublicAkshaya Srikanth BhagavathulaInstitute of Public Health, United Arab Emirates University, Al Ain, United Arab EmiratesAkshaya Srikanth BhagavathulaSchool of Public Health, University of Adelaide, Adelaide, SA, AustraliaDinesh BhandariPublic Health Research Laboratory, Tribhuvan University, Kathmandu, NepalDinesh BhandariDepartment of Anatomy, Government Medical College Pali, Pali, IndiaNikha BhardwajDepartment of Community Medicine and Family Medicine, All India Institute of Medical Sciences, Jodhpur, IndiaPankaj BhardwajSchool of Public Health, All India Institute of Medical Sciences, Jodhpur, IndiaPankaj BhardwajDepartment of Statistical and Computational Genomics, National Institute of Biomedical Genomics, Kalyani, IndiaKrittika BhattacharyyaDepartment of Statistics, University of Calcutta, Kolkata, IndiaKrittika BhattacharyyaCentre for Global Child Health, University of Toronto, Toronto, ON, CanadaZulfiqar A. BhuttaCentre of Excellence in Women & Child Health, Aga Khan University, Karachi, PakistanZulfiqar A. BhuttaSocial Determinants of Health Research Center, Babol University of Medical Sciences, Babol, IranAli BijaniPlanning, Monitoring and Evaluation Directorate, Amhara Public Health Institute, Bahir Dar, EthiopiaTesega Tesega Mengistu BirhanuNutrition Department, St. Paul’s Hospital Millennium Medical College, Addis Ababa, EthiopiaZebenay Workneh BitewSt. Paul’s Hospital Millennium Medical College, Addis Ababa, EthiopiaZebenay Workneh BitewDepartment of Internal Medicine, Manipal Academy of Higher Education, Mangalore, IndiaArchith BoloorDepartment of Infectious Disease Epidemiology, London School of Hygiene & Tropical Medicine, London, UKOliver J. BradySchool of Public Health and Health Systems, University of Waterloo, Waterloo, ON, CanadaZahid A. ButtAl Shifa School of Public Health, Al Shifa Trust Eye Hospital, Rawalpindi, PakistanZahid A. ButtCentre for Population Health Sciences, Nanyang Technological University, Singapore, SingaporeJosip CarDepartment of Primary Care and Public Health, Imperial College London, London, UKJosip Car & Salman RawafResearch Unit on Applied Molecular Biosciences (UCIBIO), University of Porto, Porto, PortugalFelix CarvalhoDepartment of Medicine, University of Toronto, Toronto, ON, CanadaVijay Kumar ChattuGlobal Institute of Public Health (GIPH), Thiruvananthapuram, IndiaVijay Kumar ChattuMaternal and Child Health Division, International Centre for Diarrhoeal Disease Research, Bangladesh, Dhaka, BangladeshMohiuddin Ahsanul Kabir ChowdhuryDepartment of Epidemiology and Biostatistics, University of South Carolina, Columbia, SC, USAMohiuddin Ahsanul Kabir ChowdhuryFaculty of Biology, Hanoi National University of Education, Hanoi, VietnamDinh-Toi ChuLaboratory of Malaria Immunology and Vaccinology, National Institutes of Health, Bethesda, MD, USACamila H. CoelhoClinical Dermatology, IRCCS Istituto Ortopedico Galeazzi, University of Milan, Milan, ItalyGiovanni DamianiDepartment of Dermatology, Case Western Reserve University, Cleveland, OH, USAGiovanni DamianiDepartment of Public Health, Ambo University, Ambo, EthiopiaJiregna Darega GelaDepartment of Pediatrics, Tanta University, Tanta, EgyptAmira Hamed DarwishToxoplasmosis Research Center, Mazandaran University of Medical Sciences, Sari, IranAhmad DaryaniDivision of Women and Child Health, Aga Khan University, Karachi, PakistanJai K. DasWellcome Trust Brighton and Sussex Centre for Global Health Research, Brighton and Sussex Medical School, Brighton, UKKebede DeribeSchool of Public Health, Addis Ababa University, Addis Ababa, EthiopiaKebede DeribeSchool of Nursing and Midwifery, Haramaya University, Harar, EthiopiaAssefa DesalewDepartment of Community Medicine, University of Peradeniya, Peradeniya, Sri LankaSamath Dhamminda DharmaratneDepartment of Epidemiology and Biostatistics, Shahroud University of Medical Sciences, Shahroud, IranMostafa DianatinasabDepartment of Epidemiology, Shiraz University of Medical Sciences, Shiraz, IranMostafa DianatinasabCenter of Complexity Sciences, National Autonomous University of Mexico, Mexico City, MexicoDaniel DiazFaculty of Veterinary Medicine and Zootechnics, Autonomous University of Sinaloa, Culiacán Rosales, MexicoDaniel DiazDevelopment of Research and Technology Center, Ministry of Health and Medical Education, Tehran, IranShirin DjalaliniaDepartment of Medical Laboratory Sciences, Iran University of Medical Sciences, Tehran, IranFariba DorostkarInstitute of Microbiology and Immunology, University of Belgrade, Belgrade, SerbiaEleonora DubljaninSchool of Public Health, Hawassa University, Hawassa, EthiopiaBereket DukoSchool of Public Health, Curtin University, Perth, WA, AustraliaBereket Duko & Ted R. MillerCentre Clinical Epidemiology and Biostatistics, University of Newcastle, Newcastle, NSW, AustraliaAndem EffiongReference Laboratory of Egyptian Universities Hospitals, Ministry of Higher Education and Research, Cairo, EgyptMaysaa El Sayed ZakiPediatric Dentistry and Dental Public Health Department, Alexandria University, Alexandria, EgyptMaha El TantawiDepartment of Microbiology and Immunology, Suez Canal University, Ismailia, EgyptShymaa EnanyResearch Center for Environmental Determinants of Health, Kermanshah University of Medical Sciences, Kermanshah, IranNazir Fattahi & Masoud MoradiNational Institute for Stroke and Applied Neurosciences, Auckland University of Technology, Auckland, New ZealandValery L. FeiginResearch Center of Neurology, Moscow, RussiaValery L. FeiginAssociated Laboratory for Green Chemistry (LAQV), University of Porto, Porto, PortugalEduarda FernandesResearch Center on Public Health, University of Milan Bicocca, Monza, ItalyPietro FerraraInstitute of Gerontological Health Services and Nursing Research, Ravensburg-Weingarten University of Applied Sciences, Weingarten, GermanyFlorian FischerInstitute of Gerontology, National Academy of Medical Sciences of Ukraine, Kyiv, UkraineNataliya A. FoigtDepartment of Child Dental Health, Obafemi Awolowo University, Ile-Ife, NigeriaMorenike Oluwatoyin FolayanDepartment of Medical Parasitology, Abadan Faculty of Medical Sciences, Abadan, IranMasoud ForoutanDepartment of Dermatology, Kobe University, Kobe, JapanTakeshi FukumotoDepartment of Community Medicine, Datta Meghe Institute of Medical Sciences, Wardha, IndiaAbhay Motiramji Gaidhane, Zahiruddin Quazi Syed & Deepak SaxenaDepartment of Pediatric Nursing, Aksum University, Aksum, EthiopiaHailemikael Gebrekidan G. K. GebrekrstosSchool of Pharmacy, Aksum University, Aksum, EthiopiaLeake GebremeskelDepartment of Pharmacy, Mekelle University, Mekelle, EthiopiaLeake GebremeskelDepartment of Reproductive Health, Mekelle University, Mekelle, EthiopiaAssefa Ayalew GebreslassieTelethon Kids Institute, Perth Children’s Hospital, Nedlands, WA, AustraliaPeter W. GethingCurtin University, Bentley, WA, AustraliaPeter W. GethingDepartment of Biostatistics, Mekelle University, Mekelle, EthiopiaKebede Embaye GezaeInfectious Disease Research Center, Kermanshah University of Medical Sciences, Kermanshah, IranKeyghobad GhadiriPediatric Department, Kermanshah University of Medical Sciences, Kermanshah, IranKeyghobad GhadiriStudent Research Committee, Iran University of Medical Sciences, Tehran, IranAhmad GhashghaeeHealth Systems and Policy Research, Indian Institute of Public Health Gandhinagar, Gandhinagar, IndiaMahaveer GolechhaDepartment of Family and Community Medicine, University Of Sulaimani, Sulaimani, IraqMohammed Ibrahim Mohialdeen GubariDepartment of Pediatrics and Child Health, Mekelle University, Mekelle, EthiopiaFikaden Berhe HadguSchool of Health and Environmental Studies, Hamdan Bin Mohammed Smart University, Dubai, United Arab EmiratesSamer HamidiDepartment of Public Health, Wachemo University, Hossana, EthiopiaDemelash Woldeyohannes HandisoDepartment of Public Health, Jigjiga University, Jijiga, EthiopiaAbdiwahab Hashi & Muktar Omer OmerCenter for International Health (CIH), University of Bergen, Bergen, NorwayShoaib HassanBergen Center for Ethics and Priority Setting (BCEPS), University of Bergen, Bergen, NorwayShoaib HassanInstitute of Pharmaceutical Sciences, University of Veterinary and Animal Sciences, Lahore, PakistanKhezar HayatDepartment of Pharmacy Administration and Clinical Pharmacy, Xian Jiaotong University, Xian, ChinaKhezar HayatSchool of Business, London South Bank University, London, UKClaudiu HerteliuDepartment of Urban Planning and Design, University of Hong Kong, Hong Kong, ChinaHung Chak HoKasturba Medical College, Mangalore, Manipal Academy of Higher Education, Manipal, IndiaRamesh Holla & Priya RathiInstitute of Research and Development, Duy Tan University, Da Nang, VietnamMehdi Hosseinzadeh & Yasser VasseghianDepartment of Computer Science, University of Human Development, Sulaymaniyah, IraqMehdi HosseinzadehCollege of Science and Engineering, Hamad Bin Khalifa University, Doha, QatarMowafa HousehSchool of Pharmaceutical Sciences, University of Science Malaysia, Penang, MalaysiaRabia HussainDepartment of Occupational Safety and Health, China Medical University, Taichung, TaiwanBing-Fang HwangDepartment of Health Promotion and Education, University of Ibadan, Ibadan, NigeriaSegun Emmanuel IbitoyeDepartment of Community Medicine, University of Ibadan, Ibadan, NigeriaOlayinka Stephen IlesanmiDepartment of Community Medicine, University College Hospital, Ibadan, Ibadan, NigeriaOlayinka Stephen IlesanmiFaculty of Medicine, University of Belgrade, Belgrade, SerbiaIrena M. IlicDepartment of Epidemiology, University of Kragujevac, Kragujevac, SerbiaMilena D. IlicResearch Institute for Endocrine Sciences, Shahid Beheshti University of Medical Sciences, Tehran, IranSeyed Sina Naghibi IrvaniDepartment of Environmental Health Engineering, Guilan University of Medical Sciences, Rasht, IranJalil JaafariHealth Informatic Lab, Boston University, Boston, MA, USATahereh JavaheriDepartment of Community Medicine, Dr. Baba Saheb Ambedkar Medical College & Hospital, Delhi, IndiaRavi Prakash JhaDepartment of Community Medicine, Banaras Hindu University, Varanasi, IndiaRavi Prakash JhaDepartment of Ophthalmology, Heidelberg University, Heidelberg, GermanyJost B. JonasBeijing Institute of Ophthalmology, Beijing Tongren Hospital, Beijing, ChinaJost B. JonasDepartment of Family Medicine and Public Health, University of Opole, Opole, PolandJacek Jerzy JozwiakMinimally Invasive Surgery Research Center, Iran University of Medical Sciences, Tehran, IranAli KabirInstitute for Prevention of Non-communicable Diseases, Qazvin University of Medical Sciences, Qazvin, IranRohollah KalhorHealth Services Management Department, Qazvin University of Medical Sciences, Qazvin, IranRohollah KalhorDepartment of Forensic Medicine and Toxicology, All India Institute of Medical Sciences, Jodhpur, IndiaTanuj KanchanInstitute for Epidemiology and Social Medicine, University of Münster, Münster, GermanyAndré KarchInternational Research Center of Excellence, Institute of Human Virology Nigeria, Abuja, NigeriaGbenga A. KayodeJulius Centre for Health Sciences and Primary Care, Utrecht University, Utrecht, NetherlandsGbenga A. KayodeOpen, Distance and eLearning Campus, University of Nairobi, Nairobi, KenyaPeter Njenga KeiyoroDepartment of Public Health, Jordan University of Science and Technology, Irbid, JordanYousef Saleh KhaderDepartment of Global Health, University of Washington, Seattle, WA, USAIbrahim A. Khalil & Sonali KochharDepartment of Population Science, Jatiya Kabi Kazi Nazrul Islam University, Mymensingh, BangladeshMd Nuruzzaman KhanEpidemiology Department, Jazan University, Jazan, Saudi ArabiaMaseer KhanDepartment of Medical Microbiology & Immunology, United Arab Emirates University, Al Ain, United Arab EmiratesGulfaraz KhanFaculty of Health and Wellbeing, Sheffield Hallam University, Sheffield, UKKhaled KhatabCollege of Arts and Sciences, Ohio University, Zanesville, OH, USAKhaled KhatabDepartment of Medical Parasitology, Cairo University, Cairo, EgyptMona M. KhaterGlobal Evidence Synthesis Initiative, Datta Meghe Institute of Medical Sciences, Wardha, IndiaMahalaqua Nazli KhatibDepartment of Public Health, Kermanshah University of Medical Sciences, Kermanshah, IranNeda KianipourSchool of Traditional Chinese Medicine, Xiamen University Malaysia, Sepang, MalaysiaYun Jin KimDepartment of Nutrition, Simmons University, Boston, MA, USARuth W. KimokotiDepartment of Nursing and Health Promotion, Oslo Metropolitan University, Oslo, NorwaySezer KisaSchool of Health Sciences, Kristiania University College, Oslo, NorwayAdnan KisaGlobal Community Health and Behavioral Sciences, Tulane University, New Orleans, LA, USAAdnan KisaDepartment of Pediatrics, University of British Columbia, Vancouver, BC, CanadaNiranjan KissoonGlobal Healthcare Consulting, New Delhi, IndiaSonali KochharDepartment of Environmental Health Engineering, Arak University of Medical Sciences, Arak, IranAli KoolivandSchool of Population and Public Health, University of British Columbia, Vancouver, BC, CanadaJacek A. KopecArthritis Research Canada, Richmond, BC, CanadaJacek A. KopecCIBERSAM, San Juan de Dios Sanitary Park, Sant Boi de Llobregat, SpainAi KoyanagiCatalan Institution for Research and Advanced Studies (ICREA), Barcelona, SpainAi KoyanagiDepartment of Anthropology, Panjab University, Chandigarh, IndiaKewal KrishanInternational Institute for Population Sciences, Mumbai, IndiaPushpendra KumarFaculty of Health and Life Sciences, Coventry University, Coventry, UKOm P. KurmiDepartment of Medicine, McMaster University, Hamilton, ON, CanadaOm P. KurmiImperial College Business School, Imperial College London, London, UKDian KusumaFaculty of Public Health, University of Indonesia, Depok, IndonesiaDian KusumaPublic Health Foundation of India, Gurugram, IndiaDharmesh Kumar LalDepartment of Community and Family Medicine, University of Baghdad, Baghdad, IraqFaris Hasan LamiUnit of Genetics and Public Health, Institute of Medical Sciences, Las Tablas, PanamaIván LandiresMinistry of Health, Herrera, PanamaIván LandiresMedical Director, HelpMeSee, New York, NY, USAVan Charles LansinghGeneral Director, Mexican Institute of Ophthalmology, Queretaro, MexicoVan Charles LansinghDepartment of Otorhinolaryngology, Father Muller Medical College, Mangalore, IndiaSavita LasradoDepartment of Clinical Sciences and Community Health, University of Milan, Milan, ItalyCarlo La VecchiaSchool of Nursing, Hong Kong Polytechnic University, Hong Kong, ChinaPaul H. LeeCentre for Tropical Medicine and Global Health, University of Oxford, Oxford, UKSonia LewyckaOxford University Clinical Research Unit, Wellcome Trust Asia Programme, Hanoi, VietnamSonia LewyckaDepartment of Sociology, Shenzhen University, Shenzhen, ChinaBingyu LiDepartment of Systems, Populations, and Leadership, University of Michigan, Ann Arbor, MI, USAXuefeng LiuDepartment of Vector Biology, Liverpool School of Tropical Medicine, Liverpool, UKJoshua LongbottomIndependent Consultant, Melbourne, VIC, AustraliaAlan D. LopezRadiology Department, Egypt Ministry of Health and Population, Mansoura, EgyptHassan Magdy Abd El RazekGrants, Innovation and Product Development Unit, South African Medical Research Council, Cape Town, South AfricaPhetole Walter MahashaEnvironmental Health, Tehran University of Medical Sciences, Tehran, IranAfshin MalekiEnvironmental Health Research Center, Kurdistan University of Medical Sciences, Sanandaj, IranAfshin Maleki & Shadieh MohammadiInstitute for Social Science Research, The University of Queensland, Indooroopilly, QLD, AustraliaAbdullah A. MamunDepartment of Epidemiology and Biostatistics, Tehran University of Medical Sciences, Tehran, IranMohammad Ali MansourniaCampus Caucaia, Federal Institute of Education, Science and Technology of Ceará, Caucaia, BrazilFrancisco Rogerlândio Martins-MeloICF International, DHS Program, Rockville, MD, USABenjamin K. MayalaDepartment of Pharmacy, Wollo University, Dessie, EthiopiaBirhanu Geta MeharieDepartment of Medical Laboratory Sciences, Bahir Dar University, Bahir Dar, EthiopiaAddisu MelesePeru Country Office, United Nations Population Fund (UNFPA), Lima, PeruWalter MendozaForensic Medicine Division, Imam Abdulrahman Bin Faisal University, Dammam, Saudi ArabiaRitesh G. MenezesDepartment of Reproductive Health and Population Studies, Bahir Dar University, Bahir Dar, EthiopiaEndalkachew Worku MengeshaCenter for Translation Research and Implementation Science, National Institutes of Health, Bethesda, MD, USAGeorge A. MensahDepartment of Medicine, University of Cape Town, Cape Town, South AfricaGeorge A. MensahBreast Surgery Unit, Helsinki University Hospital, Helsinki, FinlandTuomo J. MeretojaUniversity of Helsinki, Helsinki, FinlandTuomo J. MeretojaClinical Microbiology and Parasitology Unit, Dr. Zora Profozic Polyclinic, Zagreb, CroatiaTomislav MestrovicUniversity Centre Varazdin, University North, Varazdin, CroatiaTomislav MestrovicPacific Institute for Research & Evaluation, Calverton, MD, USATed R. MillerInternal Medicine Programme, Kyrgyz State Medical Academy, Bishkek, KyrgyzstanErkin M. MirrakhimovDepartment of Atherosclerosis and Coronary Heart Disease, National Center of Cardiology and Internal Disease, Bishkek, KyrgyzstanErkin M. MirrakhimovHeidelberg Institute of Global Health (HIGH), Heidelberg University, Heidelberg, GermanyBabak Moazen & Shafiu MohammedInstitute of Addiction Research (ISFF), Frankfurt University of Applied Sciences, Frankfurt, GermanyBabak MoazenDepartment of Biostatistics, Hamadan University of Medical Sciences, Hamadan, IranNaser Mohammad Gholi MezerjiResearch Institute for Health Development, Kurdistan University of Medical Sciences, Sanandaj City, IranShadieh MohammadiHealth Systems and Policy Research Unit, Ahmadu Bello University, Zaria, NigeriaShafiu MohammedComputer, Electrical, and Mathematical Sciences and Engineering Division, King Abdullah University of Science and Technology, Thuwal, Saudi ArabiaPaula MoragaClinical Research Development Center, Kermanshah University of Medical Sciences, Kermanshah, IranMehdi NaderiResearch and Analytics Department, Initiative for Financing Health and Human Development, Chennai, IndiaAhamarshan Jayaraman NagarajanDepartment of Research and Analytics, Bioinsilico Technologies, Chennai, IndiaAhamarshan Jayaraman NagarajanDepartment of Pediatrics, Arak University of Medical Sciences, Arak, IranJavad NazariDisease Control and Environmental Health, Makerere University, Kampala, UgandaRawlance NdejjoDepartment of General Surgery, Carol Davila University of Medicine and Pharmacy, Bucharest, RomaniaIonut NegoiDepartment of General Surgery, Emergency Hospital of Bucharest, Bucharest, RomaniaIonut NegoiDepartment of Biological Sciences, University of Embu, Embu, KenyaJosephine W. NgunjiriInstitute for Global Health Innovations, Duy Tan University, Hanoi, VietnamHuong Lan Thi Nguyen & Hai Quang PhamSouth African Medical Research Council, Cape Town, South AfricaChukwudi A. Nnaji & Charles Shey WiysongeSchool of Public Health and Family Medicine, University of Cape Town, Cape Town, South AfricaChukwudi A. Nnaji & Charles Shey WiysongeCentre for Heart Rhythm Disorders, University of Adelaide, Adelaide, SA, AustraliaJean Jacques NoubiapUnit of Microbiology and Public Health, Institute of Medical Sciences, Las Tablas, PanamaVirginia Nuñez-SamudioDepartment of Public Health, Ministry of Health, Herrera, PanamaVirginia Nuñez-SamudioDepartment of Psychiatry and Behavioural Neurosciences, McMaster University, Hamilton, ON, CanadaAndrew T. OlagunjuDepartment of Psychiatry, University of Lagos, Lagos, NigeriaAndrew T. OlagunjuCentre for Healthy Start Initiative, Lagos, NigeriaJacob Olusegun Olusanya & Bolajoko Olubukunola OlusanyaDepartment of Pharmacology and Therapeutics, University of Nigeria Nsukka, Enugu, NigeriaObinna E. OnwujekweLaboratory of Public Health Indicators Analysis and Health Digitalization, Moscow Institute of Physics and Technology, Dolgoprudny, RussiaNikita Otstavnov & Stanislav S. OtstavnovDepartment of Project Management, National Research University Higher School of Economics, Moscow, RussiaStanislav S. OtstavnovDepartment of Medicine, University of Ibadan, Ibadan, NigeriaMayowa O. OwolabiDepartment of Medicine, University College Hospital, Ibadan, Ibadan, NigeriaMayowa O. OwolabiDepartment of Respiratory Medicine, Jagadguru Sri Shivarathreeswara Academy of Health Education and Research, Mysore, IndiaMahesh P ADepartment of Forensic Medicine, Manipal Academy of Higher Education, Mangalore, IndiaJagadish Rao PadubidriDepartment of Health Metrics, Center for Health Outcomes & Evaluation, Bucharest, RomaniaAdrian PanaSchool of Global Public Health, New York University, New York, NY, USAEmmanuel K. PeprahDepartment of Parasitology and Entomology, Tarbiat Modares University, Tehran, IranMajid PirestaniUniversity Medical Center Groningen, University of Groningen, Groningen, NetherlandsMaarten J. PostmaSchool of Economics and Business, University of Groningen, Groningen, NetherlandsMaarten J. PostmaDepartment of Pharmacology, Imam Abdulrahman Bin Faisal University, Dammam, Saudi ArabiaFaheem Hyder PottooDepartment of Nutrition and Food Sciences, Maragheh University of Medical Sciences, Maragheh, IranHadi PourjafarDietary Supplements and Probiotic Research Center, Alborz University of Medical Sciences, Karaj, IranHadi PourjafarThalassemia and Hemoglobinopathy Research Center, Ahvaz Jundishapur University of Medical Sciences, Ahvaz, IranFakher RahimMetabolomics and Genomics Research Center, Tehran University of Medical Sciences, Tehran, IranFakher RahimSina Trauma and Surgery Research Center, Tehran University of Medical Sciences, Tehran, IranVafa Rahimi-MovagharDepartment of Community Medicine, Maharishi Markandeshwar Medical College & Hospital, Solan, IndiaMohammad Hifz Ur RahmanDepartment of Oral Pathology, Srinivas Institute of Dental Sciences, Mangalore, IndiaSowmya J. RaoAcademic Public Health England, Public Health England, London, UKSalman RawafWHO Collaborating Centre for Public Health Education and Training, Imperial College London, London, UKDavid Laith RawafUniversity College London Hospitals, London, UKDavid Laith RawafSchool of Health, Medical and Applied Sciences, CQ University, Sydney, NSW, AustraliaLal RawalDepartment of Computer Science, Boston University, Boston, MA, USAReza RawassizadehSchool of Public Health, Haramaya University, Harar, EthiopiaLemma Demissie RegassaSchool of Social Sciences and Psychology, Western Sydney University, Penrith, NSW, AustraliaAndre M. N. RenzahoTranslational Health Research Institute, Western Sydney University, Penrith, NSW, AustraliaAndre M. N. RenzahoNetwork of Immunity in Infection, Malignancy and Autoimmunity (NIIMA), Universal Scientific Education and Research Network (USERN), Tehran, IranNima RezaeiPediatric Infectious Diseases Research Center, Mazandaran University of Medical Sciences, Sari, IranMohammad Sadegh RezaiEpidemiology Research Unit Institute of Public Health (EPIUnit-ISPUP), University of Porto, Porto, PortugalAna Isabel RibeiroDepartment of Surgery, University of Minnesota, Minneapolis, MN, USAJennifer RickardDepartment of Surgery, University Teaching Hospital of Kigali, Kigali, RwandaJennifer RickardFaculty of Medical Sciences, Research Department, National University of Caaguazu, Cnel. Oviedo, ParaguayCarlos Miguel Rios-GonzálezDepartment of Research and Publications, National Institute of Health, Asunción, ParaguayCarlos Miguel Rios-GonzálezDepartment of Health Statistics, National Institute for Medical Research, Dar es Salaam, TanzaniaSusan Fred RumishaDepartment of Epidemiology, Shahid Beheshti University of Medical Sciences, Tehran, IranSiamak SabourDepartment of Phytochemistry, Soran University, Soran, IraqS. Mohammad SajadiDepartment of Nutrition, Cihan University-Erbil, Kurdistan Region, IraqS. Mohammad SajadiCenter for Health Policy & Center for Primary Care and Outcomes Research, Stanford University, Stanford, CA, USAJoshua A. SalomonDrug Applied Research Center, Tabriz University of Medical Sciences, Tabriz, IranHossein Samadi KafilDepartment of Entomology, Ain Shams University, Cairo, EgyptAbdallah M. SamyDepartment of Surgery, Marshall University, Huntington, WV, USAJuan SanabriaDepartment of Nutrition and Preventive Medicine, Case Western Reserve University, Cleveland, OH, USAJuan SanabriaFaculty of Infectious and Tropical Diseases, London School of Hygiene & Tropical Medicine, London, UKBenn SartoriusDepartment of Epidemiology, Indian Institute of Public Health, Gandhinagar, IndiaDeepak SaxenaGlobal Programs, Medical Teams International, Seattle, WA, USALauren E. SchaefferDepartment of Pediatric Newborn Medicine, Brigham and Women’s Hospital, Boston, MA, USALauren E. SchaefferEmergency Department, Manian Medical Centre, Erode, IndiaSubramanian SenthilkumaranCenter for Biomedical Information Technology, Shenzhen Institutes of Advanced Technology, Shenzhen, ChinaFeng ShaPublic Health Division, An-Najah National University, Nablus, PalestineAmira A. ShaheenIndependent Consultant, Karachi, PakistanMasood Ali ShaikhUniversity School of Management and Entrepreneurship, Delhi Technological University, Delhi, IndiaRajesh SharmaCentre for Medical Informatics, University of Edinburgh, Edinburgh, UKAziz SheikhDivision of General Internal Medicine, Harvard University, Boston, MA, USAAziz SheikhInstitute for Population Health, King’s College London, London, UKKenji ShibuyaNational Institute of Infectious Diseases, Tokyo, JapanMika ShigematsuCollege of Medicine, Yonsei University, Seoul, South KoreaJae Il ShinDepartment of Law, Economics, Management and Quantitative Methods, University of Sannio, Benevento, ItalyBiagio SimonettiWSB University in Gdańsk, Gdansk, PolandBiagio SimonettiSchool of Medicine, University of Alabama at Birmingham, Birmingham, AL, USAJasvinder A. SinghMedicine Service, US Department of Veterans Affairs (VA), Birmingham, AL, USAJasvinder A. SinghNursing Care Research Center, Semnan University of Medical Sciences, Semnan, IranAmin SoheiliDepartment of Infectious Diseases, Kharkiv National Medical University, Kharkiv, UkraineAnton SokhanDivision of Community Medicine, International Medical University, Kuala Lumpur, MalaysiaChandrashekhar T. SreeramareddyDepartment of Community Medicine, Ahmadu Bello University, Zaria, NigeriaMu’awiyyah Babale SufiyanSchool of Medicine, University of California San Francisco, San Francisco, CA, USAScott J. SwartzJoint Medical Program, University of California Berkeley, Berkeley, CA, USAScott J. SwartzDepartment of Nursing, Aksum University, Aksum, EthiopiaDegena Bahrey TadesseDepartment of Midwifery, University of Gondar, Gondar, EthiopiaAnimut Tagele TamiruDepartment of Clinical Pharmacy, University of Gondar, Gondar, EthiopiaYonas Getaye TeferaDepartment of Epidemiology and Biostatistics, University of Gondar, Gondar, EthiopiaZemenu Tadesse TessemaK.A. Timiryazev Institute of Plant Physiology, Russian Academy of Sciences, Moscow, RussiaMariya Vladimirovna TitovaLaboratory of Public Health Indicators Analysis and Health Digitalization, Moscow Institute of Physics and Technology, Moscow, RussiaMariya Vladimirovna TitovaDepartment of Health Economics, Hanoi Medical University, Hanoi, VietnamBach Xuan TranFaculty of Geo-Information Science and Earth Observation, University of Twente, Enschede, NetherlandsPhuong N. TruongKasturba Medical College, Manipal Academy of Higher Education, Mangalore, IndiaBhaskaran UnnikrishnanAmity Institute of Biotechnology, Amity University Rajasthan, Jaipur, IndiaEra UpadhyayUKK Institute, Tampere, FinlandTommi Juhani VasankariDepartment of Medical and Surgical Sciences, University of Bologna, Bologna, ItalyFrancesco S. ViolanteOccupational Health Unit, Sant’Orsola Malpighi Hospital, Bologna, ItalyFrancesco S. ViolanteCenter of Excellence in Behavioral Medicine, Nguyen Tat Thanh University, Ho Chi Minh City, VietnamGiang Thu VuFoundation University Medical College, Foundation University Islamabad, Islamabad, PakistanYasir WaheedCultures, Societies and Global Studies, & Integrated Initiative for Global Health, Northeastern University, Boston, MA, USARichard G. WamaiSchool of Public Health, University of Nairobi, Nairobi, KenyaRichard G. WamaiDepartment of Human Nutrition and Food Sciences, Debre Markos University, Debre Markos, EthiopiaEmebet Gashaw WassieDepartment of Midwifery, Adigrat University, Adigrat, EthiopiaFissaha Tekulu WelayDepartment of Community Medicine, Rajarata University of Sri Lanka, Anuradhapura, Sri LankaNuwan Darshana WickramasingheDepartment of Epidemiology, Johns Hopkins University, Baltimore, MD, USAKirsten E. WiensDepartment of Neurology, University of Melbourne, Melbourne, VIC, AustraliaTissa WijeratneDepartment of Medicine, University of Rajarata, Saliyapura Anuradhapuraya, Sri LankaTissa WijeratneDepartment of Public Health, Samara University, Samara, EthiopiaTemesgen Gebeyehu WondmenehDepartment of Diabetes and Metabolic Diseases, University of Tokyo, Tokyo, JapanTomohide YamadaSchool of International Development and Global Studies, University of Ottawa, Ottawa, ON, CanadaSanni YayaThe George Institute for Global Health, University of Oxford, Oxford, UKSanni YayaDepartment of Nursing, Arba Minch University, Arba Minch, EthiopiaYordanos Gizachew YeshitilaCentre for Suicide Research and Prevention, University of Hong Kong, Hong Kong, ChinaPaul YipDepartment of Social Work and Social Administration, University of Hong Kong, Hong Kong, ChinaPaul YipDepartment of Neuropsychopharmacology, National Center of Neurology and Psychiatry, Kodaira, JapanNaohiro YonemotoDepartment of Public Health, Juntendo University, Tokyo, JapanNaohiro YonemotoDepartment of Epidemiology and Biostatistics, Wuhan University, Wuhan, ChinaChuanhua YuCancer Institute, Hacettepe University, Ankara, TurkeyDeniz YuceDepartment of Health Care Management and Economics, Urmia University of Medical Science, Urmia, IranHasan YusefzadehDepartment of Medicine, University Ferhat Abbas of Setif, Sétif, AlgeriaZoubida ZaidiSocial Development and Health Promotion Research Center, Kermanshah University of Medical Sciences, Kermanshah, IranAlireza ZangenehSchool of Medicine, Wuhan University, Wuhan, ChinaZhi-Jiang ZhangSchool of Public Health, Wuhan University of Science and Technology, Wuhan, ChinaYunquan ZhangHubei Province Key Laboratory of Occupational Hazard Identification and Control, Wuhan University of Science and Technology, Wuhan, ChinaYunquan ZhangDepartment of Health Education and Health Promotion, Kermanshah University of Medical Sciences, Kermanshah, IranArash ZiapourManaging the estimation or publication process. L.B.D. T.B.C. Writing the first draft of the manuscript. R.C.R.J. Primary responsibility for this manuscript focused on: applying analytical methods to produce estimates. L.B.D. T.B.C. Primary responsibility for this manuscript focused on: seeking, cataloguing, extracting, or cleaning data; production or coding of figures and tables. L.B.D. T.B.C. Providing data or critical feedback on data sources. L.B.D. T.B.C. and S.I.H. Development of methods or computational machinery. R.C.R.J. and L.B.D. T.B.C. Providing critical feedback on methods or results. L.B.D. T.B.C. and S.I.H. Drafting the manuscript or revising it critically for important intellectual content. R.C.R.J., L.B.D. T.B.C., and S.I.H. Management of the overall research enterprise (for example, through membership in the Scientific Council). L.B.D. T.B.C. and S.I.H. Consortia author contributions Managing the estimation or publication process. B.F.B., M.K.M.P. Writing the first draft of the manuscript. R.C.R.J. Primary responsibility for this manuscript focused on: applying analytical methods to produce estimates. C.A.W. Primary responsibility for this manuscript focused on: seeking, cataloguing, extracting, or cleaning data; production or coding of figures and tables. M.M.B. Providing data or critical feedback on data sources D.J.W., A.D., C.E.T., H.A., A.E.A., V.A., O.O.A., M.A., B.A., F.A., S.A., H.A., S.M.A., A.A.-H., N.A.-G., A.T.A., S.A., C.L.A., F.A., D.A., S.C.Y.A., J.A., O.A., M.A., F.A., Y.A.A., A.B., M.B., N.B., A.S.B., A.B., V.K.C., D.-T.C., G.D., J.D.G., A.D., S.D.D., M.D., A.E., M.El.S.Z., S.E., T.F., A.M.G., L.G., P.W.G., K.G., A.G., M.G., A.H., S.H., K.H., C.H., H.C.H., M.H., M.H., S.S.N.I., T.J., J.B.J., J.J.J., A.K., G.A.K., Y.S.K., I.A.K., M.N.K., M.K., K.K., M.N.K., Y.J.K., S.K., A.K., N.K., K.K., P.K., D.K., D.K.L., F.H.L., V.C.L., S.L., A.L.-A., K.E.L., S.S.L., P.A.L., X.L., H.M.A.E.R., M.A.M., B.K.M., W.M., R.G.M., E.M.M., B.M., N.M.G.M., S.M., S.M., A.H.M., M.M., A.J.N., J.N., I.N., J.W.N., Q.P.N., H.L.T.N., C.A.N., J.J.N., A.T.O., J.O.O., B.O.O., O.E.O., N.O., S.S.O., M.O.O., M.P.A., J.R.P., A.P., E.K.P., H.Q.P., M.P., M.J.P., H.P., Z.Q.S., F.R., V.R.-M., S.J.R., P.R., S.R., D.L.R., L.R., R.R., A.M.N.R., N.R., J.R., C.M.R.-G., S.S., S.M.S., A.M.S., B.S., D.S., A.A.S., M.A.S., J.I.S., J.A.S., A.S., E.S., C.T.S., S.J.S., D.B.T., A.T.T., B.X.T., P.N.T., B.U., E.U., T.J.V., Y.V., G.T.V., Y.W., R.G.W., T.W., C.S.W., T.G.W., S.Y., Y.G.Y., N.Y., C.Y., H.Y., Z.Z., A.Z., and S.I.H. Development of methods or computational machinery R.C.R.J., C.A.W., M.M.B., A.D., L.E., S.B., C.E.T., H.A., D.A., Y.A.A., A.S.B., D.C.C., V.K.C., F.D. A.D., M.D., M.E.S.Z., N.F., J.J.F., P.W.G., M.H., K.B.J., S.K., A., A.D.L., S.M., A.H.M., J.W.N., Q.P.N., S.F.R., A.M.S., E.E.S., S.J.S., E.U., Y.V., K.E.W., Y.G.Y., and N.Y. Providing critical feedback on methods or results C.A.W., A.D., C.E.T., H.A., A.E.A., E.A.-G., V.A., O.O.A., M.A., B.A., F.A., Z.A.-A., R.K.A., S.A., H.A., A.A.-H., H.M.A.M., K.A.A., N.A.-Gu., A.T.A., S.A., D.A.A., C.L.A., F.A., D.A., S.C.Y.A., J.A., O.A., M.M.W.A., M.A., F.A., Y.A.A., Z.N.A., A.B., M.B., A.S.B., D.B., N.B., P.B., K.B., O.J.B., Z.A.B., A.B., Z.W.B., A.B., Z.A.B., V.C., M.A.K.C., D.-T.C., C.H.C., G.D., J.D.G., A.H.D., A.D., J.K.D., K.D., A.D., S.D.D., M.D., D.D., S.D., F.D., B.D., L.D.-L., A.E., V.L.F., F.F., N.A.F., M.O.F., M.F., T.F., A.M.G., H.G.G.K.G., L.G., A.A.G., K.E.G., A.G., M.G., F.B.H., S.H., A.H., S.H., C.H., H.C.H., R.H., M.H., M.H., R.H., B.-F.H., S.E.I., O.S.I., I.M.I., M.D.I., S.S.N.I., T.J., R.P.J., J.B.J., J.J.J., A.K., R.K., T.K., A.K., G.A.K., P.N.K., Y.S.K., I.A.K., M.N.K., M.K., K.K., M.M.K., M.N.K., Y.J.K., R.W.K., S.K., A.K., N.K., S.K., A.K., J.A.K., A.K., K.K., P.K., O.P.K., D.K., D.K.L., S.L., K.E.L., S.L., B.L., X.L., A.D.L., H.M.A.E.R., P.W.M., A.A.M., M.A.M., L.B.M., F.R.M.-M., B.K.M., W.M., R.G.M., E.W.M., T.J.M., T.R.M., E.M.M., B.M., N.M.G.M., S.M., S.M., A.H.M., R.M., J.F.M., M.N., A.J.N., J.N., R.N., I.N., J.W.N., H.L.T.N., C.A.N., J.J.N., A.T.O., J.O.O., B.O.O., M.O.O., O.E.O., N.O., S.S.O., M.O.O., M.P.A., J.R.P., A.P., E.K.P., H.Q.P., M.J.P., F.H.P., H.P., Z.Q.S., F.R., V.R.-M., S.J.R., P.R., S.R., D.L.R., L.R., R.R., L.D.R., A.M.N.R., N.R., M.S.R., A.I.R., J.R., C.M.R.-G., S.S., S.M.S., J.A.S., H.S.K., A.M.S., J.S., B.S., D.S., L.E.S., S.S., F.S., A.A.S., M.A.S., A.S., K.S., M.S., J.I.S., B.S., J.A.S., D.L.S., A.S., E.E.S., C.T.S., M.B.S., D.B.T., A.T.T., Y.G.T., M.-H.T., Z.T.T., M.V.T., B.X.T., P.N.T., B.U., E.U., Y.V., F.S.V., G.T.V., Y.W., R.G.W., E.G.W., F.T.W., N.D.W., K.E.W., T.W., C.S.W., T.G.W., T.Y., S.Y., Y.G.Y., P.Y., N.Y., C.Y., D.Y., Z.Z., M.Z., Z.-J.Z., Y.Z., and S.I.H. Drafting the manuscript or revising it critically for important intellectual content R.C.R.J., C.A.W., M.K.M.-P., L.E., H.A., E.A.-G., V.A., O.O.A., M.A., B.A, F.A., R.K.A., H.A., A.A.-H., N.A.-G., A.T.A., S.A., D.A.A., R.A., C.L.A., J.A., O.A., M.M.W.A., M.A., F.A., M.A.A., Z.N.A., A.B., A.A.B., M.B., N.B. A.S.B., D.B., K.B., T.T.M.B., O.J.B., J.C., F.C., V.K.C., G.D., A.D., N.D.W., K.D., S.D.D., D.D., E.D., A.E., M.E.S.Z., M.E.T., S.E., V.L.F., E.F., P.F., F.F., N.A.F., M.O.F., M.F., T.F., A.M.G., L.G., A.G., M.I.M.G., D.W.H., A.H., S.H., C.H., H.C.H., R.H., M.H., S.E.I., O.S.I., I.M.I., M.D.I., S.S.N.I., J.J., R.P.J., J.B.J., J.J.J., A.K., A.K., G.A.K., M.N.K., M.K., G.K., K.K., M.M.K., M.N.K., A.K., N.K., A.K., A.K., K.K., P.K., O.P.K., D.K., I.L., S.L., C.L.V., P.H.L., K.E.L., J.L., A.D.L., H.M.A.E.R., P.W.M., A.M., A.A.M., M.A.M., L.B.M., F.R.M.-M., B.G.M., W.M., R.G.M., E.W.M., G.A.M., T.J.M., T.M., T.R.M., B.M., S.M., S.M., A.H.M., R.M., P.M., J.F.M., A.J.N., J.N., I.N., J.W.N., H.L.T.N., V.N.-S., A.T.O., J.O.O., B.O.O., M.O.O., O.E.O., N.O., S.S.O., M.O.O., M.P.A., J.R.P., A.P., H.Q.P., M.J.P., Z.Q.S., F.R., V.R.-M., M.H.U.R., S.J.R., S.R., D.L.R., L.R., N.R., A.I.R., J.R., C.M.R.-G., S.F.R., S.S., J.A.S., H.S.K., A.M.S., J.S., D.S., R.S., M.S., J.A.S., A.S., C.T.S., M.B.S., D.B.T., A.T.T., M.V.T., B.X.T., B.U., E.U., T.J.V., Y.V., F.S.V., G.T.V., R.G.W., N.D.W., K.E.W., T.W., .C.S.W., S.Y., Y.G.Y., Z.Z., M.Z., Z.-J.Z., and S.I.H. Management of the overall research enterprise (for example, through membership in the Scientific Council) B.F.B., A.J.C., P.W.G., J.A.K., A.H.M., C.J.L.M., P.C.R., J.A.S., B.S., and S.I.H. More

Zachos, J. C., Pagani, M., Sloan, L., Thomas, E. & Billups, K. Trends, rhythms, and aberrations in global climate 65 Ma to present. Science 292, 686–693 (2001).Article 
CAS 
PubMed 
ADS 

Google Scholar 
Erwin, D. H. Climate as a driver of evolutionary change. Curr. Biol. 19, R575–R583 (2009).Article 
CAS 
PubMed 

Google Scholar 
Mayhew, P. J., Jenkins, G. B. & Benton, T. G. A long-term association between global temperature and biodiversity, origination and extinction in the fossil record. Proc. R. Soc. Lond. B 275, 47–53 (2008).
Google Scholar 
Raia, P. et al. Past extinctions of Homo species coincided with increased vulnerability to climatic change. One Earth 3, 480–490 (2020).Article 
ADS 

Google Scholar 
deMencoal, P. Climate and human evolution. Science 331, 540–542 (2011).Article 
ADS 

Google Scholar 
Stroud, J. T. & Losos, J. B. Ecological opportunity and adaptive radiation. Annu. Rev. Ecol. Evol. Syst. 47, 507–532 (2016).Article 

Google Scholar 
Potts, R. & Faith, J. T. Alternating high and low climate variability: The context of natural selection and speciation in Plio-Pleistocene hominin evolution. J. Hum. Evol. 87, 5–20 (2015).Article 
PubMed 

Google Scholar 
Clavel, J. & Morlon, H. Accelerated body size evolution during cold climatic periods in the Cenozoic. Proc. Natl Acad. Sci. USA 114, 4183–4188 (2017).Article 
CAS 
PubMed 
PubMed Central 
ADS 

Google Scholar 
Mihlbachler, M. C., Rivals, F., Solounias, N. & Semprebon, G. M. Dietary change and evolution of horses in North America. Science 331, 1178–1181 (2011).Article 
CAS 
PubMed 
ADS 

Google Scholar 
Mennecart, B. et al. Bony labyrinth morphology clarifies the origin and evolution of deer. Sci. Rep. 7, 13176 (2017).Article 
PubMed 
PubMed Central 
ADS 

Google Scholar 
Ponce, deLeón et al. Human bony labyrinth is an indicator of population history and dispersal from Africa. Proc. Natl Acad. Sci. USA 115, 4128–4133 (2018).Article 
ADS 

Google Scholar 
Luo, Z.-X. The inner ear and its bony housing in tritylodontids and implications for the evolution of the mammalian ear. Bull. Mus. Comp. Zool. 156, 81–97 (2001).
Google Scholar 
Ekdale, E. G. Comparative anatomy of the bony labyrinth (inner ear) of placental mammals. PLoS ONE 8, e66624 (2013).Article 
CAS 
PubMed 
PubMed Central 
ADS 

Google Scholar 
O’Leary, M. A. An anatomical and phylogenetic study of the osteology of the petrosal of extant and extinct artiodactylans (Mammalia) and relatives. Bull. Am. Mus. Nat. Hist. 335, 1–206 (2010).Article 

Google Scholar 
Costeur, L. et al. The bony labyrinth of toothed whales reflects both phylogeny and habitat preferences. Sci. Rep. 8, 7841 (2018).Article 
PubMed 
PubMed Central 
ADS 

Google Scholar 
Spoor, F., Bajpai, S., Hussain, S. T., Kumar, K. & Thewissen, J. G. M. Vestibular evidence for the evolution of aquatic behavior in early cetaceans. Nature 417, 163–166 (2002).Article 
CAS 
PubMed 
ADS 

Google Scholar 
Davies, K. T. J., Bates, P. J. J., Maryanto, I., Cotton, J. A. & Rossiter, S. J. The evolution of bat vestibular systems in the face of potential antagonistic selection pressures for flight and echolocation. PLoS ONE 8, e61998 (2013).Article 
CAS 
PubMed 
PubMed Central 
ADS 

Google Scholar 
Park, T., Mennecart, B., Costeur, L., Grohé, C. & Cooper, N. Convergent evolution in toothed whale cochleae. BMC Evol. Biol. 1, 195 (2019).Article 

Google Scholar 
Benoit, J. et al. A test of the lateral semicircular canal correlation to head posture, diet and other biological traits in “ungulate” mammals. Sci. Rep. 10, 19602 (2020).Article 
CAS 
PubMed 
PubMed Central 
ADS 

Google Scholar 
Morimoto, N. et al. Variation of bony labyrinthine morphology in Mio-Plio-Pleistocene and modern anthropoids. Am. J. Phys. Anthropol. 2020, 1–17 (2020).
Google Scholar 
DeMiguel, D., Azanza, B. & Morales, J. Key innovations in ruminant evolution: A paleontological perspective. Int. Zool. 9, 412–433 (2014).Article 

Google Scholar 
Gunz, P., Ramsier, M., Kuhrig, M., Hublin, J. & Spoor, F. The mammalian bony labyrinth reconsidered, introducing a comprehensive geometric morphometric approach. J. Anat. 220, 529–543 (2012).Article 
PubMed 
PubMed Central 

Google Scholar 
Grohe, C., Tseng, Z. J., Lebrun, R., Boistel, R. & Flynn, J. J. Bony labyrinth shape variation in extant Carnivora: a case study of Musteloidea. J. Anat. 228, 366–383 (2015).Article 
PubMed 
PubMed Central 

Google Scholar 
Urciuoli, A. et al. A comparative analysis of the vestibular apparatus in Epipliopithecus vindobonensis: Phylogenetic implications. J. Hum. Evol. 151, 102930 (2021).Article 
PubMed 

Google Scholar 
IUCN. The IUCN red list of threatened species. Version 2021-1. https://www.iucnredlist.org. Accessed 17 June 2021.Kingdon, J. & Hoffmann. M. Mammals of Africa. Volume VI pigs, hippopotamuses, chevrotains, Giraffes, deer and bovids 704 (Bloomsbury Publishing, 2013).Chen, L. et al. Large-scale ruminant genome sequencing provides insights into their evolution and distinct traits. Science 364 eaav6202 (2019).Hassanin, A. et al. Pattern and timing of diversification of Cetartiodactyla (Mammalia, Laurasiatheria), as revealed by a comprehensive analysis of mitochondrial genomes. C. R. Biol. 335, 32–50 (2012).Article 
PubMed 

Google Scholar 
Wang, Y. et al. Genetic basis of ruminant headgear and rapid antler regeneration. Science 364, 1153 (2019).Article 

Google Scholar 
Myers, E. A. & Bubrink, F. T. Ecological opportunity: Trigger of adaptative radiation. Nat. Educ. Knowl. 3, 23 (2012).
Google Scholar 
Gentry, A. W. Bovidae. In Cenozoic mammals of Africa (eds Werdelin, L. & Sanders, W. J.) 741–796 (University of California Press, 2010).Harris, J. M., Solounias, N. & Geraads, D. Giraffoidea. In Werdelin, L. & Sanders, W. J. Cenozoic mammals of Africa. 797–812 (University of California Press, 2010).Clauss, M. & Rössner, G. E. Old world ruminant morphophysiology, life history, and fossil record: exploring key innovations of a diversification sequence. Ann. Zool. Fenn. 51, 80–94 (2014).Article 

Google Scholar 
Johnston, A. R. & Anthony, N. M. A multi-locus species phylogeny of African forest duikers in the subfamily Cephalophinae: evidence for a recent radiation in the Pleistocene. BMC Evol. Biol. 12, 120 (2012).Article 
PubMed 
PubMed Central 

Google Scholar 
Cooney, C. R. & Thomas, G. H. Heterogeneous relationships between rates of speciation and body size evolution across vertebrate clades. Nat. Ecol. Evol. 5, 101–110 (2020).Article 
PubMed 

Google Scholar 
Köhler, M. & Moyà-Solà, S. Physiological and life history strategies of a fossil large mammal in a resource-limited environment. Proc. Natl Acad. Sci. USA 106, 20354–22035 (2009).Article 
PubMed 
PubMed Central 
ADS 

Google Scholar 
Bibi, F. A multi-calibrated mitochondrial phylogeny of extant Bovidae (Artiodactyla, Ruminantia) and the importance of the fossil record to systematics. BMC Evol. Biol. 13, 1–15 (2013).Article 

Google Scholar 
Geraads, D. A reassessment of the Bovidae (Mammalia) from the Nawata Formation of Lothagam, Kenya, and the late Miocene diversification of the family in Africa. J. Syst. Palaeontol. 17, 1–14 (2017).
Google Scholar 
Mennecart, B., Aiglstorfer, M., Li, Y., Li, C. & Wang, S. Ruminants reveal Eocene Asiatic palaeobiogeographical provinces as the origin of diachronous mammalian Oligocene dispersals into Europe. Sci. Rep. 11, 17710 (2021).Article 
CAS 
PubMed 
PubMed Central 
ADS 

Google Scholar 
Rössner, G. E. Family tragulidae. In: The evolution of artiodactyls (eds Prothero, D. R. & Foss S. C.) (The Johns Hopkins University Press, Baltimore, 2007).Sánchez, I. M., Quiralte, V., Morales, J. & Pickford, M. A new genus of tragulid ruminant from the early Miocene of Kenya. Acta Palaeontol. Pol. 55, 177–187 (2010).Article 

Google Scholar 
Sánchez, I. M., Quiralte, V., Ríos, M., Morales, J. & Pickford, M. First African record of the Miocene Asian mouse-deer Siamotragulus (Mammalia, Ruminantia, Tragulidae): implications for the phylogeny and evolutionary history of the advanced selenodont tragulids. J. Syst. Palaeontol. 13, 543–556 (2015).Article 

Google Scholar 
Mennecart, B. et al. The first French tragulid skull (Mammalia, Ruminantia, Tragulidae) and associated tragulid remains from the Middle Miocene of Contres (Loir-et-Cher, France). C. R. Palevol 17, 189–200 (2018).Article 

Google Scholar 
Bobe, R. & Eck, G. C. Responses of African bovids to Pliocene climatic change. Paleobiology 27, 1–47 (2001).Article 

Google Scholar 
Strömberg, C. A. E. Decoupled taxonomic radiation and ecological expansion of open-habitat grasses in the Cenozoic of North America. Proc. Natl Acad. Sci. USA 102, 11980–11984 (2005).Article 
PubMed 
PubMed Central 
ADS 

Google Scholar 
Kaya, F. et al. The rise and fall of the Old World savannah fauna and the origins of the African savannah biome. Nat. Ecol. Evol. 2, 241–246 (2017).Article 

Google Scholar 
Gravilets, S. & Losos, J. B. Adaptive radiation: contrasting theory with data. Science 323, 732–737 (2009).Article 
ADS 

Google Scholar 
Moen, D. & Morlon, H. Why does diversification slow down? Trends Ecol. Evol. 29, 190–197 (2014).Article 
PubMed 

Google Scholar 
Couvreur, T. L. P. et al. Tectonics, climate and the diversification of the tropical African terrestrial flora and fauna. Biol. Rev. 96, 16–51 (2020).Article 
PubMed 

Google Scholar 
Fontoura, E., Darival Ferreira, J., Bubadué, J., Ribeiro, A. M. & Kerber, L. Virtual brain endocast of Antifer (Mammalia: Cervidae), an extinct large cervid from South America. J. Morphol. 281, 1–18 (2020).Article 

Google Scholar 
Trauth M. A. et al. Recurring types of variability and transitions in the ∼620 kyr record of climate change from the Chew Bahir basin, southern Ethiopia Quaternary. Sci. Rev. https://doi.org/10.1016/j.quascirev.2020.106777 (2021).Janis, C. M. & Manning, E. Antilocapridae. In Evolution of tertiary mammals of North America (eds Janis, C. M., Scott, K. M. & Jacobs, L. L.) 491–507 (Cambridge University Press, 1998).Klimova, A., Munguia-Vega, A., Hoffman, J. I. & Culver, M. Genetic diversity and demography of two endangered captive pronghorn subspecies from the Sonoran Desert. J. Mammal. 95, 1263–1277 (2014).Article 

Google Scholar 
Evin, A., et al. Size and shape of the semicircular canal of the inner ear: A new marker of pig domestication? J. Exp. Zool. B Mol. Dev. Evol. https://doi.org/10.1002/jez.b.23127 (2022).Sánchez, I. M., Cantalapiedra, J. L., Ríos, M., Quiralte, V. & Morales, J. Systematics and evolution of the Miocene three-horned Palaeomerycid ruminants (Mammalia, Cetartiodactyla). PLoS ONE 10, e0143034 (2015).Article 
PubMed 
PubMed Central 

Google Scholar 
Wiley, D. Landmark Editor 3.6 (Institute for Data Analysis and Visualization, Davis, CA, University of California, 2006).R Core Team. R: A language and environment for statistical computing (R Foundation for Statistical Computing, Vienna, Austria, 2022). https://www.R-project.org/.Gunz, P. & Mitteroecker, P. Semilandmarks: a method for quantifying curves and surfaces. Hystrix 24, 103–109 (2013).
Google Scholar 
Adams, D. C. & Otárola-Castillo, E. geomorph: an R package for the collection and analysis of geometric morphometric shape data. Methods Ecol. Evol. 4, 393–399 (2013).Article 

Google Scholar 
Adams, D. C., Collyer, M. L., Kaliontzopoulou, A. geomorph: software for geometric morphometric analyses. R package version 3.2.1 software (2020).Gunz, P., Mitteroecker, P., Bookstein, F. L. Semilandmarks in three dimensions. In Modern morphometrics in physical anthropology. Springer, pp. 73–98 (2005).Maddison, W. P., Maddison, D. R. Mesquite: a modular system for evolutionary analysis. Version 3.04. (2010).Gromolard, C. & Guérin, C. Mise au point sur Parabos cordieri (de Christol), un Bovidé (Mammalia, Artiodactyla) du Pliocène d’Europe occidentale. Géobios 13, 741–755 (1980).Article 

Google Scholar 
Duvernois, M.-P. Mise au point sur le genre Leptobos (Mammalia, Artiodactyla, Bovidae); implications biostratigraphiques et phylogénétiques. Géobios 25, 155–166 (1992).Article 

Google Scholar 
Janis, C. M., Manning, E. Dromomerycidae. In Evolution of Tertiary mammals of North America Volume1: Terrestrial carnivores, ungulates, and ungulatelike mammals (eds. Janis, C. M., Scott, K. M., Jacobs L. L.) 477–490 (Cambridge University Press, 1998).Birungi, J. & Arctander, P. Molecular systematics and phylogeny of the reduncini (artiodactyla: bovidae) inferred from the analysis of mitochondrial cytochrome b gene sequences. J. Mamm. Evol. 8, 125–147 (2001).Article 

Google Scholar 
Lalueza-Fox, C. et al. Molecular dating of caprines using ancient DNA sequences of Myotragus balearicus, an extinct endemic Balear mammal. BMC Evol. Biol. 5, 1–11 (2005).Article 

Google Scholar 
Marot, J. D. Molecular phylogeny of terrestrial artiodactyls, conflict and resolution. In The evolution of artiodactyls (eds Prothero, D. R., Foss, S. C.) 4–18 (The Johns Hopkins University Press, 2007).Webb, D. S. Hornless ruminants. In Evolution of Tertiary mammals of North America Volume1: Terrestrial carnivores, ungulates, and ungulatelike mammals (eds Janis, C. M., Scott, K. M., Jacobs, L. L.) 463–476 (Cambridge University Press, 1998).Mennecart, B. & Métais, G. Mosaicomeryx gen. nov., a ruminant mammal from the Oligocene of Europe and the significance of ‘gelocids’. J. Syst. Palaeontol. 13, 581–600 (2015).Article 

Google Scholar 
Sánchez, I. M., DeMiguel, D., Quiralte, V. & Morales, J. The first known Asian Hispanomeryx (Mammalia, Ruminantia, Moschidae.). J. Vert. Paleontolo. 31, 1397–1403 (2011).Heckeberg, N. S., Erpenbeck, D., Wörheide, G. & Rössner, G. Systematic relationships of five newly sequenced cervid species. PeerJ 4, e2307 (2016).Article 
PubMed 
PubMed Central 

Google Scholar 
Ríos, M., Sánchez, I. M. & Morales, J. A new giraffid (Mammalia, Ruminantia, Pecora) from the late Miocene of Spain, and the evolution of the sivathere-samothere lineage. PLoS ONE 12, e0185378 (2017).Article 
PubMed 
PubMed Central 

Google Scholar 
Vislobokova, I. New data on late Miocene mammals of Kohfidisch, Austria. Paleontol. J. 41, 451–460 (2007).Article 

Google Scholar 
Aiglstorfer, M., Rössner, G. E. & Böhme, M. Dorcatherium naui and pecoran ruminants from the late Middle Miocene Gratkorn locality (Austria). Palaeobiodivers. Palaeoenviron. 94, 83–123 (2014).Article 

Google Scholar 
Janis, C. M. & Scott, K. M. The interrelationships of higher ruminant families with special emphasis on the members of the Cervoidea. Am. Mus. Novit. 2893, 1–85 (1987).
Google Scholar 
Leinders, J. Hoplitomerycidae fam. nov. (Ruminantia, Mammalia) from Neogene fissure fillings in Gargano (Italy). Scr. Geol. 70, 1–68 (1984).
Google Scholar 
Hassanin, A. & Douzery, E. Molecular and morphological phylogenies of Ruminantia, and the alternative position of the Moschidae. Syst. Biol. 52, 206–228 (2003).Article 
PubMed 

Google Scholar 
Métais, G. & Vislobokova, I. Basal ruminants. In The evolution of artiodactyls (eds Prothero, D. R. & Foss, S. C.) 189–212 (The Johns Hopkins University Press, 2007).Mennecart, B., Zoboli, D., Costeur, L. & Pillola, G. L. On the systematic position of the oldest insular ruminant Sardomeryx oschiriensis (Mammalia, Ruminantia) and the early evolution of the Giraffomorpha. J. Syst. Palaeontol. 17, 691–704 (2019).Article 

Google Scholar 
Aiglstorfer, M. et al. Musk Deer on the Run – Dispersal of Miocene Moschidae in the Context of Environmental Changes. In Evolution of Cenozoic land mammal faunas and ecosystems: 25 years of the NOW database of fossil mammals. (eds Casanovas-Vilar, I., van den Hoek Ostende, L. W., Janis, C. M. & Saarinen J.) (Cham: Springer, in press).Klingenberg, C. P. MorphoJ: an integrated software package for geometric morphometrics. Mol. Ecol. Resour. 11, 353–357 (2011).Article 
PubMed 

Google Scholar 
Schlager, S. Morpho and Rvcg – Shape analysis in R. In Zheng, G., Li, S., Szekely, G. Statistical shape and deformation analysis, 217–256 (MA: Academic Press, 2017).Klingenberg, C. P. & Gidaszewski, N. A. Testing and quantifying phylogenetic signals and homoplasy in morphometric data. Syst. Biol. 59, 245–261 (2010).Article 
CAS 
PubMed 

Google Scholar 
Marriott, F. H. C. Barnard’s monte carlo tests: how many simulations? Appl. Stat. 28, 75–77 (1979).Article 

Google Scholar 
Edgington, E. S. Randomization tests (Marcel Dekker, 1987).Tzeng, T. D. & Yeh, S. Y. Permutation tests for difference between two multivariate allometric patterns. Zool. Stud. 38, 10–18 (1999).
Google Scholar 
Revell, L. J. phytools: an R package for phylogenetic comparative biology (and other things). Methods Ecol. Evol. 3, 217–223 (2012).Article 

Google Scholar 
Renaud, S., Dufour, A.-B., Hardouin, E. A., Ledevin, R. & Auffray, C. Once upon multivariate analyses: when they tell several stories about biological evolution. PLoS ONE 10, e0132801 (2015).Article 
PubMed 
PubMed Central 

Google Scholar 
Mitteroecker, P. & Bookstein, F. Linear discrimination, ordination, and the visualization of selection gradients in modern morphometrics. Evol. Biol. 38, 100–114 (2011).Article 

Google Scholar 
Raia, P., Castiglione, S., Serio, C., Mondanaro, A. & Raia, M. P. Package ‘RRphylo’. CRAN Repos. 4, 1–31 (2018).
Google Scholar 
Castiglione, S. et al. A new method for testing evolutionary rate variation and shifts in phenotypic evolution. Methods Ecol. Evol. 9, 974–983 (2018).Article 

Google Scholar 
Morlon, H. et al. “RPANDA: an R package for macroevolutionary analyses on phylogenetic trees.”. Methods Ecol. Evol. 7, 589–597 (2016).Article 

Google Scholar 
Costeur, L., Mennecart, B., Müller, B., Schulz, G. Observations on the scaling relationship between bony labyrinth, skull size and body mass in ruminants. Proc. SPIE 11113, https://doi.org/10.1117/12.2530702 (2019).Costeur, L., Mennecart, B., Müller, B. & Schulz, G. Prenatal growth stages show the development of the ruminant bony labyrinth and petrosal bone. J. Anat. 230, 347–353 (2017).Article 
PubMed 

Google Scholar 
Mennecart, B. & Costeur, L. Shape variation and ontogeny of the ruminant bony labyrinth, an example in Tragulidae. J. Anat. 229, 422–435 (2016).Article 
PubMed 
PubMed Central 

Google Scholar 
Clauss, M., Steuer, P., Müller, D. W. H., Codron, D. & Hummel, J. Herbivory and body size: allometries of diet quality and gastrointestinal physiology, and implications for herbivore ecology and dinosaur gigantism. PLoS One 8, e68714 (2013).Article 
CAS 
PubMed 
PubMed Central 
ADS 

Google Scholar 
du Toit, J. T. & Owen-Smith, N. Body size, population metabolism, and habitat specialization among large African herbivores. Am. Nat. 133, 736–740 (1989).Article 

Google Scholar 
Mennecart B., Becker D., & Berger J. -P. Mandible shape of ruminants: between phylogeny and feeding habits. In: Ruminants: Anatomy, behavior, and diseases, (ed. Mendes R. E.) 205–226 (Nova Science Publishers, 2012).Bokma, F. et al. Testing for Depéret’s rule (body size increase) in mammals using combined extinct and extant data. Syst. Biol. 65, 98–108 (2016).Article 
PubMed 

Google Scholar 
Besiou, E., Choupa, M. N., Lyras, G. & van der Geer, A. Body mass divergence in sympatric deer species of Pleistocene Crete (Greece). Palaeontol. Electron. 25, a23 (2022).
Google Scholar 
Mennecart B., Métais G., Tissier J., Rössner G. E., & Costeur L. 3D models related to the publication: Reassessment of the enigmatic ruminant Miocene genus Amphimoschus Bourgeois, 1873 (Mammalia, Artiodactyla, Ruminantia, Pecora). MorphoMuseuM 7, e131 (2021).Mennecart, B., Perthuis de, A. D. & Costeur, L. 3D models related to the publication: The first French tragulid skull (Mammalia, Ruminantia, Tragulidae) and associated tragulid remains from the Middle Miocene of Contres (Loir-et-Cher, France). MorphoMuseuM 3, e4 (2018).Article 

Google Scholar 
Aiglstorfer, M., Costeur, L., Mennecart, B. & Heizmann, E. P. J. Micromeryx? eiselei – a new moschid species from Steinheim am Albuch, Germany, and the first comprehensive description of moschid cranial material from the Miocene of Central Europe. MorphoMuseuM 3, e4 (2107).Article 

Google Scholar 
Costeur, L. & Mennecart, B. 3D models related to the publication: Prenatal growth stages show the development of the ruminant bony labyrinth and petrosal bone. MorphoMuseuM 2, e3 (2016).Article 

Google Scholar 
Mennecart, B. & Costeur, L. 3D models related to the publication: a Dorcatherium (Mammalia, Ruminantia, Middle Miocene) petrosal bone and the tragulid ear region. MorphoMuseuM 2, e2 (2016).Article 

Google Scholar 
Mennecart, B. et al. Allometric and phylogenetic aspects of stapes morphology in ruminantia (Mammalia, Artiodactyla). Front. Earth Sci. 8, 176 (2020). More

The rise of animals (metazoans) is a seminal event in the history of life. The Cambrian Radiation ~540 Ma marks the appearance of abundant and diverse metazoans and increasing ecosystem complexity in the fossil record1. A causal relationship between the redox and fossil records is proposed, where oxygen provision reached a threshold, or series of thresholds, which allowed for the diversification of metazoans with increasing metabolic demands2. Global geochemical data, however, suggest that oxygenation was not a simple, linear process, but rather occurred episodically via a series of short-lived pulses (1–3 Myr), or ‘oceanic oxygenation events’ (OOEs)3,4. Early and even later Cambrian seas likely had shallower, and more dynamic, oxygen minimum zones (OMZs) than modern oceans5,6. Such pulses of increased oxygenation (and related changes in productivity) are hypothesised to have increased the extent of shallow-ocean oxygenation and hence to have promoted diversification7. But what remains unquantified is the community-wide response of metazoans to such redox cycles, an insight into the evolutionary processes involved, and hence whether these pulses were indeed a driving force for the Cambrian Radiation.In order to test the hypothesis that oxic pulses led to diversification and potentially ecological development, a correlation between increased oxygenation, rates of origination, and metrics of metazoan ecosystem complexity needs to be demonstrated. Early Cambrian marine environments were heterogeneous with respect to oxygen provision and nutrient load at a regional scale, so in order to investigate potential correlations, we require the integration of global and local redox proxies, and biotic records in the same stratigraphically well-constrained geological successions.During the early Cambrian, the Siberian Platform was a vast isolated, tropical continent almost entirely covered by an epicontinental sea (Fig. 1a)8,9. The platform supported a single metacommunity, i.e. a species pool with many local, interacting communities e.g.10, representing a third of total early Cambrian metazoan benthic diversity with widespread metazoan (archaeocyath sponge) reefs that formed bioherms (Fig. 1b)7,11. Dynamic and synchronous changes of body size in archaeocyath sponges, hyoliths, and helcionelloid molluscs through the early Cambrian on the Siberian Platform have been quantified, which coincide with elevated biodiversity and rates of origination: these have been proposed to follow OOEs12. Here we consider temporal changes in both the position of archaeocyath sponge reefs as a function of relative water depth, and in individual reef size (diameter), as well as the ecological complexity of the reef-building and dwelling communities by quantification of changing reef community membership of sessile archaeocyath sponge, coralomorph, and cribricyath species, on the Siberian Platform.Fig. 1: Palaeogeographic and stratigraphic position of the early Cambrian archaeocyath reefs of the Lena-Aldan area on the Siberian Platform.a Early Cambrian palaeofacies zonation map of the Siberian Platform. b Cross section to show relative positions of sampled transects along the Lena River11,40,66,67,68. c Lithostratigraphy, biostratigraphy, carbon isotope (δ13C)29,31,32 and carbonate-associated sulfate sulfur isotope (δ34SCAS)7 data for sections from the middle Lena River (Isit’, Zhurinsky Mys, Achchagy-Kyyry-Taas, and Achchagy-Tuoydakh). S.E.—Sinsk Event; Tolb.—Tolba Formation; ATD., BOT., N.-D., TOM.—Atdabanian, Botoman, Nemakit-Daldynian, and Tommotian local stages, respectively.Full size imageTo quantify ecological complexity, we used metacommunity analyses, which compare the structure between communities in terms of taxa (generally species) compositions spatially and temporally10 (see Methods). The ‘Elements of Metacommunity Structure’ framework used here is a hierarchical analysis that identifies properties in site-by-species presence/absence matrices that are related to the underlying processes, such as species interactions, dispersal, and environmental filtering that shape species distributions10. Application to various marine and terrestrial palaeocommunities has demonstrated the robustness of these methods to fossil data and sample size variations13,14. There are fourteen different types of metacommunity structure which are determined by the calculation of three metacommunity metrics: Coherence, Turnover, and Boundary Clumping, which reveal different controlling processes of underlying metacommunity structure10,15,16,17,18.The most ecologically complex metacommunities are classified as Clementsian, and have positive coherence, turnover and boundary clumping16. Clementsian metacommunities contain groups of taxa with similar range boundaries that respond to the environment synchronously as taxa have physiological or evolutionary trade-offs within the communities associated with environmental thresholds19. By contrast, when taxa respond individualistically to the underlying environment, without accounting for other taxa within the community, the structure is Gleasonian, and is defined by positive coherence and turnover but no significant boundary clumping16. When coherence is positive, but turnover is not significantly different from random, then the resultant metacommunity structures are known as quasi-structures (e.g. quasi-Clementsian), which reflect weaker underlying structuring processes.We determined the metacommunity structure for archaeocyath sponge species on the Siberian Platform throughout their early Cambrian record using an entire previously published data set11 then on a sub-set of metacommunities which had a sufficient number of reef sites to be suitable for analyses, i.e. with a sufficient number of sites to be statistically significant. Further, to investigate the effects of water depth on metacommunity structure, we used Spearman rank correlations to test whether the metacommunity ranking (as determined by reciprocal averaging, a type of correspondence analysis which ordinates the sites based on their species composition17), is significantly correlated to water depth. Finally, to quantify how pairwise associations between taxa change between the three temporally different metacommunities, we determined which pairwise taxa co-occurrences are significantly non-random using a combinatorics approach, and whether any non-random co-occurrences are positive or negative20.Species richness estimates are highly sensitive to differences in sampling. When comparing species richness of assemblages from several time intervals, it is advisable to standardise sampling across those assemblages to ensure that changes in species richness are not attributable to sampling differences. One approach is to subsample each time interval down to a standardised number of individuals (size-based rarefaction), but this approach can underestimate changes in richness because it tends to sample low-richness assemblages more completely than high-richness ones21. Coverage-based rarefaction, where each sample is down-sampled to a standardised level of taxonomic completeness, avoids this potential issue. The coverage of a sample is the proportion of species in the assemblage which are represented in that sample, and it can be estimated by subtracting the proportion of singletons in a sample from 1 (e.g.22; see also21 for details). We used the estimateD function from R package iNEXT23 to produce coverage-standardised species richness estimates with 95% confidence intervals, by repeatedly down-sampling the sampled assemblage from each time interval to match the coverage of the lowest-coverage interval. We did this by setting datatype = “abundance”, base = “coverage” and leaving all other arguments as default.In sum, we test the biotic response to OOEs by compiling metrics of archaeocyath reef size, location, and metacommunity complexity, integrated with existing data on archaeocyath individual size, species richness and origination and extinction rates12 and high-resolution geochemistry4,7 recalculated to the same stratigraphic scale, on the Siberian Platform over 11 Myr through Cambrian stages 2–3 (mid-Tommotian to early Botoman on the Siberian stratigraphic scale; 525–514 Ma). These results are used to quantify the community-wide response of metazoans to extrinsic redox cycles, and hence gain insight into the evolutionary processes involved.Geological setting and evolution of redoxDuring the early Cambrian shallow marine carbonates associated with evaporites and siliciclastics dominated the inner Siberian Platform, passing to shallow marginal carbonates of transitional facies known as the transitional zone (or the Anabar-Sinsk), thence to deep ramp and slope settings that accumulated organic-rich limestone and shale (Fig. 1a)24,25,26. Archaeocyathan reefs or bioherms were almost entirely restricted to the transitional facies. Such reefs appeared and proliferated during Cambrian stages 2 and 3 (Tommotian, Atdabanian and earliest Botoman), disappeared at the beginning of Stage 4 (middle Botoman) and re-appeared briefly at the end of this stage (Toyonian).We integrate palaeontological (archaeocyath species number and individual size), palaeoecological (reef size and palaeodepth location) and chemostratigraphic information (carbon isotope cycles 5p, 6p, and II–VII) for sections of the Aldan, Selinde and Lena rivers with sub-metre-scale lithostratigraphic subdivisions27,28,29,30,31,32,33 (Figs. 1c, 2a–c, 3a). This results in negligible uncertainty associated with sample heights, which are fixed relative to a consistent datum within each section.Fig. 2: Lithostratigraphy, biostratigraphy and carbon isotope (δ13C) data for sections of the Aldan and Selinde rivers bearing the earlierst archaeocyath reef communities of the Siberian Platform.a Dvortsy27,28,30 b Ulakhan-Sulugur33,34, and c Selinde69,70.Full size imageFig. 3: Summary of geochemical and biotic changes through the early Cambrian, Siberian Platform, and uranium isotope data representing a global record.a International and Siberian timescale, within age model C of 57. ND—Nemakit-Daldynian regional stage; U’-Y—Ust’-Yudoma Formation. b Summary of carbon and sulphur isotopes (from the Lena River, Siberia7). c Uranium isotopes from Siberia (grey; Sukharikha and Bol’shaya Kuonamka rivers), South China (blue), and Morocco (orange) (all data points are larger than 2SE)4. d Archaeocyath sponge species diversity and maximum diameter12. Plotted richness values are the species richness estimator21 with accompanying 95% confidence interval, calculated using the estimated function from R package iNEXT62. e Rates of archaeocyath sponge species origination and extinction12. f Reef location as a function of relative water depth (Supplementary Table 1). FWWB—Fair weather wave base. SWB—Storm weather wave base. g Reef/bioherm diameter, coloured by relative water depth (see column f, and Supplementary Table 2). h Number of reef community types (Supplementary Table 3). i Archaeocyath reef ecosystem complexity, with percentage of species co-occurrence as changing proportions of total non-random and positive and negative. G = Gleasonian, QG = Quasi-Gleasonian, C = Clementsian.Full size imageThroughout Cambrian stages 2 and 3, high-amplitude positive δ13C carbon isotope excursions show a strong positive covariation with the sulphur isotope composition of carbonate-associated sulphate (δ34SCAS) in sections from the Lena River (Fig. 3b)7. The rising limbs of these excursions are interpreted as intervals of progressive burial of reductants under anoxic bottom water conditions, and a progressive increase in atmospheric oxygen7. Coincident δ13C and δ34SCAS peaks (numbered II–VII) correspond with a pulse of atmospheric oxygen into the shallow marine environment (creating an OOE), followed by a corresponding decrease in reductant burial under more widespread marine oxia (falling limbs of δ13C and δ34SCAS), and leading to gradual de-oxygenation over Myr7. In addition, phosphorous retention might have occurred under oxic shallow marine conditions, acting to reduce primary productivity and further oxygenate the shallow marine environment in the short-term ( More

Fossil record of insectsWe compiled all species-level fossil occurrences of insects using https://paleobiodb.org/ (PBDB) as a starting point (downloaded October 12, 2021). The dataset obtained from PBDB contained initially 5808 occurrences for a period ranging from the Asselian to the Rhaetian. The dataset was cleaned of synonyms, outdated combinations, nomina dubia, and other erroneous and doubtful records, based on revisions provided in the literature and/or on the expertise of the authors. After correction, including data addition from the literature, our dataset was composed of 3636 species (1784 genera, and 418 families) for 17,250 occurrences resulting from an in-depth study and curation of the entire bibliography of fossil insects, spanning from the Asselian (lowermost Permian) to the Rhaetian (uppermost Triassic). Although most of the taxa included in the datasets are nominal taxa (published and named), a few unnamed taxa (genera or species) that are considered separate from others were also included, although not formally named in the literature or not published yet. These unpublished taxa are identifiable by the notation ‘fam. nov.’ or ‘gen. nov.’ following their names.Occurrences used here are specimens originating from a given stratigraphic horizon assigned to a given taxon. The age of each occurrence is based on data from PBDB, corrected with a more precise age (generally stage, sometimes substage), and the age of each time bin boundaries relies on the stratigraphic framework proposed in the International Chronostratigraphic Chart (updated to correspond with the ICS 2022/0295). Similarly, the ages of some species assigned to the wrong stage were corrected. In fact, some species from the French Permian deposit of Lodève were initially considered to be of Artinskian age in PBDB but most species from this deposit originate from the Merifons member, which is of Kungurian age96.Our data compilation allows a robust integration of data before and after our period of interest (i.e. the lower Permian and all geologic stages after the Carnian) to encompass occurrences of genera that may survive until the Late Triassic and to generate a sufficient background for the model to correctly estimate the extinction events around the P/T boundary. Since we used different datasets, the differences between genus-level or family-level occurrence numbers are explained by the systematic placement of some specimens that can only be placed confidently in a family but not in a genus (Supplementary Table 1). Tentative species identifications originally placed with uncertainty (reported as ‘aff.’ or ‘?’) were always included at a higher taxonomic level. Uncertain generic attributions were integrated as occurrences at the family level (e.g. a fossil initially considered Tupus? is recorded as an occurrence of Meganeuridae). Our total dataset was subdivided into smaller datasets, which represent orders or other subclades of insects (e.g. Mecoptera, Holometabola and Polyneoptera). Note that all the ichnospecies—a species name assigned to trace fossils (e.g. resting trace, nest and leaf damage)— and insect eggs (e.g. Clavapartus latus, Furcapartus exilis and Monilipartus tenuis) were not included in the analyses97. To prevent potential issues regarding the diversification estimates for clades with poor delineation, we refrained from analysing several orders that serve as taxonomic ‘wastebaskets’ (e.g. Grylloblattodea). These groups are poorly defined, likely polyphyletic or paraphyletic, and not supported by apomorphic characters—e.g. the monophyly of the ‘Grylloblattodea’ (Grylloblattida Walker, 1914 plus numerous fossil families and genera of uncertain affinities) is not supported by any synapomorphy, nor the relationships within this group. The occurrences assigned to these orders were rather included in analyses conducted at a higher taxonomic level (at the Polyneoptera level in the case of the ‘Grylloblattodea’). The detail of the composition of all the datasets is given in Supplementary Table 14, and each dataset is available in Supplementary Data 1.Studying extinction should, when possible, rely on species-level diversity to better circumscribe extinction events at this taxonomic rank, which is primarily affected by extinction98,99,100. However, in palaeoentomology, species-level occurrence data may contain less information than genus-level data, mainly because species are most of the time only known from one deposit, resulting in reduced life span, and are also sometimes poorly defined. Insects are also less prone to long-lasting genera or species than other lineages, maybe because of the relatively short time between generations (allowing for rapid evolution) or because morphological characters are better preserved or more diagnostic than in other lineages (i.e. wing venation), allowing easier differentiation. Another argument for the use of genus-level datasets is the possibility to add occurrences represented by fossils that cannot be assigned at the species level because of poor preservation or an insufficient number of specimens/available characters. By extension, the genus life span provides clues as to survivor taxa and times of origination during periods of post-extinction or recovery. A genus encompassing extinction events indicates that at least one species of this genus crossed the extinction. To get the best signal and infer a robust pattern of insect dynamics around the P/T events, we have chosen to analyse our dataset at different taxonomic ranks (e.g. genus, family and order levels) to extract as much evidence as possible.To further support our choice to work at these different levels, most recent works aiming to decipher the diversification and extinction in insect lineages have worked using a combination of analyses21,22,26; this also applies to non-insect clades51,101,102. This multi-level approach should maximise our understanding of the Permo–Triassic events.Assessing optimal parameters and preliminary testsPrior to choosing the settings for the final analyses (see detail in Dynamics of origination and extinction), a series of tests were carried out to better evaluate the convergence of our analyses. First, we analysed our genus-level dataset with PyRate36 running for 10 million generations and sampling every 10,000 generations, on ten randomly replicated datasets using the reversible-jump Markov Chain Monte Carlo (RJMCMC) model37 and the parameters of PyRate set by default. As the convergence was too low, new settings were used, notably increasing the number of generations to 50 million generations and monitoring the MCMC mixing and effective sample size (ESS) each 10 million generations. We modified the minimal interval between two shifts (-min_dt option, testing 0.5, 1.5 and 2), and found no major difference in diversification patterns between our tests. We have opted for 50 million generations with a predefined time frame set for bins corresponding to the Permian and Triassic stages, and a minimum interval between two shifts of two Ma. These parameters allow for maintaining a short bin frame and high convergence values while correctly identifying periods of diversification and extinction. For each analysis, ten datasets were generated using the extract.ages function to randomly resample the age of fossil occurrences within their respective temporal ranges (i.e. resampled ages are randomly drawn between the minimum and the maximum ages of the geological stratum). We monitored chain mixing and ESS by examining the log files in Tracer 1.7.1103 after excluding the first 10% of the samples as a burn-in period. The parameters are considered convergent when their ESS are greater than 200.Dynamics of origination and extinctionWe carried out the analyses of the fossil datasets based on the Bayesian framework implemented in the programme PyRate36. We analysed the fossil datasets under two models: the birth–death model with constrained shifts (BDCS38) and the RJMCMC (-A 4 option37). These models allow for a simultaneous estimate for each taxon: (1) the parameters of the preservation process (Supplementary Fig. 17), (2) the times of origination (Ts) and extinction (Te) of each taxon, (3) the origination and extinction rates and their variation through time for each stage and (4) the number and magnitude of shifts in origination and extinction rates.All analyses were set with the best-fit preservation process after comparing (-PPmodeltest option) the homogeneous Poisson process (-mHPP option), the non-homogeneous Poisson process (default option), and the time-variable Poisson process (-qShift option). The preservation process infers the individual origination and extinction times of each taxon based on all fossil occurrences and on an estimated preservation rate, denoted q, expressed as expected occurrences per taxon per Ma. The time-variable Poisson process assumes that preservation rates are constant within a predefined time frame but may vary over time (here, set for bins corresponding to stages). This model is thus appropriate when rates over time are heterogeneous.We ran PyRate for 50 million MCMC generations and a sampling every 50,000 generations for the BDCS and RJMCMC models with time bins corresponding to Permian and Triassic stages (-fixShift option). All analyses were set with a time-variable Poisson process (-qShift option) of preservation and accounted for varying preservation rates across taxa using the Gamma model (-mG option), that is, with gamma-distributed rate heterogeneity with four rate categories36. As explained above, the minimal interval between two shifts (-min_dt option) was modified and a value of 2 was used. The default prior to the vector of preservation rates is a single gamma distribution with shape = 1.5 and rate = 1.5. We reduced the subjectivity of this parameter, and favoured a better adequation to the data, allowing PyRate to estimate the rate parameter of the prior from the data by setting the rate parameter to 0 (-pP option). Therefore, PyRate assigns a vague exponential hyper-prior to the rate and samples the rate along with all other model parameters. Similarly, because our dataset does not encompass the entire fossil record of insects, we assumed that a possible edge effect may interfere with our analyses, with a strong diversification during the lowermost Permian and, conversely a strong extinction during the uppermost Triassic. Because the RJMCMC and BDCS algorithms look for rate shifts, we constrained the algorithm to only search for shifts (-edgeShift option) within the following time range 295.0 to 204.5 Ma. We monitored chain mixing and ESS by examining the log files in Tracer 1.7.1103 after excluding the first 10% of the samples as a burn-in period. The parameters are considered convergent when their ESS are greater than 200.We then combined the posterior estimates of the origination and extinction rates across all replicates to generate rates through-time plots (origination, extinction, and net diversification). Shifts of diversification were considered significant when log Bayes factors were >6 in the RJMCMC model, while we considered shifts to be significant in the BDCS model when mean rates in a time bin did not overlap with the 95% credibility interval (CI) of the rates of adjacent time bins.We replicated all the analyses on ten randomly generated datasets of each clade and calculated estimates of the Ts and the Te as the average of the posterior samples from each replicate. Thus, we obtained ten posterior estimates of the Ts and Te for all taxa and we used these values to estimate the past diversity dynamics by calculating the number of living taxa at each time point. For all the subsequent analyses, we used the estimated Ts and Te of all taxa to test whether or not the origination and the extinction rate dynamics were correlated with particular abiotic factors, as suggested by the drastic changes in environmental conditions known during the Permo–Triassic. We used proxies for abiotic factors, such as global continental fragmentation or the dynamic of major clades of plants, and for biotic factors via species interaction within and between ecological guilds. This approach avoids re-modelling preservation and re-estimating times of origination and extinction, which reduces drastically the computational burden, while still allowing to account for the preservation process and the uncertainties associated with fossil ages. Similarly, the times of origination and extinction used in all the subsequent analyses were obtained while accounting for the heterogeneity of preservation, origination and extinction rates. To discuss the magnitude of the periods of extinction and diversification, we compared the magnitude of these events to the background origination and extinction rates (i.e. not during extinction or diversification peaks).The PyRate approach has proven to be robust following a series of tests and simulations that reflect commonly observed biases when modelling past diversity dynamics31,38. These simulations were based on datasets simulated under a range of potential biases (i.e. violations of the sampling assumptions, variable preservation rates, and incomplete taxon sampling) and reflecting the limitations of the fossil record. Simulation results showed that PyRate is able to correctly estimate the dynamics of origination and extinction rates, including sudden rate changes and mass extinction, even if the preservation levels are low (down to 1–3 fossil occurrences per species on average), the taxon sampling is partial (up to 80% missing) or if the datasets have a high proportion of singletons (exceeding 30% of the taxa in some cases). The strongest bias in birth–death rate estimates is caused by incomplete data (i.e. missing lineages) altering the distribution of taxa; a pervasive effect often mentioned for phylogeny-based models104,105,106. However, in the case of PyRate, the simulations confirm the absence of consistent biases due to an incomplete fossil record36. Finally, the recently implemented RJMCMC model was shown to be very accurate for estimating origination and extinction rates (i.e. more accurate than the BDCS model, the boundary-crossing and three-time methods) and is able to recover sudden extinction events regardless of the biases in the fossil dataset37.The severity of extinctions and survivorsFor each event—the Roadian–Wordian, the LPME, and the Ladinian–Carnian—we quantified the percentage of extinctions and survivors at the genus level. We used the Te and Ts from our RJMCMC analysis and computed the mean for the Te (Tem) and for the Ts (Tsm) of each genus. We then filtered our dataset to keep only the genera with a Tsm older than the upper boundary of the focal event, i.e., we only kept the genera that appeared before the end of the event. Then, we discarded the genera that have disappeared before the lower boundary of the focal event, i.e. Tem older that the lower boundary of the event. The remaining genera, which corresponds to all the genera (total) present during the crisis (Ttgen), can be classified into two categories, ‘survivor genera’ (Sgen), i.e. those that survived the crisis, and those that died: ‘extinct genera’ (Egen). The survivors have a Tem younger than the upper boundary of the focal event, while the ‘extinct genera’ died out during the event and have a Tem between the lower and upper boundaries of the event of interest. To obtain the percentage of survivors, we used the following formula: (Sgen/Ttgen) × 100. Similarly, the percentage of extinction is calculated as: (Egen/Ttgen) × 100.Age-dependent extinction modelWe assessed the effect of taxon age on the extinction probability by fitting the age-dependent extinction (ADE; -ADE 1 option) model50. This model estimates the probability for a lineage to become extinct as a function of its age, also named longevity, which is the elapsed time since its origination. It is recommended to run the ADE model over time windows with roughly constant origination and extinction rates, as convergence is difficult—but not impossible—to reach in extinction or diversification contexts50. We ran PyRate for 50 million MCMC generations with a sampling every 50,000 generations, with a time-variable Poisson process of preservation (-qShift option), while accounting for varying preservation rates across taxa using the Gamma model (-mG option). We replicated the analyses on ten randomised datasets and combined the posterior estimates across all replicates. We estimated the shape (Φ) and scale (Ψ) parameters of the Weibull distribution, and the taxon longevity in a million years. According to ref. 50, there is no evidence of age-dependent extinction rates if Φ = 1. However, the extinction rate is higher for young species and decreases with species age if Φ  1. Although ADE models are prone to high error rates when origination and extinction rates increase or decrease through time, simulations with PyRate have shown that fossil-based inferences are robust50. We investigated the effect of ADE during three different periods (-filter option) as follows: (1) between 264.28 Ma and 255 Ma (pre-decline), (2) between 254.5 Ma and 251.5 Ma (decline) and (3) between 234 Ma and 212 Ma (post-crisis). We monitored chain mixing and ESS by examining the log files in Tracer 1.7.1103 and considered the convergence of parameters sufficient when their ESS were greater than 200.Selection of abiotic and biotic variablesTo test correlations of insect diversification with environmental changes, we examined the link between a series of environmental variables and origination/extinction rates over a period encompassing the GEE, the LPME and the CPE but also for each extinction event. We focused on the role of nine variables, also called proxies, which have been demonstrated or assumed to be linked to extinctions and changes in insect diversity26,67.The variations in the atmospheric CO2 and O2 concentrations are thought to be correlated with the diversification of several insect lineages, including the charismatic giant Meganeuridae65,66,67. Because the increase of O2 concentration has likely driven the diversification of some insects, its diminution may have resulted in the extinction or decline of some lineages. Therefore, we investigated the potential correlation of the variations of this variable with insect dynamics using data from ref. 55. We extracted the data, with 1-million-year time intervals, spanning the Permo–Triassic.Similarly, the modification of CO2 concentration, notably its increase, is known to promote speciation in some modern insect groups107. Therefore, a similar effect may have occurred during the Permian and Triassic but remains to be tested. We based our analyses on the dataset of ref. 108. We used their cleaned dataset and extracted all verified values for the Permo–Triassic interval. Because the initial data (i.e. independent estimates) were made in various locations for the same age, different values of the CO2 concentration are provided. We incorporated all these values in our analysis, allowing PyRate to search for a correlation for each value of the CO2 concentration. We obtain a final correlation independent of the sampling location, in line with our large-scale analysis.The continental fragmentation, as approximated by plate tectonic change over time, has recently been proposed as a driver of Plecoptera dynamics26. Because the period studied encompasses a major geological event, the fragmentation of the supercontinent Pangea, we investigated the effect of continental fragmentation on insect diversification dynamics. We retrieved the index of continental fragmentation developed by ref. 69 using paleogeographic reconstructions for 1-million-year time intervals. This index approaches 1 when all plates are disjoined (complete plate fragmentation) and approaches 0 when the continental aggregation is maximal.Climate change (variations in warming and cooling periods) is a probable driver of diversification changes over the history of insects21,109. Temperature is likely directly linked with insect dynamics109 but also with their food sources, notably plants110. Because it was demonstrated that modification of temperature impacted floral assemblages110, we tested the correlation between temperature variations and the diversification dynamic of insects. Major trends in global climate change through time are typically estimated from relative proportions of different oxygen isotopes (δ18O) in samples of benthic foraminiferan shells111. We used the data from ref. 112, converted to absolute temperatures following the methodology described in Condamine et al.113 (see their section Global temperature variations through time). The resulting temperature data reflects planetary-scale climatic trends, with time intervals inferior to 1-million-year, which can be expected to have led to temporally coordinated diversification changes in several clades rather than local or seasonal fluctuations.The fluctuation in relative diversity of gymnosperms, non-Polypodiales ferns, Polypodiales ferns, spore-plants, and later the rise of angiosperms has likely driven the diversification of numerous insects57,60,61,114. Close interactions between insects and plants are well-recorded during the Permian and Triassic57,60,61. In fact, herbivorous insects are known to experience high selection pressure from bottom-up forces, resulting from interactions with their hosts or feeding plants30,72. Therefore, it appears crucial to investigate the effect of these modifications on the insects’ past dynamics. We used the data from ref. 38 for the different plant lineages (all with 1-million-year time intervals). All the datasets for these variables are available in the publications cited aside from each variable or in Supplementary Data 1.Multivariate birth–death modelWe used the multivariate birth–death (MBD) model to assess to what extent biotic and abiotic factors can explain temporal variation in origination and extinction rates55. The model is described in ref. 55, where origination and extinction rates can change through time in relation to environmental variables so that origination and extinction rates depend on the temporal variations of each factor. The strength and sign (positive or negative) of the correlations are jointly estimated for each variable. The sign of the correlation parameters indicates the sign of the resulting correlation. When their value is estimated around zero, no correlation is estimated. An MCMC algorithm combined with a horseshoe prior, controlling for over-parameterisation and for the potential effects of multiple testing, jointly estimates the baseline origination (λ0) and extinction (µ0) rates and all correlation parameters (Gλ and Gµ)55. The horseshoe prior is used to discriminate which correlation parameters should be treated as noise (shrunk around 0) and which represent a true signal (i.e. significantly different from 0). In the MBD model, a correlation parameter is estimated to quantify independently the role of each variable on the origination and the extinction.We ran the MBD model using 20 (for short intervals) or 50 million MCMC generations and sampling every 20,000 or 50,000 to approximate the posterior distribution of all parameters (λ0, µ0, nine Gλ, nine Gµ and the shrinkage weights of each correlation parameter, ωG). The MBD analyses used the Ts and the Te derived from our previous analyses under the RJMCMC model. The results of the MBD analyses were summarised by calculating the posterior mean and 95% CI of all correlation parameters and the mean of the respective shrinkage weights (across ten replicates), as well as the mean and 95% CI of the baseline origination and extinction rates. We carried out six analyses, over: (1) the Permo–Triassic (between 298.9 and 201.3 Ma); (2) the Roadian–Wordian (R/W) boundary (between 270 and 265 Ma), (3) the LPME (between 254.5 and 250 Ma), (4) the Ladinian–Carnian (L/C) boundary (between 240 to 234 Ma), (5) the Permian period (between 298.9 and 251.902 Ma) and (6) the Triassic period (between 251.902 and 201.3 Ma). We monitored chain mixing and ESS by examining the log files in Tracer 1.7.1103 and considered the convergence of parameters sufficient when their ESS were greater than 200.Multiple clade diversity-dependence modelTo assess the potential effect of diversity-dependence on the diversity dynamics of three or four insect guilds, we used the multiple clade diversity-dependence (MCDD) model in which origination and extinction rates are correlated with the diversity trajectory of other clades31. This model postulates that competitive interactions linked with an increase in diversity results in decreasing origination rates and/or increasing extinction rates. The MCDD model allows for testing diversity-dependence between genera of a given clade or between genera of distinct clades sharing a similar ecology.We estimated the past diversity dynamics for three (i.e. herbivores, predators, and a guild composed of generalists + detritivores/fungivores dubbed ‘others’) or four insect groups or guilds (i.e. herbivores, predators, generalists and detritivores/fungivores) by calculating the number of living species at every point in time based on the times of origination (Ts) and extinction (Te) estimated under the RJMCMC model (see above) (Supplementary Figs. 19–24). We defined our four insect groups with a cautious approach i.e. insect genera, families or orders for which nothing is known about the ecology or about the ecology of their close relatives were not considered for the analysis. For example, no diet was assigned to Diptera, Mecoptera or Glosselytrodea. The ecology of the Triassic Diptera and Permo–Triassic Mecoptera is difficult to establish because extant Diptera and Mecoptera have a wide diversity of ecology. Fossil Mecoptera are also putatively involved in numerous interactions with plants (species with elongated mouthparts), suggesting a placement in the herbivore group, while other species were likely predators. Therefore, we cannot decide to which group each species belongs. Similarly, nothing is known about the body and mouthparts of the Glosselytrodea, most of the time described based on isolated wings; we did not assign the order to any group. The definition and delineation of insect clades have also challenged the placement of several orders (e.g. ‘Grylloblattodea’) in one of our four groups. The order ‘Grylloblattodea’ is poorly delineated and mostly serves as a taxonomic ‘wastebasket’ to which it is impossible to assign a particular ecology. Finally, genera, species, or families not placed in a higher clade (e.g. Meshemipteron, Perielytridae) were not included in the analysis. Oppositely, the guilds ‘herbivores’ and ‘predators’ are well defined, and their ecology is evidenced by the morphology of their representatives and the principle of actualism. For example, the ecology of Meganeurites gracilipes (Meganeuridae) has been deeply studied, and its enlarged compound eyes, its sturdy mandibles with acute teeth, its tarsi and tibiae bearing strong spines, and the presence of a pronounced thoracic skewness are specialisations today found in dragonflies that capture their prey while in flight115. All Odonatoptera are well-known predator insects. The raptorial forelegs of the representatives of the order Titanoptera and their mouthparts with strong mandibles are linked with predatory habits81. The Palaeodictyopteroidea were herbivorous insects with long, beak-like, piercing mouthparts, and probably a sucking organ81,82. Most Hemiptera are confidently considered herbivorous insects by comparison with their extant representatives. For example, the Cicadomorpha or Sternorrhyncha are known to feed on plants and their fossil representatives likely possessed the same ecology because of similar morphologies116. Some hemipteran families (e.g. Nabidae) are predators and we cautiously distinguished herbivorous and carnivorous taxa among Hemiptera. The detail of the ecological assignations for the 1009 genera included in our analyses can be found in Supplementary Data 1 (Table MCCD).We calculated ten diversity trajectories from the ten replicated analyses under the RJMCMC model. The estimation of past species diversity might be biased by low preservation rates or taxonomic uncertainties. However, such trajectory curves are likely to provide a reasonably accurate representation of the past diversity changes in the studied clades, notably because the preservation during the Permian and Triassic period is relatively good for insects (i.e. no gaps).Our MCDD analyses comprise all the insect genera spanning from the lowermost Permian to the uppermost Triassic and were run and repeated on ten replicates (using the Te and Ts estimated under the RJMCMC model) with 50 million MCMC generations and a sampling frequency of 50,000. For each of the four insect groups, we computed the median and the 95% CI of the baseline origination and extinction rates (λi and µi), the within-group diversity-dependence parameters gλi and gµi, and the between-groups diversity dependence parameters gλij and gµij. The mean of the sampled diversity dependence parameters (e.g. gλij) was used as a measure of the intensity of the negative (if positive) or positive interactions (if negative) between each pair of groups. The interactions were considered significant when their median was different from 0 and the 95% CI did not overlap with 0. We monitored chain mixing and ESS by examining the log files in Tracer 1.7.1103 and considered the convergence of parameters sufficient when their ESS were greater than 200.We cross-validated the result of the MCDD model using the MBD model. The MBD model can be used to run a multiple clade diversity-dependence analysis by providing the diversity trajectories of insect guilds as a continuous variable. These data are directly generated by PyRate using the lineages-through-time generated by the RJMCMC analyses (-ltt option). We ran the MBD model using 50 million MCMC generations and sampling every 50,000 to approximate the posterior distribution of all parameters (λ0, µ0, four Gλ, four Gµ and the shrinkage weights of each correlation parameter, ωG). We carried out three analyses, over the period encompassing the three extinction events (between 275 and 230 Ma): (1) for herbivores; (2) for predators; and (3) for ‘others’. For each analysis, the lineages-through-time data of the two other guilds are used as continuous variables to investigate a diversity dependence effect. We monitored chain mixing and ESS by examining the log files in Tracer 1.7.1103 and considered the convergence of parameters sufficient when their ESS were greater than 200.Reporting summaryFurther information on research design is available in the Nature Portfolio Reporting Summary linked to this article. More

SamplesWe examined 1,022 avian species (~ 10% recorded species) in this study, representing 418 genera, from 91 families (37% recorded families) and 29 orders (73% of all orders). Specimens were from the skin collection of the World Museum Liverpool, Tring Natural History Museum, Manchester Museum and Wollaton Hall Museum, all situated in the United Kingdom. All work was carried out in accordance with ethical regulations at Manchester Metropolitan University and with the permission of all aforementioned museums. Only the best-preserved adult specimens (no signs of cut off feathers or holes in the skin near the beak) were chosen for this study to ensure accurate measurements of bristle length, shape and presence, which should not be affected by the process of skin removal and specimen conservation. Species were randomly chosen, without targeting our sampling towards species known a priori to have bristles. Where possible, two specimens per species were measured (occurring in 82% of all species examined). Specimens of each sex were measured when present; however, this was not always possible since labelling was often inaccurate or missing. In total, the sample included 508 males, 412 females and 374 individuals of unknown sex. Both sexes were examined in 274 species and there was no difference whatsoever between the presence of bristles on male or female species (n = 97 with bristles present and n = 180 with bristles absent for both males and females). Length (Mann–Whitney U test, W = 37,962, N = 552, P = 0.94) and shape (Chi-square test, χ2 = 0, N = 552, df = 3, P = 1) of rictal bristles also did not significantly differ between males and females. Therefore, rictal bristles are likely to be sexually monomorphic and data for males and females was pooled for further analyses. Overall, rictal bristles were absent in 64% of species examined (n = 656) and just over a third of species (n = 366) had bristles present.Bristle descriptionsFacial bristles were initially identified by sight and touch in each specimen. Bristles were recorded as either present or absent from the upper rictal, lorial, lower rictal, narial and interramal regions (Fig. 1a). We use the term ‘rictal bristle’ here for bristles on both the upper rictal and/or the lorial region, since there was no clear differentiation and morphological differences between the bristles found in these regions forming a continuum of bristles above the edge of the beak. When present, rictal bristle shape was recorded as: (i) unbranched rictal bristles, (ii) rictal bristles with barbs only at the base (“Base”) and (iii) branched rictal bristles (“Branched”), i.e. barbs and barbules present along the bristle rachis (Fig. 1b). The three longest rictal bristles were measured on both sides of the head of each specimen using digital callipers, and these lengths were averaged to provide a mean length of rictal bristles per species. In species lacking rictal bristles, a length of “0” and a shape category of “Absent” was recorded.Ancestral reconstruction of facial bristle presenceFollowing Felice et al.19, a single consensus phylogenetic tree was generated from the Hackett posterior distribution of trees from Birdtree.org20 with a sample size of 10,000 post burn-in, using the TreeAnnotator utility in BEAST software21 with a burn-in of 0. Maximum Clade Credibility (MCC) with the option “-heights ca” was selected as the method of reconstruction. The common ancestor trees option (-heights ca) builds a consensus tree by summarising clade ages across all posterior trees. Both the consensus tree and posterior distribution of 10,000 trees were imported into RStudio v. 1.2.5 for R22,23 and pruned so that only species present in the dataset of this study remained in the phylogeny. Taxon names were modified where necessary to match those from the Birdtree.org (http://birdtree.org) species record. Negative terminal branches in our consensus tree were slightly lengthened to be positive using ‘edge.length[tree$edge.length  More