Ecology

Download PDF

Mya Breitbart has hunted novel viruses in African termite mounds, Antarctic seals and water from the Red Sea. But to hit pay dirt, she has only to step into her back garden in Florida. Hanging around her swimming pool are spiny-backed orbweavers (Gasteracantha cancriformis) — striking spiders with bulbous white bodies, black speckles and six scarlet spikes that make them look like a piece of medieval weaponry. Even more striking for Breitbart, a viral ecologist at the University of South Florida in St Petersburg, was what was inside. When she and her colleagues collected a few spiders and ground them up, they found two viruses previously unknown to science1.Although we humans have been focused on one particularly nasty virus since early 2020, there are legions of other viruses out there waiting to be discovered. Scientists estimate that there are about 1031 individual viral particles inhabiting the oceans alone at any given time — 10 billion times the estimated number of stars in the known Universe.It’s becoming clear that ecosystems and organisms rely on viruses. Tiny but mighty, they have fuelled evolution for millions of years by shuttling genes between hosts. In the oceans, they slice open microorganisms, spilling their contents into the sea and flooding the food web with nutrients. “Without viruses,” says Curtis Suttle, a virologist at the University of British Columbia in Vancouver, Canada, “we would not be alive.”There are just 9,110 named species listed by the International Committee on Taxonomy of Viruses (ICTV), but that’s obviously a pitiful fraction of the total. In part, that’s because officially classifying a virus used to require scientists to culture a virus in its host or host cells — a time-consuming if not impossible process. It’s also because the search has been biased towards viruses that cause diseases in humans or organisms we care about, such as farm animals and crop plants. Yet, as the COVID-19 pandemic has reminded us, it’s important to understand viruses that might jump from one host to another, threatening us, our animals or our crops.
The new scope of virus taxonomy: partitioning the virosphere into 15 hierarchical ranks
Over the past ten years, the number of known and named viruses has exploded, owing to advances in the technology for finding them, plus a recent change to the rules for identifying new species, to allow naming without having to culture virus and host. One of the most influential techniques is metagenomics, which allows researchers to sample the genomes in an environment without having to culture individual viruses. Newer technologies, such as single-virus sequencing, are adding even more viruses to the list, including some that are surprisingly common yet remained hidden until now. It’s an exciting time to be doing this kind of research, says Breitbart. “I think, in many ways, now is the time of the virome.”In 2020 alone, the ICTV added 1,044 species to its official list, and thousands more await description and naming. This proliferation of genomes prompted virologists to rethink the way they classify viruses and helped to clarify their evolution. There is strong evidence that viruses emerged multiple times, rather than sprouting from a single origin.Even so, the true range of the viral world remains mostly uncharted, says Jens Kuhn, a virologist at the US National Institute of Allergy and Infectious Diseases facility at Fort Detrick, Maryland. “We really have absolutely no idea what’s out there.”Here, there and everywhereAll viruses have two things in common: each encases its genome in a protein-based shell, and each relies on its host — be it a person, spider or plant — to reproduce itself. But beyond that general pattern lie endless variations.There are minuscule circoviruses with only two or three genes, and massive mimiviruses that are bigger than some bacteria and carry hundreds of genes. There are lunar-lander-looking phage that infect bacteria and, of course, the killer spiky balls the world is now painfully familiar with. There are viruses that store their genes as DNA, and others that use RNA; there’s even a phage that uses an alternative genetic alphabet, replacing the chemical base A in the standard ACGT system with a different molecule, designated Z.

Studies of the spiny-backed orbweaver found two viruses previously unknown to science.Credit: Scott Leslie/Minden Pictures/Alamy

Viruses are so ubiquitous that they can turn up even when scientists aren’t looking for them. Frederik Schulz did not intend to study viruses as he pored over genome sequences from waste water. As a graduate student at the University of Vienna, in 2015 he was using metagenomics to hunt for bacteria. This involves isolating DNA from a whole mix of organisms, chopping it into bits and sequencing all of them. A computer program then assembles the bits into individual genomes; it’s like solving hundreds of jigsaw puzzles whose pieces have been jumbled up.Among the bacterial genomes, Schulz couldn’t help but notice a whopper of a virus genome — obvious because it carried genes for a viral shell — with a remarkable 1.57 million base pairs2. It turned out to be a giant virus, part of a group whose members are large in terms of both genome size and absolute size (typically, 200 nanometres or more across). These viruses infect amoebae, algae and other protists, putting them in a position to influence ecosystems both aquatic and terrestrial.
Profile of a killer: the complex biology powering the coronavirus pandemic
Schulz, now a microbiologist at the US Department of Energy Joint Genome Institute in Berkeley, California, decided to search for related viruses in metagenome data sets. In 2020, in a single paper3, he and his colleagues described more than 2,000 genomes from the group that contains giant viruses; before that, just 205 such genomes had been deposited in public databases.Virologists have also looked inwards to find new species. Viral bioinformatician Luis Camarillo-Guerrero worked with colleagues at the Wellcome Sanger Institute in Hinxton, UK, to analyse metagenomes from the human gut, and built a database containing more than 140,000 kinds of phage. More than half of these were new to science. Their study4, published in February, matched others’ findings that one of the most common viruses to infect the bacteria in our guts is a group known as crAssphage (named after the cross-assembly software that picked it up in 2014). Despite its abundance, not much is known about how it contributes to our microbiome, says Camarillo-Guerrero, who now works at DNA-sequencing company Illumina in Cambridge, UK.Metagenomics has turned up a wealth of viruses, but it ignores many, too. RNA viruses aren’t sequenced in typical metagenomes, so microbiologist Colin Hill at University College Cork, Ireland, and his colleagues looked for them in databases of RNAs, called metatranscriptomes. Scientists normally use these data to understand the genes in a population that are actively being turned into messenger RNA in to make proteins, but RNA virus genomes can show up, too. Using computational techniques to pull sequences out of the data, the team found 1,015 viral genomes in metatrancriptomes from sludge and water samples5. Again, they’d massively increased the number of known viruses with a single paper.

The giant tupanvirus,found in amoebae, is more than 1,000 nanometres long and has the largest set of protein-coding genes of any known virus.Credit: J. Abrahão et al./Nature Commun.

Although it’s possible for these techniques to accidentally assemble genomes that aren’t real, researchers have quality-control techniques to guard against this. But there are other blind spots. For instance, viral species whose members are very diverse are fiendishly difficult to find because it’s hard for computer programs to piece together the disparate sequences.The alternative is to sequence viral genomes one at a time, as microbiologist Manuel Martinez-Garcia does at the University of Alicante, Spain. He decided to try trickling seawater through a sorting machine to isolate single viruses, amplified their DNA, and got down to sequencing.On his first attempt, he found 44 genomes. One turned out to represent some of the most abundant viruses in the ocean6. This virus is so diverse — its genetic jigsaw pieces so varied from one virus particle to the next — that its genome had never popped up in metagenomics studies. The team calls it 37-F6, for its location on the original laboratory dish, but Martinez-Garcia jokes that, given its ability to hide in plain sight, it should have been named 007, after fictional superspy James Bond.Virus family treesThe James Bond of ocean viruses lacks an official Latin species name, and so do most of the thousands of viral genomes discovered by metagenomics over the past decade. Those sequences presented the ICTV with a dilemma: is a genome enough to name a virus? Until 2016, proposing a new virus or taxonomic group to the ICTV required scientists to have that virus and its host in culture, with rare exceptions. But that year, after a contentious but cordial debate, virologists agreed that a genome was sufficient7.Proposals for new viruses and groups poured in (see ‘Adding to the family’). But the evolutionary relationships between these viruses were often unclear. Virologists usually categorize viruses on the basis of their shapes (long and thin, say, or a head with a tail) or their genomes (DNA or RNA, single- or double-stranded), but this says surprisingly little about shared ancestry. For example, viruses with double-stranded DNA genomes seem to have arisen on at least four separate occasions.

Source: ICTV

The original ICTV viral classification, which is entirely separate from the tree of cellular life, included only the lower rungs of the evolutionary hierarchy, from species and genus up to the order level — a tier equivalent to primates or trees with cones in the classification of multicellular life. There were no higher levels. And many viral families floated alone, with no links to other kinds of virus. So in 2018, the ICTV added higher-order levels: classes, phyla and kingdoms8.At the very top, it invented ‘realms’, intended as counterparts to the ‘domains’ of cellular life — Bacteria, Archaea and Eukaryota — but using a different word to differentiate between the two trees. (Several years ago, some scientists suggested that certain viruses might fit into the cell-based evolutionary tree, but that idea has not gained widespread favour.)The ICTV outlined the branches of the tree, and grouped RNA-based viruses into a realm called Riboviria. SARS-CoV-2 and other coronaviruses, which have single-stranded RNA genomes, are part of this realm. But then it was up to the broader community of virologists to propose further taxonomic groups. As it happened, Eugene Koonin, an evolutionary biologist at the National Center for Biotechnology Information in Bethesda, Maryland, had assembled a team to analyse all the viral genomes, as well as the latest research on viral proteins, to create a first-draft taxonomy9.They reorganized Riboviria and proposed three more realms (see ‘Virus realms’). There was some quibbling over the details, Koonin says, but the taxonomy was ratified without much trouble by ICTV members in 2020. Two further realms got the green light in 2021, but the original four realms will probably remain the largest, he says. Eventually, Koonin speculates, the realms might number up to 25.

Source: ICTV (talk.ictvonline.org/taxonomy); ICTV Coronaviridae Study Group. Nature Microbiol. 5, 536–544 (2020)

That number supports many scientists’ suspicion that there’s no one common ancestor for virus-kind. “There is no single root for all viruses,” says Koonin. “It simply does not exist.” That means that viruses probably arose several times in the history of life on Earth — and there’s no reason to think such emergence can’t happen again. “The de novo origin of new viruses, it’s still ongoing,” says Mart Krupovic, a virologist at the Pasteur Institute in Paris who was involved in both the ICTV decisions and Koonin’s taxonomy team.As to how the realms arose, virologists have several ideas. Perhaps they descended from independent genetic elements at the dawn of life on Earth, before cells even took shape. Maybe they escaped or ‘devolved’ from whole cells, ditching most of the cellular machinery for a minimal lifestyle. Koonin and Krupovic favour a hybrid hypothesis in which those primordial genetic elements stole genes from cellular life to build their virus particles. Because there are multiple origins for viruses, it’s possible there are multiple ways they’ve originated, says Kuhn, who also served on the ICTV committee and worked on the new taxonomy proposal.Thus, although the viral and cellular trees of life are distinct, the branches touch, and genes pass between the two. Whether viruses count as being ‘alive’ depends on your personal definition of life. Many researchers do not consider them to be living things, but others disagree. “I tend to believe that they are living,” says Hiroyuki Ogata, a bioinformatician working on viruses at Kyoto University in Japan. “They are evolving, they have genetic material composed of DNA and RNA, and they are very important in the evolution of all life.”The current classification is widely recognized as just the first attempt, and some virologists say it’s a bit of a mess. A score of families still lack links to any realm. “The good point is, we are trying to put some order in that mess,” says Martinez-Garcia.World changers With the total mass of viruses on Earth equivalent to that of 75 million blue whales, scientists are certain they make a difference to food webs, ecosystems and even the planet’s atmosphere. The accelerating discovery of new viruses “has revealed a watershed of new ways viruses directly impact ecosystems”, says Matthew Sullivan, an environmental virologist at Ohio State University in Columbus. But scientists are still struggling to quantify how much of an impact they have.“We don’t have a very simple story around here at the moment,” says Ogata. In the ocean, viruses can burst out of their microbial hosts, releasing carbon to be recycled by others that eat the host’s innards and then produce carbon dioxide. But, more recently, scientists have also come to appreciate that popped cells often clump together and sink to the bottom of the ocean, sequestering carbon away from the atmosphere.

Viral genomes collected from thawing permafrost at Stordalen Mire in Sweden have genes that could help break down and release carbon.Credit: Bob Gibbons/Alamy

On land, thawing permafrost is a major source of carbon, says Sullivan, and viruses seem to be instrumental in carbon release from microbes in that environment. In 2018, he and his colleagues described 1,907 viral genomes and fragments collected from thawing permafrost in Sweden, including genes for proteins that might influence how carbon compounds break down and, potentially, become greenhouse gases10.Viruses can also influence other organisms by stirring up their genomes. For example, when viruses transfer antibiotic-resistance genes from one bacterium to another, drug-resistant strains can take over. Over time, this kind of transfer can create major evolutionary shifts in a population, says Camarillo-Guerrero. And not just in bacteria — an estimated 8% of human DNA is of viral origin. For example, our mammalian ancestors acquired a gene essential for placental development from a virus.For many questions about viral lifestyles, scientists will need more than just genomes. They will need to find the virus’s hosts. A virus itself might carry clues: it could be toting a recognizable bit of host genetic material in its own genome, for example.Martinez-Garcia and his colleagues used single-cell genomics to identify the microbes that contained the newly discovered 37-F6 virus. The host, too, is one of the most abundant and diverse organisms in the sea, a bacterium known as Pelagibacter11. In some waters, Pelagibacter makes up half the cells present. If just this one type of virus were to suddenly disappear, says Martinez-Garcia, ocean life would be thrown wildly off balance.To understand a virus’s full impact, scientists need to work out how it changes its host, says Alexandra Worden, an evolutionary ecologist at the GEOMAR Helmholtz Centre for Ocean Research in Kiel, Germany. She’s studying giant viruses that carry genes for light-harvesting proteins called rhodopsins. Theoretically, these genes could be useful to the hosts — for purposes such as energy transfer or signalling — but the sequences can’t confirm that. To find out what’s going on with these rhodopsin genes, Worden plans to culture the host and virus together, and study how the pair function in the combined, ‘virocell’ state. “Cell biology is the only way you can say what that true role is, how does this really affect the carbon cycle,” she says.Back in Florida, Breitbart hasn’t cultured her spider viruses, but she’s learnt some more about them. The pair of viruses belong to a category Breitbart calls mind-boggling for their tiny, circular genomes, encoding just one gene for their protein coat and one for their replication protein. One of the viruses is found only in the spider’s body, never its legs, so she thinks it’s actually infecting some creature the spider eats. The other virus is found throughout the spider’s body, and in its eggs and spiderlings, so she thinks it’s transmitted from parent to offspring12. It doesn’t seem to be doing them any harm, as far as Breitbart can tell.With viruses, “finding them’s actually the easy part”, she says. Picking apart how viruses influence host life cycles and ecology is much trickier. But first, virologists must answer one of the toughest questions of all, Breitbart says: “How do you pick which one to study?”

Nature 595, 22-25 (2021)
doi: https://doi.org/10.1038/d41586-021-01749-7

References1.Rosario, K. et al. PeerJ 6, e5761 (2018).PubMed 
Article 

Google Scholar 
2.Schulz, F. et al. Science 356, 85–85 (2017).PubMed 
Article 

Google Scholar 
3.Schulz, F. et al. Nature 578, 432–436 (2020).PubMed 
Article 

Google Scholar 
4.Camarillo-Guerrero, L. F., Almeida, A., Rangel-Pineros, G., Finn, R. D. & Lawley, T. D. Cell 184, 1098–1109.e9 (2021).PubMed 
Article 

Google Scholar 
5.Callanan, J. et al. Sci. Adv. 6, eaay591 (2020).Article 

Google Scholar 
6.Martinez-Hernandez, F. et al. Nature Commun. 8, 15892 (2017).PubMed 
Article 

Google Scholar 
7.Simmonds, P. et al. Nature Rev. Microbiol. 15, 161–168 (2017).PubMed 
Article 

Google Scholar 
8.International Committee on Taxonomy of Viruses Executive Committee. Nature Microbiol. 5, 668–674 (2020).PubMed 
Article 

Google Scholar 
9.Koonin, E. V. et al. Microbiol. Mol. Biol. Rev. 84, e00061-19 (2020).PubMed 
Article 

Google Scholar 
10.Emerson, J. B. et al. Nature Microbiol. 3, 870–770 (2018).PubMed 
Article 

Google Scholar 
11.Martinez-Hernandez, F. et al. ISME J. 13, 232–236 (2019).PubMed 
Article 

Google Scholar 
12.Rosario, K., Mettel, K. A., Greco, A. M. & Breitbart, M. J. Gen. Virol. 100, 1253–1265 (2019).PubMed 
Article 

Google Scholar 
Download references

Related Articles

Profile of a killer: the complex biology powering the coronavirus pandemic

Should virus-naming rules change during a pandemic? The question divides virologists

The new scope of virus taxonomy: partitioning the virosphere into 15 hierarchical ranks

Subjects

Virology

Evolution

Ecology

Latest on:

Virology

Age-related immune response heterogeneity to SARS-CoV-2 vaccine BNT162b2
Article 30 JUN 21

Untangling introductions and persistence in COVID-19 resurgence in Europe
Article 30 JUN 21

SARS-CoV-2 mRNA vaccines induce persistent human germinal centre responses
Article 28 JUN 21

Evolution

Untangling introductions and persistence in COVID-19 resurgence in Europe
Article 30 JUN 21

From the archive
News & Views 29 JUN 21

Mysterious skull fossils expand human family tree — but questions remain
News 25 JUN 21

Ecology

China’s economic approach to protecting its ecology
Spotlight 29 JUN 21

Red light, green light: both signal ‘go’ to deadly algae
Research Highlight 24 JUN 21

Limited potential for bird migration to disperse plants to cooler latitudes
Article 23 JUN 21

Jobs from Nature Careers

All jobs

Advanced Research Assistant
Wellcome Trust Sanger Institute
Cambridge, United Kingdom

JOB POST

Assistant or Associate Professor of Molecular Pharmacology and Therapeutics
Columbia University Medical Center (CUMC), CU
New York, United States

JOB POST

154308 Tenure Track Assistant Professor or Associate Professor at the Department of Biomedical Sciences
University of Copenhagen (UCPH)
Copenhagen, Denmark

JOB POST

Research Associate in Synthetic Genomics
The University of Manchester
Manchester, United Kingdom

JOB POST

Nature Briefing
An essential round-up of science news, opinion and analysis, delivered to your inbox every weekday.

Email address

Yes! Sign me up to receive the daily Nature Briefing email. I agree my information will be processed in accordance with the Nature and Springer Nature Limited Privacy Policy.

Sign up More

Angeloni F, Ouborg NJ, Leimu R (2011) Meta-analysis on the association of population size and life history with inbreeding depression in plants. Biol Conserv 144:35–43Article 

Google Scholar 
Armbruster P, Reed DH (2005) Inbreeding depression in benign and stressful environments. Heredity 95:235–242CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
Barrett SCH (2010) Understanding plant reproductive diversity. Philos Trans R Soc B 365:99–109Article 

Google Scholar 
Bierne N, Tsitrone A, David P (2000) An inbreeding model of associative overdominance during a population bottleneck. Genetics 155:1981–1990CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Bijlsma R, Bundgaard J, van Putten WF (1999) Environmental dependence of inbreeding depression and purging in Drosophila melanogaster. J Evol Biol 12:1125–1137Article 

Google Scholar 
Brown KE, Kelly JK (2019) Severe inbreeding depression is predicted by the “rare allele load” in Mimulus guttatus. Evolution 74:587–596PubMed 
Article 
CAS 
PubMed Central 

Google Scholar 
Byers DL, Waller DM (1999) Do plant populations purge their genetic load? Effects of population size and mating history on inbreeding depression. Annu Rev Ecol Syst 30:479–513Article 

Google Scholar 
Campbell SA, Halitschke R, Thaler JS, Kessler A (2014) Plant mating systems affect adaptive plasticity in response to herbivory. Plant J 78:481–490CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
Carleial S, van Kleunen M, Stift M (2017) Relatively weak inbreeding depression in selfing but also in outcrossing populations of North American Arabidopsis lyrata. J Evol Biol 30:1994–2004CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
Charlesworth B (2018) Mutational load, inbreeding depression and heterosis in subdivided populations. Mol Ecol 27:4991–5003PubMed 
Article 

Google Scholar 
Charlesworth D, Willis JH (2009) The genetics of inbreeding depression. Nat Rev Genet 10:783–796CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Cheptou PO, Imbert E, Lepart J, Escarre J (2000) Effects of competition on lifetime estimates of inbreeding depression in the outcrossing plant Crepis sancta (Asteraceae). J Evol Biol 13:522–531Article 

Google Scholar 
Cheptou P-O, Donohue K (2011) Environment-dependent inbreeding depression. Its ecological and evolutionary significance. N. Phytol 189:395–407Article 

Google Scholar 
Cheptou P-O, Lepart J, Escarré J (2001) Inbreeding depression under intraspecific competition in a highly outcrossing population of Crepis sancta (Asteraceae). Evidence for frequency‐dependent variation. Am J Bot 88:1424–1429CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
Colmer TD, Voesenek LACJ (2009) Flooding tolerance. Suites of plant traits in variable environments. Funct Plant Biol 36:665CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
Crow JF (1958) Some possibilities for measuring selection intensities in man. Hum Biol 30:1–13CAS 
PubMed 
PubMed Central 

Google Scholar 
Darwin C (1878) The effects of cross and self fertilisation in the vegetable kingdom, 2nd edn. John Murray, London, UKBook 

Google Scholar 
Diekmann M (2003) Species indicator values as an important tool in applied plant ecology—a review. Basic Appl Ecol 4:493–506Article 

Google Scholar 
Dudash MR (1990) Relative fitness of selfed and outcrossed progeny in a self-compatible, protandrous species, Sabatia angularis L. (Gentianaceae). A comparison in three environments. Evolution 44:1129–1139PubMed 
Article 
PubMed Central 

Google Scholar 
Eckert CG, Barrett SCH (1994) Inbreeding depression in partially self-fertilizing Decodon verticillatus (Lythraceae). Population-genetic and experimental analyses. Evolution 48:952–964PubMed 
Article 
PubMed Central 

Google Scholar 
Ellenberg H, Weber HE, Dull R, Wirth V, Werner W, Paulissen D (2001) Zeigerwerte von Pflanzen in Mitteleuropa, 3rd edn. Erich Goltze, Göttingen
Google Scholar 
Fox CW, Reed DH (2011) Inbreeding depression increases with environmental stress. An experimental study and meta-analysis. Evolution 65:246–258PubMed 
Article 
PubMed Central 

Google Scholar 
Glémin S (2003) How are deleterious mutations purged? Drift versus nonrandom mating. Evolution 57:2678–2687PubMed 
Article 
PubMed Central 

Google Scholar 
González AM, Ron AM, de, Lores M, Santalla M (2014) Effect of the inbreeding depression in progeny fitness of runner bean (Phaseolus coccineus L.) and it is implications for breeding. Euphytica 200:413–428Article 

Google Scholar 
Hedrick PW, Kalinowski ST (2000) Inbreeding depression in conservation biology. Annu Rev Ecol Syst 31:139–162Article 

Google Scholar 
Hilbig W (1975) Übersicht über die Pflanzengesellschaften des südlichen Teiles der DDR XII – Die Großseggenrieder. Hercynia 12:341–356
Google Scholar 
Husband BC, Schemske DW (1996) Evolution of the magnitude and timing of inbreeding depression in plants. Evolution 50:54–70PubMed 
Article 
PubMed Central 

Google Scholar 
Ivey CT, Carr DE (2005) Effects of herbivory and inbreeding on the pollinators and mating system of Mimulus guttatus (Phrymaceae). Am J Bot 92:1641–1649PubMed 
Article 
PubMed Central 

Google Scholar 
Keller LF, Grant PR, Grant BR, Petren K (2002) Environmental conditions affect the magnitude of inbreeding depression in survival of Darwin’s finches. Evolution 56:1229–1239PubMed 
Article 
PubMed Central 

Google Scholar 
Keller LF, Waller DM (2002) Inbreeding effects in wild populations. Trends Ecol Evol 17:230–241Article 

Google Scholar 
Kéry M, Matthies D, Spillmann H-H (2000) Reduced fecundity and offspring performance in small populations of the declining grassland plants Primula veris and Gentiana lutea. J Ecol 88:17–30Article 

Google Scholar 
Kondrashov AS, Houle D (1997) Genotype—environment interactions and the estimation of the genomic mutation rate in Drosophila melanogaster. Proc R Soc Lond B Biol Sci 258:221–227
Google Scholar 
Li Y, van Kleunen M, Stift M (2019) Sibling competition does not magnify inbreeding depression in North American Arabidopsis lyrata. Heredity 123:723–732PubMed 
PubMed Central 
Article 

Google Scholar 
Lynch M, Walsh B (1998) Genetics and analysis of quantitative traits. Sinauer Assoc., Sunderland, MA,
Google Scholar 
O’Grady JJ, Brook BW, Reed DH, Ballou JD, Tonkyn DW, Frankham R (2006) Realistic levels of inbreeding depression strongly affect extinction risk in wild populations. Biol Conserv 133:42–51Article 

Google Scholar 
Oberdorfer E, Schwabe A (2001) Pflanzensoziologische Exkursionsflora. Für Deutschand und angrenzende Gebiete, 8th edn. E. Ulmer, Stuttgart
Google Scholar 
Reed DH, Fox CW, Enders LS, Kristensen TN (2012) Inbreeding-stress interactions. Evolutionary and conservation consequences. Ann N. Y Acad Sci 1256:33–48PubMed 
Article 
PubMed Central 

Google Scholar 
Rehling F, Matthies D, Sandner TM (2019) Responses of a legume to inbreeding and the intensity of novel and familiar stresses. Ecol Evol 9:1255–1267PubMed 
PubMed Central 
Article 

Google Scholar 
Rosche C, Hensen I, Lachmuth S (2018) Local pre-adaptation to disturbance and inbreeding-environment interactions affect colonisation abilities of diploid and tetraploid Centaurea stoebe. Plant Biol 20:75–84CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
Sandner TM, Matthies D (2016) The effects of stress intensity and stress type on inbreeding depression in Silene vulgaris. Evolution 70:1225–1238PubMed 
Article 
PubMed Central 

Google Scholar 
Sandner TM, Matthies D (2017) Interactions of inbreeding and stress by poor host quality in a root hemiparasite. Ann Bot 119:143–150PubMed 
Article 
PubMed Central 

Google Scholar 
Sandner TM, Matthies D (2018) Inbreeding limits responses to environmental stress in Silene vulgaris. Environ Exp Bot 147:86–94Article 

Google Scholar 
Schmitt J, Eccleston J, Ehrhardt DW (1987) Dominance and suppression, size-dependent growth and self-thinning in a natural Impatiens capensis population. J Ecol 75:651Article 

Google Scholar 
Schmitt J, Ehrhardt DW (1990) Enhancement of inbreeding depression by dominance and suppression in Impatiens capensis. Evolution 44:269–278PubMed 
Article 
PubMed Central 

Google Scholar 
Schrieber K, Schweiger R, Kröner L, Müller C, Campbell H (2019) Inbreeding diminishes herbivore-induced metabolic responses in native and invasive plant populations. J Ecol 107:923–936Article 

Google Scholar 
Truscott A-M, Soulsby C, Palmer SCF, Newell L, Hulme PE (2006) The dispersal characteristics of the invasive plant Mimulus guttatus and the ecological significance of increased occurrence of high-flow events. J Ecol 94:1080–1091Article 

Google Scholar 
Walisch TJ, Colling G, Poncelet M, Matthies D (2012) Effects of inbreeding and interpopulation crosses on performance and plasticity of two generations of offspring of a declining grassland plant. Am J Bot 99:1300–1313PubMed 
Article 
PubMed Central 

Google Scholar 
Waller DM (1985) The genesis of size hierarchies in seedling populations of Impatiens capensis Meerb. N. Phytol 100:243–260Article 

Google Scholar 
Waller DM, Dole J, Bersch AJ (2008) Effects of stress and phenotypic variation on inbreeding depression in Brassica rapa. Evolution 62:917–931PubMed 
Article 
PubMed Central 

Google Scholar 
Waples RS (2020) An estimator of the Opportunity for Selection that is independent of mean fitness. Evolution 74:1942–1953PubMed 
Article 
PubMed Central 

Google Scholar 
Weigel KA (2001) Controlling inbreeding in modern breeding programs. J Dairy Sci 84:E177–E184CAS 
Article 

Google Scholar 
Weiner J (1985) Size hierarchies in experimental populations of annual plants. Ecology 66:743–752Article 

Google Scholar 
Willi Y, Dietrich S, van Kleunen M, Fischer M (2007) Inter-specific competitive stress does not affect the magnitude of inbreeding depression. Evol Ecol Res 9:959–974
Google Scholar 
Willis JH (1993) Partial self-fertilization and inbreeding depression in two populations of Mimulus guttatus. Heredity 71:145–154Article 

Google Scholar 
Wright S (1922) Coefficients of inbreeding and relationship. Am Nat 56:330–338Article 

Google Scholar 
Yun L, Agrawal AF (2014) Variation in the strength of inbreeding depression across environments. Effects of stress and density dependence. Evolution 68:3599–3606PubMed 
Article 

Google Scholar  More

Genome skimming provides robust phylogenyPioneering molecular phylogenetic studies in Sepiolida that used short regions of a few mitochondrial and nuclear genes failed to resolve the relationship of major clades9,22,23. To increase the number of phylogenetically informative sites, Sanchez et al.11 sequenced and analyzed the transcriptomes of multiple species of Euprymna Steenstrup, 1887, related bobtail squids including Sepiola parva Sasaki, 1913 and Sepiola birostrata Sasaki, 1918, and several bottletail squids. They found that S. parva grouped with the Euprymna species to the exclusion of S. birostrata, and further morphological analysis led to the formal redefinition of the genus Euprymna and the reassignment of S. parva Sasaki, 1913 to Euprymna parva11. The following year, in an exhaustive study of hectocotylus structure, Bello19 proposed that Euprymna be split back into the original Euprymna Steenstrup, 1998 and a newly defined genus, Eumandya Bello 2020 that contains E. parva Sasaki, 1913 and E. pardalota Reid 2011, two taxa whose arms have two rows of suckers rather than four as in other Euprymna species. Similarly, Bello introduced a new genus, Lusepiola Bello, 2020 that has the effect of renaming Sepiola birostrata Sasaki, 1918 as Lusepiola birostrata. For clarity, we adopt the finer-grained nomenclature of Bello below, but happily note that E. parva and E. pardalota have the same abbreviations in both the notation of Sanchez et al.11 and Bello19.Sanchez et al.11 also emphasized the need for more taxon sampling, careful species assignment, and the inclusion of more informative sites when studying this group of cephalopods. However, the distribution and lifestyle of many lineages of Sepiolida makes the collection of fresh tissue for RNA sequencing very challenging. To overcome this limitation, we sequenced the genomic DNA of several Sepiolida species at shallow coverage up to 3.6× and accessed by this way several mitochondrial and nuclear loci. Most of our samples were carefully identified at the species level based on morphological characters.We recovered the mitochondrial genomes of the species targeted in this study and annotated the 13 protein-coding genes, 22 tRNAs, and two rRNAs (although only the conserved region of the large and small rRNA was obtained for Rondeletiola minor Naef, 1912).Additionally, we also downloaded the complete mitochondrial genomes of S. austrinum and Idiosepius sp., and the transcriptome of E. tasmanica available in the NCBI database. The transcriptome of E. tasmanica was used to extract its complete set of mitochondrial protein-coding genes. We could reconstruct the mitochondrial gene order for all species with complete mtDNA genomes, but we observed no re-arrangement for members of Sepiolidae, and only Sepiadarium austrinum deviated from the arrangement seen in all other Sepiadariidae (Fig. S1).To complement the mitochondrial-based evolutionary history, we also annotated several nuclear loci. As ribosomal gene clusters are present in numerous copies, they were successfully retrieved for almost all the species, except for 28S of the Sepiadariidae sp. specimen, which appeared problematic and was excluded.By mapping reads to the reference genome of E. scolopes, we obtained 3,279,410 loci shared between at least two species and further selected 5215 loci presented in most of our Sepiolidae species, but allowed some missing data in the Euprymna + Eumandya clade. This was done because the phylogenetic relationships of the Euprymna + Eumandya species were previously described in detail in Sanchez et al.11 using transcriptome data. Out of the 5215 loci, 5164 loci had a per-site coverage ranging between two and five. After trimming and removal of regions without informative sites, 577 loci remained. These ultraconserved loci had lengths ranging between 10 and 690 base pairs (bp), with an average of 65 bp. Our alignment matrix had a length of 37,512 bp and consisted of 16,495 distinct site patterns, and variable sites between 1 and 130 bp with an average value of 7 bp. We expected a low value of variable sites because these regions are highly conserved.We considered resolved nodes to be those with the ultrafast bootstrap support and posterior probability larger than 95% and 0.9, respectively. Only the very unresolved nodes were found based on the mito_nc matrix (Fig. 1). However, among the species in these nodes, Adinaefiola ligulata Naef, 1912 was well supported with amino acid sequences from mitochondrial genes (posterior probability of 1 and 94% bootstrap support) and partially by the ultraconserved loci (posterior probability of 1, but only 85% bootstrap support) as sister to the Sepiola clade (Figs. 2 and S2). Moreover, compared to the mito_nc matrix and with identical topology, mito_aa and UCEbob fully resolved the relationship of the Indo-Pacific and Mediterranean Sea Sepiolinae. The tree generated by the nuclear_rRNA produces a topology with most nodes unsupported (Fig. S3), suggesting these markers are too conserved for assessing the relationships among this group.Fig. 1: Phylogeny of Sepiolida based on nucleotide sequences from the mitochondria (mt_nc matrix).The topology of the maximum likelihood tree is shown. Numbers by the nodes indicate bootstrap support and the Bayesian posterior probabilities. Values of bootstrap support and posterior probabilities above 95% and 0.95, respectively, are not shown. (*) indicates that the node was resolved with the mito_aa and UCEbob matrices. (+) indicate that A. ligulata is sister to Sepiola using mito_aa with ultrafast bootstrap support of 94% and a posterior probability of 1. Abbreviations: IP, Indo-Pacific Ocean; MA, Mediterranean Sea, and the Atlantic Ocean.Full size imageFig. 2: Phylogenetic tree of Sepiolida based on conserved nuclear loci (UCEbob matrix).The topology of the maximum likelihood tree is shown. Numbers in by the nodes indicate the bootstrap support and the Bayesian posterior probability. Values of bootstrap support and posterior probabilities above 95% and 0.95, respectively, are not shown. IP, Indo-Pacific Ocean; MA, Mediterranean Sea, and the Atlantic Ocean.Full size imageUsing the UCEbob matrix, the topology and supported relationships of Euprymna + Eumandya species resemble those reported in Sanchez et al.11 using transcriptome sequences, proving our protocol valid when using low coverage sequencing and when a reference genome of the closest related species is available.The position of R. minor showed discordance between mitochondrial and nuclear datasets. Using the mitochondrial matrices, R. minor rendered the Sepietta Naef, 1912 clade paraphyletic, whereas using the UCEbob and rRNA_nc matrices, R. minor appeared sister to the Sepietta clade. These relationships were resolved in both mitochondrial and nuclear-based trees and require further investigation with more DNA markers and a wider population sampling.Molecular systematics of Sepiolida cladesUsing the complete mitochondrial genome, ribosomal nuclear genes, and ultraconserved loci, we recovered the monophyly of the two families of the order Sepiolida—Sepiadariidae and Sepiolidae9,24—and the monophyly of the three described subfamilies of the family Sepiolinae. However, contrary to what is proposed based on morphology in Young24, the Rossinae is not sister to all the remaining Sepiolidae but rather is sister to Heteroteuthinae, although this is unresolved in the UCE phylogeny. With the lack of systematic work on these subfamilies, our robust phylogenetic backbone in Sepiolida using new samples carefully identified by morphology and with museum vouchers, represents a notable advance to clarify the evolution of morphological traits in major clades within the family.Based on morphological characters of the hectocotylus, Bello19 recently split the polyphyletic Sepiola Leach 1817 into Lusepiola, Adinaefiola, and Boletzkyola, reserving Sepiola for the S. atlantica group sensu Naef 1923. These newly defined clades are consistent with our molecular phylogeny here and in Sanchez11, who also noted the polyphyly of Sepiola in the Indo-Pacific lineage.We find that Sepiolinae can be robustly split into two geographically distinct tribes: one that comprises species with known distribution in the Indo-Pacific region (tribe Euprymmini new tribe, defined as Sepiolinae with a closed bursa copulatrix, type genus Euprymna) and the other including all the Mediterranean and Atlantic species (tribe Sepiolini Appellof, 1989, defined here as Sepiolinae with an open bursa copulatrix, type genus Sepiola). Our molecular relationship is consistent with 13 of the 15 apomorphies used in the cladogram shown in Fig. 21 in Bello19. The other two proposed apomorphies in Bello (his apomorphic characters 4 and 6) group two IP lineages, Lusepiola and Inioteuthis, in a clade with species from the Mediterranean and Atlantic. Such relationships contradict our Euprymmini-Sepiolini sister relationship. Moreover, according to our phylogeny, apomorphy 6 of Bello, characterized by the participation of ventral and dorsal pedicels in the formation of the hectocotylus copulatory apparatus, implies that the male ancestor of Sepiolinae had a more developed hectocotylus that was simplified in the Euprymna and Eumandya clades.Among euprymins, we confirmed the monophyly of Euprymna Steenstrup 1887 as found previously by transcriptome analysis11. We also support the monophyly of Eumandya Bello, 2020 (Figs. 1 and 2), grouping the type species E. pardalota with E. parva along with the unnamed “Type 1” Ryukyuan species of Sanchez et al.11, for which only hatchlings were available. The phylogenomic grouping of Ryukyuan “Type 1” with Eumandya suggests that when its adults are found (or hatchlings are raised to maturity), its arms will carry two rows of suckers. We also found an adult of a Ryukyuan “Type 4” (extending the notation of Sanchez et al.11 in the coastal waters of Kume Island, that groups with E. scolopes from Hawaii, suggesting a divergence based on geographic isolation in the North Pacific. We also find that Lusepiola birostrata (formerly Sepiola birostrata) is grouped with Inioteuthis japonica as sister to a clade containing Euprymna, Eumandya, and an unnamed sepioline from Port Kembla, at the northeast of Martin Island in Australian waters.Among the sepiolins, we confirm the monophyly of Sepietta (only for nuclear-genome-based trees, see below). Adinaefiola, another genus erected by Bello19, with Sepiola ligulata Naef 1912 as its type species; was found sister to the Sepiola clade, but only in the tree based on amino acid mitochondrial sequences (mito_aa matrix) with a bootstrap value of 94% and a posterior probability of 1 (Fig. S2).Outside the sepiolines, members of the subfamily Heteroteuthinae are the most elusive and underrepresented in studies of cephalopod systematics due to their oceanic lifestyles. The placement of several heteroteuthin remains controversial. Lindgren et al.9, with six nuclear and four mitochondrial genes downloaded from GenBank found that Sepiolina Naef, 1912 was sister to Heteroteuthis Gray, 1849 + Rossia Owen, 1834+ Stoloteuthis Verril, 1881; rendering the subfamily Heteroteuthinae polyphyletic. In contrast, our work supports the monophyly of Heteroteuthinae by including Stoloteuthis and Heteroteuthis in this subfamily, while Rossia was placed within the Rossinae (Figs. 1, 2, S2). Members of Heteroteuthinae included in this study formed a sister group to a monophyletic Rossinae (Figs. 1, 2, S2). Semirossia, however, rendered the Rossia clade paraphyletic. Further discussion about the position of Semirossia is difficult because of the lack of information about the original source of this specimen in Kawashima et al.25.The light organ and luminescence evolutionBobtail squids are thought to use the bioluminescence of their light organ to camouflage them from predators while foraging and swimming at night through a mechanism called counter-illumination. This has been researched extensively using E. scolopes as a model system26,27,28. Unfortunately, the limited number of sequences available and the misidentification of bobtail squids in the GenBank database11,29,30 have hindered our understanding of the light organ evolution in the whole taxon.Our robust phylogeny and Bayesian reconstruction of ancestral bioluminescence clarify how the light organ and its luminescence have evolved in the family Sepiolidae. Members of Sepiolinae comprise neritic and benthic adults with bilobed light organs, except for two genera: Inioteuthis from the Indo-Pacific region, and the Sepietta species from the Mediterranean Sea and the Atlantic waters. The ancestor of the Sepiolinae very likely possessed a bilobed light organ that harbored luminescent symbiotic bacteria (Fig. 3). This character persisted until the ancestor of the euprymnins and sepiolins. Assuming that R. minor is sister to the Sepietta clade (as shown with the nuclear-based dataset, Fig. 2), it is clear that the bilobed light organ was lost once in Inioteuthis and Sepietta, and simplified to a rounded organ in R. minor. The alternative scenario, where R. minor renders the Sepietta clade paraphyletic (based on mitochondrial matrices, Fig. 1), is less plausible as it implies that the light organ was lost twice in the Sepietta group, once in S. obscura and then in the ancestor of S. neglecta and S. oweniana; or alternatively that it was lost in the ancestor of Sepietta-Rondelentiola followed by a reversion of this character in the lineage of Rondelentiola.Fig. 3: Ancestral character reconstruction (ASR).ASR of (a) the shape of the light organ and (b) the origin of luminescence in the Sepiolida. The posterior probability of each state is shown as a pie chart, mapped tree generated in BEAST (based on mito_nt matrix, see below), with the outgroups removed.Full size imageThe light organ is also present in all members of Heteroteuthinae. These bobtails are pelagic as adults, and their light organ appears as a single visceral organ rather than the bilobed form found in nektobenthic Sepiolinae. In contrast to the bacteriogenic luminescence of the light organ in E. scolopes31, previous studies in H. dispar3 failed to detect symbiotic bacteria and suggested that the luminescence has an autogenic origin. Thus, it seems plausible that the monophyly of Heteroteuthinae found in our study supports the findings in Lindgren et al.9 for convergent evolution of autogenic light organs associated with pelagic lifestyle in many squid, octopus, and Vampyroteuthis Chun, 19039,32.Divergence time of SepiolidaThe absence of fossils for this group limited our calculations of divergence time to the use of secondary calibrations. These calibrations can provide more accurate estimates depending on the type of primary calibrations that are used33. We retrieved secondary calibrations from previous estimations in Tanner et al.15, who used eleven fossil records spanning from coleoids to gastropods in transcriptome-based phylogenetic trees. Specifically, we used the time for the splits of Sepia esculenta and S. officinalis (~91 Mya), Idiosepiidae, and Sepiolida (~132 Mya) and the origin of the Decapodiformes (root age, ~174 Mya) (Fig. 4). These calibrations and our robust phylogenetic trees allow us to investigate the events that shape the divergence of some clades of the order Sepiolida (Figs. 4,  S4).Fig. 4: A chronogram of sepiolids using complete mitochondrial genes.Red dots indicate the nodes with secondary calibrations. K-Pg, refers to the Cretaceous-Paleogene boundary and MSC, to the Messinian salinity crisis.Full size imageSepiolida appeared before the Cretaceous-Paleogene extinction event34, during the middle Mesozoic around 94 Mya (95% HPD = 60.61–130.72). This time frame coincides with the rapid diversification of several oegopsida lineages15,35. Our molecular estimates also indicate that radiation of Sepiolidae and Sepidariidae occurred around the Cretaceous-Paleogene boundaries and is concurrent with the rapid diversification of modern marine percomorph fishes around the globe, after the extinction of Mesozoic fishes36,37.Among the species of Sepiolinae collected in the Mediterranean Sea for this study, only Sepiola robusta Naef, 1912, and Sepiola affinis Naef, 1912 are endemic to the Mediterranean Sea38. The distribution of the other species includes the Mediterranean Sea, North Atlantic Ocean, East Atlantic Ocean, and/or up to the Gulf of Cadiz. The confidence intervals for the split between the Mediterranean-Atlantic and Indo-Pacific lineages, and their diversification, overlap during the early Eocene to the beginning of the Oligocene (Figs. 4 and  S4). This time interval coincides with the end of the Tethys Sea, which separated the Indo-Pacific from the Mediterranean and Atlantic region through the Indian-Mediterranean Seaway39,40. This separation also influenced the divergence of loliginid clades, coinciding with the split between the Eastern Atlantic plus Mediterranean clade (Loligo, Afrololigo, Alloteuthis) and Indo-Pacific clade (Uroteuthis and Loliolus) (~55 Mya based on Fig. 2 in Anderson and Marian41).Our chronogram indicates that the ancestor of Sepiolinae arose prior to the early Eocene around 46 Mya (95% HPD = 25.16–69.49) (Fig. 4), already possessing a bilobed light organ hosting luminescent bacteria (Fig. 3). We estimate that the split between S. affinis and S. intermedia occurred around 2.62 Mya (95% HPD = 0.3–7.4) (Fig. 4) during the end of the Zanclean period, when the Atlantic Ocean refilled the Mediterranean after the Messinian salinity crisis42,43. While S. affinis is a coastal species with a narrow depth limit, S. intermedia inhabits a wider range of deeper waters. It is possible that two populations of their ancestor, each adapted to a different ecological niche and diverged sympatrically in Mediterranean waters, and, after the speciation, S. intermedia extended its distribution outside the Mediterranean to the Gulf of Cadiz44.We also estimate that the split between H. dispar Rüppell, 1844 and H. hawaiiensis (Berry, 1909) occurred around 2.4 Mya (95% HPD = 0.46–5.88), coinciding with the closure of the Isthmus of Panama around 2.8 Mya45. Surveys of these species found H. hawaiiensis in the North Pacific and H. dispar in the North Atlantic Ocean and Mediterranean Sea46. A recent speciation event might be the reason for the lack of morphological differences between the two species46. Thus, these species may be rendered as cryptic species, a phenomenon increasingly reported in oceanic cephalopods47. The sister species of this cryptic species complex, H. dagamensis Robson, 1924, appeared before, around 6 Mya, and is reported with broad distribution in the South Atlantic Ocean off South Africa, the Gulf of Mexico, North Atlantic Ocean between Ireland and Newfoundland in Canada, and the South Pacific Ocean off New Zealand48,49,50.The origin of the Heteroteuthis ancestor of H. dispar, H. hawaiiensis, and H. dagamensis can be placed in the Pacific Ocean. After the formation of the Isthmus of Panama, the northern population of Heteroteuthis might have split into H. hawaiiensis in North Pacific and H. dispar in the Atlantic Ocean (from where it also migrated to the Mediterranean Sea). Meanwhile, the formation of the equatorial currents isolated the southern population of Heteroteuthis and gave rise to H. dagamensis. Then, H. dagamensis extended its distribution from the Southern Pacific to the South Atlantic Ocean, the North Atlantic waters, and the Gulf of Mexico. Analysis of molecular species delimitation, however, suggests that H. dagamensis includes cryptic lineages among Atlantic and New Zealand populations30.While the origin of Heteroteuthis might also be in the Atlantic Ocean, the higher diversity of heteroteuthins in the Pacific (H. hawaiiensis, H. dagamensis, H. ryukyuensis Kubodera, Okutani and Kosuge, 2009, H. nordopacifica Kubodera and Okutani, 2011, and an unknown H. sp. KER (only known from molecular studies49)) than at the Atlantic (H. dispar and H. dagamensis), make its origin at the Atlantic less plausible. Moreover, the Atlantic Heteroteuthis were found nested within Heteroteuthinae species from the Pacific, supporting Pacific Ocean origin (Figs. 1, 4).By sequencing the genomic DNA of sepiolids at low coverage, we recovered complete mitochondrial genomes and nuclear ribosomal genes for most of our collections. Furthermore, mapping reads to the reference genome of E. scolopes allowed us to retrieve additional nuclear-ultraconserved regions. We demonstrate that these nuclear and mitochondrial loci are useful to reconstruct robust phylogenetic trees, especially when the transcriptomes of specimens are difficult to collect, as for sepiolids inhabiting oceanic environments. Finally, our study integrated genomic DNA sequencing with confident morphological identification, which helped to reconstruct the ancestral character of the light organ and its luminescence in sepiolids, and clarify how major lineages have evolved, establishing the existence of distinct Indo-Pacific and Mediterranean-Atlantic subfamilies of Sepiolinae. Our collections and genomically anchored phylogenies will provide a reliable foundation classification of sepiolids for future studies. More

To determine whether exposure to (Z)-11-16:Ald induced detectable defence-related responses in B. nigra, we investigated the amount of feeding damage to plants and the volatile emissions of plants exposed to volatilised pure compound. Contrary to our hypothesis, P. xylostella larvae fed more on plants previously exposed to (Z)-11-16:Ald at biologically relevant levels than on controls, suggesting that exposure to (Z)-11-16:Ald might increase the susceptibility of plants to future herbivory. Earlier work investigating the responses of plants to insect sex pheromone showed that exposure primed defences14,16, resulting in higher volatile emissions13, and lower feeding12. Our results are different to these findings and suggest alternative ecological effects of detecting insect pheromone to those reported earlier.The greater feeding on (Z)-11-16:Ald-exposed plants compared to controls could relate to the differences in volatile emissions of the differently treated plants. Although there were no detectable differences in volatile emissions between exposed plants and controls immediately after (Z)-11-16:Ald-exposure, there were differences after a subsequent 24 h of feeding. The difference could be due to (Z)-11-16:Ald-induced changes in the plant directly affecting herbivore-induced volatile emission, or could be related to altered plant nutrition or defences leading to an increase in the leaf area consumed by herbivores and a consequent difference in induction of volatile emissions. Thus (Z)-11-16:Ald could potentially alter plant defences and in doing so increase the survival of eggs and future hatched larvae. It has earlier been shown in Arabidopsis thaliana that application of Pieris brassicae or Spodoptera littoralis egg extract onto leaves reduced the induction of genes related to defence against insects after caterpillar feeding23. However, at this stage our data does not allow us to further explore these ecological hypotheses. Additional studies would need to test in greater detail the deposition of eggs, the hatching of eggs, and performance and survival of developing larvae of plants exposed to (Z)-11-16:Ald.The observation of changes in feeding amount and plant volatile responses prompted us to examine the early and late signaling events in response to (Z)-11-16:Ald exposure. The results showed that exposure to (Z)-11-16:Ald induced a transmembrane depolarisation by plants. The depolarisation of the plasma membrane potential is known to be the first detectable event in the detection of a biotic or abiotic stress19. We estimated the detection threshold of the pheromone to be a concentration between 25 and 10 ppm. A lower level of transmembrane depolarisation was observed when plants were exposed to (Z)-11-16:Ac, hence (Z)-11-16:Ald appears to be the most phytoactive component of the P. xylostella sex pheromone. However, it remains unclear how specific the plant response is to the (Z)-11-16:Ald, which should be further elucidated by comparison with aldehydes that have a more similar chemical structures.The plasma membrane is the only cellular structure in direct contact with the environment, which makes it critical for sensing environmental stimuli and initiating a cascade of events that eventually leads to a specific response21. As demonstrated in Arabidopsis, Vm depolarisation depends on a cascade of events that include changes in [Ca2+]cyt and the production of ROS24. We found that 50 ppm and 100 ppm of (Z)-11-16:Ald increased both the [Ca2+]cyt and the production of ROS. Moreover, we demonstrated that the increased ROS detected by CLSM were associated with the increased expression of reactive oxygen species (ROS)-mediated genes and ROS-scavenging enzyme activity. ROS participate in cell oxidation, during which H2O2 is produced and later regulated by ROS-scavenging enzymes involved in its degradation to protect cells from oxidative stress25,26. The production of H2O2 is potentially harmful and can result in oxidation in the cells of aerobic living organisms27. It is also an important component of the signalling network in plants28,29 and takes part in plant defence in response to environmental stress30,31. Several plant species trigger localized cell death by producing ROS at oviposition sites on leaves, which has been shown to be associated with an increase in egg mortality or a reduction of larval survival rate32. Recently, Bittner and colleagues demonstrated that when plants are previously exposed to the female pine sawfly (D. pini) sex pheromone, they produce H2O2 and induce defence-related genes faster after egg deposition on leaves, compared to plants that have not been exposed to the pheromone15. They suggested that the sex pheromone acted as an environmental cue indicating to plants that there would be future egg deposition on needles and subsequent herbivory. Plants then responded to it by producing H2O2, which formed necrotic tissue and reduced survival of eggs. Taken together with our observations, it is possible that the (Z)-11-16:Ald also primes B. nigra plant defences by producing H2O2 as a defensive mechanism to limit future egg deposition. It was shown that B. nigra plants induce the necrosis of cells located at the oviposition site of Pieris rapae and Pieris napi33, which can support this hypothesis. Contrary to the hypothesis of defence induction, B. nigra plants exposed to (Z)-11-16:Ald received more herbivore-feeding damage than control plants, but further investigations are needed to determine whether the production of H2O2 following exposure to (Z)-11-16:Ald would reduce subsequent egg deposition or hatching.Interestingly, early signalling events following the exposure to (Z)-11-16:Ald are analogous to the responses induced by biotic stress such as herbivore-wounding20,22, and in response to HIPVs21. For example, exposure to HIPVs resulted in a depolarisation of the plasma membrane (Vm)34,35 due to the entrance of calcium (Ca2+) into the cytosol of cells34. The detection of HIPVs results in transcriptional36,37, metabolic and physiological changes in plants38. Several studies have shown that plants perceiving HIPVs deploy faster and stronger chemical defences upon subsequent stress8,9, which can negatively affect herbivorous insects39. The observed action of (Z)-11-16:Ald is typical to that of insect and mite elicitors21,40,41. However, it is notable that to determine if (Z)-11-16:Ald induced detectable responses in B. nigra plants, whole plants were exposed to vapourised (Z)-11-16:Ald for 24 h, while to determine early and late signalling events following exposure to (Z)-11-16:Ald we applied 100 ppm in aqueous state for 30 min to 1 h (Table S2). While we used a biologically realistic scenario with a realistic concentration of pheromone for the whole plant responses the in vitro experiments focussed on mechanisms do not represent biologically accurate scenarios and utilised high concentrations of pheromone. Future studies should bridge this methodological gap by utilising more biologically realistic scenarios in mechanism elucidation.(Z)-11-16:Ald has been found to be a main constituent of pheromones in many moth species from the Noctuidae family including the corn earworm Helicoverpa zea42, which is the second-most important economic pest species in North America43. Many plants that have co-evolved with the Noctuidae could have the ability to detect (Z)-11-16:Ald. Prior to this study, the ability of plants to detect insect-emitted volatiles had been reported for two species: a perennial plant, S. altissima, and a tree, P. sylvestris. We can now tentatively add an annual plant, B. nigra, to the list. These studies suggest that the ability to detect insect-emitted volatiles has widely evolved in a large variety of plant families, and highlight the need to determine how widespread this trait is. Further studies should also determine the threshold and distance of detection and the ecological consequences of this detection.In summary our results indicate for the first time that exposing B. nigra plants to volatile (Z)-11-16:Ald increases the susceptibility of plants to feeding by P. xylostella larvae and induces an alteration in herbivore-induced volatile emissions. Further mechanistic experiments conducted in vitro using high doses of pheromone indicated that exposure to (Z)-11-16:Ald induces responses in receiver plants that are characterised by a depolarisation of Vm, an increase in [Ca2+]cyt and production of H2O2 leading to an increase in ROS-mediated gene expression and ROS scavenging-enzyme activity, which are typical responses to insect elicitors. This study supports recent findings showing that plants can detect insect-emitted volatiles. However, further research should be conducted to determine an accurate dose response of whole plants to volatile pheromone and the specificity of the response to this particular aldehyde.Materials and methodsThe study complied with local and national regulations.PlantsBrassica nigra (black mustard) seeds were collected from wild populations in the Netherlands (supplied by E. Poelman, Wageningen University). For experiment 1, conducted in Kuopio (Finland), plants were grown in plastic pots (8 × 8 cm) containing a mix of peat, soil and sand (3:1:1), in plant growth chambers (Weiss Bio 1300 m, Germany) with a 16L: 8D and light intensity of 250 µmol m−2 s−1. The temperature was maintained at 21 °C with 60% relative humidity during the day and decreased to 16 °C with 80% RH during the night time. Plants were watered every day and fertilized twice per week with a 0.1% solution containing nitrogen, phosphorous and potassium in a 19:4:20 ratio (Kekkilä Oyj, Finland).For experiments conducted in Turin (Italy) (Experiments 2 to 4), plants were grown in plastic pots (8 × 8 cm) containing a mix of peat, soil (Klasmann-Deilmann, Germany), sand (Vimark, Italy) and vermiculite (Unistara, Italy), in a climate-controlled room at 22 ± 1 °C, 16L:8D with light intensity of 250 µmol m−2 s−1. Three and four-week-old plants were used for all the experiments.Insects and synthetic compoundsPlutella xylostella were reared on broccoli (B. oleracea var. italica) with an artificial light–dark cycle of 16L:8D at 22 ± 0.5 °C.Synthetic (Z)-11-hexadecenal [(Z)-11-16:Ald] and (Z)-11-hexadecenyl acetate [(Z)-11-16:Ac] (isomeric purity 93%), were purchased from Pherobank (Wageningen, The Netherlands).Exposure of plants to treatmentsA 100 ppm solution of synthetic (Z)-11-16:Ald diluted in dichloromethane was prepared and 100 µl of the solution was injected into a rubber septum (7 mm O.D; Sigma-Aldrich) and left for 30 min for the dichloromethane to evaporate. As a solvent control, 100 µl of dichloromethane was deposited on a rubber septum, and left to evaporate for 30 min. A second control, without the rubber septum and dichloromethane, was also set up. Either treatment or control septa were enclosed in 0.5 L glass jars connected with Teflon tubing to plastic bags (Polyethylene terephthalate; overall dimensions 28 × 35 cm; Look Isopussi Eskimo oy, Finland). Plants were enclosed in the PET bags, that were previously baked for 1 h at 120 °C. For 24 h, a clean air flow was passed into the treatment or control glass jars and then into the bags containing the plants (Fig. S1 and Table S2).Volatile collections and feeding assaysAfter 24 h of exposure to the treatments, jars containing the rubber septa were disconnected from the plants, and a first VOC sampling was made before adding 22 first and second instar P. xylostella larvae were added to each plant for 24 h. After 24 h of feeding, volatile compounds were collected by dynamic headspace sampling (Fig. S1). VOCs were trapped by pulling clean air at 0.22 L min−1 for one hour through steel tubes filled with 200 mg Tenax TA 60/80 adsorbent (Markes International Ltd, UK) using a vacuum pump (KNF, Germany). The collected volatiles were thermally desorbed and analysed by gas chromatography-mass spectrometry (Agilent 7890A GC, and Agilent MS model 5975C VL MSD; New York, USA). Trapped compounds were desorbed with a thermal desorption unit (TD-100; Markes International Ltd, Llantrisant, UK) at 250 °C for 10 min, and cryofocused at − 10 °C in splitless mode onto an HP‐5 capillary column (50 m, 0.2 mm i.d, 0.5 μm film thickness; Hewlett‐Packard). The oven temperature programme was held at 40 °C for 1 min and then ramped 5 °C min−1 to 210 °C, and then further ramped at 20 °C min−1 to 290 °C. The carrier gas was helium, the transfer line temperature to the MSD was 300 °C, the ionization energy was 70 eV, and the full scan range was 29–355 m/z. We identified volatiles by comparison with a series of analytical standards (Sigma-Aldrich, Germany) and by comparison of their mass spectra to those in the NIST and Wiley 275 mass spectral libraries. Compound quantification was based on Total Ion Chromatograms (TIC) and according to the responses of analytical standards. Emission rates (ER) were calculated with the formula ER = X*Ai/Dw*t*Ao. ER was expressed in ng gDW−1 h−1, X was the compound quantity (ng), Ai and Ao were the inlet and outlet air flows (mL min-1), respectively, t was the sampling time of one hour and Dw is the dry weight of the plant sampled (g).After the volatile collection, we placed leaves of plants on A4 paper for scaling and digitally photographed them. Plants were then dried in paper bags in an oven at 60 °C for 3 days. The leaf area consumed by larvae was calculated using the LeafAreaAnalyzer software (https://github.com/EmilStalvinge/LeafAreaAnalyser; emilstalvinge@gmail.com).Transmembrane potential (Vm) measurementsLeaf segments (0.5 cm2) of three individual B. nigra plants were placed in Eppendorf tubes for 1 h with either 1 ml of (Z)-11-16:Ald solution, (Z)-11-16:Ac solution or control (Table S2). Vm was determined by inserting an electrode into a plant leaf segment following a method detailed earlier44. Vm variations were measured upon perfusion of four concentrations: 10, 25, 50, and 100 ppm diluted in MES buffer (pH 6.5) + 0.1% Tween 20 (V/V). A solution of MES buffer (pH 6.5) + 0.1% Tween 20 (V/V) was used as control. Vm values were recorded every five seconds.Determination of intracellular calcium variations using confocal laser scanning microscopy (CLSM) and Calcium OrangeCalcium Orange dye (stock solution in DMSO, Molecular Probes, Leiden, The Netherlands) was diluted in 5 mM MES-Na buffer (pH 6.0) to a final concentration of 5 µM. This solution was applied to B. nigra leaves attached to the plant as previously reported20,24,45. After 1 h incubation with Calcium Orange, the leaf was mounted on a Leica TCS SP2 (Leica Microsystems Srl, Milan, Italy) multiband confocal laser scanning microscope (CLSM) stage, without separating the leaf from the plant, in order to assess the basic fluorescence levels as control. A 50 µl application of either 50 or 100 ppm (Z)-11-16:Ald was made and after 30 min the calcium signature was observed. The microscope operates with a Krypton/Argon laser at 543 nm and 568 nm wavelengths: the first wavelength excites Calcium Orange, resulting in green fluorescence and the second mainly excites chlorophyll, resulting in red fluorescence. All images were obtained with an objective HCX APO 40 × immersed in water with a numeric aperture (NA) of 0.8. The scan speed was set at 400 Hz (Hz). The microscope pinhole was set at 0.064 mm and the average size depth of images was between 65 and 70 µm; the average number of sections per image was 25. The image format was 1024 × 1024 pixels, 8 bits per sample and 1 sample per pixel.CLSM Subcellular localization of H2O2 and active peroxidases using 10-acetyl-3,7-dihydroxyphenoxazine (Amplex Red)B. nigra leaves from rooted potted plants were treated with 50 µl of either 50 or 100 ppm of (Z)-11-16:Ald (Table S2) after incubation with the dye 10-acetyl-3,7-dihydroxyphenoxazine (Amplex Red) as described earlier22. The Molecular Probes Amplex Red Hydrogen Peroxide/Peroxidase Assay kit (A-22188) was used and dissolved in MES-Na buffer 50 mM (pH6.0) containing 0.5 mM calcium sulfate to obtain a 50 μM solution. Leaves where then mounted on a Leica TCS SP2 miscroscope as described above. Scannings were recorded after 180 min using the HCX PL APO 63x/1.20 W Corr/0.17CS objective. The microscope was operated with a Laser Ar (458 nm/5 mW; 476 nm/5 mW; 488 nm/20 mW; 514 nm/20 mW); a Laser HeNe 543 nm/1,2 mW and a Laser HeNe 633 nm/10 mW.ROS-scavenging enzyme activities and soluble protein determinationLeaves were collected immediately after 30 min of exposure to either 100 ppm of (Z)-11-16:Ald or control (Table S2). Intact leaves of two pooled plants were frozen in liquid N2 and stored at -80ºC before enzyme extraction. Frozen leaves were used for extraction of ROS scavenger enzymes following the method described in Maffei et al.22. All steps were carried out at 4 °C. Plant material was ground with a mortar and pestle under liquid nitrogen in cold 50 mM sodium phosphate buffer, pH 7.5, containing 250 mM sucrose, 1.0 mM EDTA, 10 mM KCl, 1 mM MgCl2, 0.5 mM phenylmethylsulfonyl fluoride (PMSF), 0.1 mM dithiothreitol (DTT), and 1% (w/v) polyvinylpolypyrrolidone (PVPP) in a 1:10 proportion (weight of plant material to buffer volume). The homogenate was then centrifuged at 25,000 g for 20 min at 4 °C and the supernatant was used directly for measurement of enzyme activity. The soluble protein concentration was measured using the method established by Bradford46 using bovine serum albumin as a standard.Catalase (CAT) activity was assayed spectrophotometrically by monitoring the absorbance change at 240 nm due to the decreased absorption of H2O2 (ɛ = 39.4 mM−1 cm−1). The reaction mixture in 1 mL final volume contained 50 mM Na-P, pH 7.0, 15 mM H2O2, and the enzyme extract. The reaction was initiated by addition of H2O2.Peroxidase (POX) activity was measured by detecting the oxidation of guaiacol (ɛ = 26.6 mM−1 cm−1) in the presence of H2O2. The reaction mixture contained 50 mM Na-P, pH 7.0, 0.33 mM guaiacol, 0.27 mM H2O2, and the enzyme extract in a 1.0 mL final volume. The reaction was started by addition of guaiacol and measured spectrophotometrically at 470 nm.Superoxide dismutase (SOD) activity was measured by reduction of nitro blue tetrazolium due to a photochemically generated superoxide anion. One ml of assay mixture consisted of 50 mM Na-P buffer, pH 7.8, 13 mM methionine, 75 μM nitro blue tetrazolium (NBT), 2 μM riboflavin, 0.1 mM EDTA, and the enzyme extract. Riboflavin was added as the last reagent. Samples were placed 30 cm below a light source (60 µmol m−1 s−1), and the reaction was allowed to run for 15 min. The reaction was stopped by switching off the light. A non irradiated reaction mixture, which was run in parallel, did not develop colour and served as a control. The absorbance was read at 560 nm.Quantitative gene expression analysis by Real-time PCRTotal RNA was isolated from control or treatment leaf tissues using RNA Isolation mini Kit (Machery-Nagel, Germany), and RNase- Free DNase according to the manufacturer’s protocols. The quality of RNA was checked in 1% agarose gel and the final yield was checked with a Spectrophotometer (Pharmacia Biotech Ultrospec 3000, United States). The cDNA synthesis was performed starting from 1 µg RNA using the High Capacity cDNA Reverse Transcription kit (Applied Biosystem, United States). Primers for real-time PCR were designed using the Primer 3.0 software47 and the relative sequences are listed in Supplementary Table S3. The real-time PCR was performed on an Mx3000P Real-Time System (Agilent Technologies, United States) using SYBR green I with the dye ROX as an internal loading standard. The reaction mixture was 10 µl in volume and comprised 5 µl of 2 × Maxima SYBR Green qPCR Master Mix (Thermo Fisher Scientific), 0.5 ml of cDNA, and 100 nM of primers (Integrated DNA Technologies, United States). The thermal conditions were as follows: 10 min at 95 °C, 40 cycles of 15 s at 95 °C, 20 s at 57 °C, and 30 s at 72 °C. Fluorescence was read after each annealing and extension phase. All runs were followed by a melting curve analysis from 55 to 95 °C. Two reference genes, ACT1 and eEF1Balpha2, were used to normalize the results. The sequences of the primers used in this work for CAT1, CuZnSOD1, PER4, ACT1 and eEF1Balpha2 are reported in Table S3. All amplification plots were analyzed with the MX3000P software (Agilent Technologies, United States) to obtain Ct values. Real-time PCR data are expressed as fold change of the treatment with respect to the control.Statistical analysisArea consumed by larvae, Vm measurements, enzyme activity and gene expression data were analysed using SPSS 25 software (IBM Corp. Armonk, USA). The normality of data and homogeneity of variances were checked and log transformed when the data did not meet assumptions for parametric analyses. Because we observed no significant differences between the control and the solvent control, we directly compared the control with (Z)-11-16:Ald for all analyses. Differences between treatments were analysed using T-tests for the fed area, gene expression and enzyme activity. Volatile emission rates were log transformed and auto-scaled (mean-centered and divided by the standard deviation of each variable). Partial Least Squares – Discriminant Analysis (PLS-DA) was performed on emission rates with the R software (v. 3.4.3) with the package vegan and RVAideMemoire with cross validation based on 50 submodels (fivefold outer loop and fourfold inner loop). A pairwise test was performed, based on PLS-DA with 999 permutations, to highlight the differences between treatments. The PLS-DA graphics were done with metaboanalyst (https://www.metaboanalyst.ca/). More

1.Weishampel, D. B., Dodson, P. & Osmólska, H. The Dinosauria 2nd edn (University of California Press, 2004).2.Fastovsky, D. E. & Weishampel, D. B. The Evolution and Extinction of the Dinosaurs (Cambridge University Press, 2005).3.Brusatte, S. L. et al. The extinction of the dinosaurs. Biol. Rev. 90, 628–642 (2015).PubMed 
Article 
PubMed Central 

Google Scholar 
4.Alvarez, L. W., Alvarez, W., Asaro, F. & Michel, H. V. Extraterrestrial cause for the Cretaceous-Tertiary extinction. Science 208, 1095–1108 (1980).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
5.Chiarenza, A. A. et al. Asteroid impact, not volcanism, caused the end-Cretaceous dinosaur extinction. Proc. Natl Acad. Sci. USA 117, 17084–17093 (2020).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
6.Schulte, P. et al. The Chicxulub asteroid impact and mass extinction at the Cretaceous-Paleogene boundary. Science 327, 1214–1218 (2010).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
7.Russell, D. A. The gradual decline of the dinosaurs—fact or fallacy? Nature 307, 360–361 (1984).ADS 
Article 

Google Scholar 
8.Sloan, R. E., Rigby, J. K., Van Valen, L. M. & Gabriel, D. Gradual dinosaur extinction and simultaneous ungulate radiation in the Hell Creek Formation. Science 232, 629–633 (1986).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
9.Sheehan, P. M., Fastovsky, D. E., Hoffmann, R. G., Berghaus, C. B. & Gabriel, D. L. Sudden extinction of the dinosaurs: Latest Cretaceous, upper Great Plains, USA. Science 254, 835–839 (1991).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
10.Sakamoto, M., Benton, M. J. & Venditti, C. Dinosaurs in decline tens of millions of years before their final extinction. Proc. Natl Acad. Sci. USA 113, 5036–5040 (2016).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
11.Chiarenza, A. A. et al. Ecological niche modelling does not support climatically-driven dinosaur diversity decline before the Cretaceous/Paleogene mass extinction. Nat. Commun. 10, 1091 (2019).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
12.Russell, L. S. Body temperature of dinosaurs and its relationships to their extinction. J. Paleontol. 39, 497–501 (1965).
Google Scholar 
13.Brusatte, S. L., Butler, R. J., Prieto-Márquez, A. & Norell, M. A. Dinosaur morphological diversity and the end-Cretaceous extinction. Nat. Commun. 3, 804 (2012).ADS 
PubMed 
Article 
CAS 
PubMed Central 

Google Scholar 
14.Benson, R. B. J. et al. Rates of dinosaur body mass evolution indicate 170 million years of sustained ecological innovation on the avian stem lineage. PLoS Biol. 12, e1001853 (2014).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
15.Rezende, E. L., Bacigalupe, L. D., Nespolo, R. F. & Bozinovic, F. Shrinking dinosaurs and the evolution of endothermy in birds. Sci. Adv. 6, eaaw4486 (2020).ADS 
PubMed 
PubMed Central 
Article 

Google Scholar 
16.Lloyd, G. T. et al. Dinosaurs and the Cretaceous Terrestrial Revolution. Proc. R. Soc. B Biol. Sci. 275, 2483–2490 (2008).Article 

Google Scholar 
17.Gates, T. A., Prieto-Márquez, A. & Zanno, L. E. Mountain building triggered Late Cretaceous North American megaherbivore dinosaur radiation. PLoS ONE 7, e42135 (2012).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
18.Loewen, M. A., Irmis, R. B., Sertich, J. J. W., Currie, P. J. & Sampson, S. D. Tyrant dinosaur evolution tracks the rise and fall of late Cretaceous oceans. PLoS ONE 8, e79420 (2013).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
19.Archibald, J. D. et al. Cretaceous extinctions: Multiple causes. Science 328, 973 (2010).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
20.Mitchell, J. S., Roopnarine, P. D. & Angielczyk, K. D. Late Cretaceous restructuring of terrestrial communities facilitated the end-Cretaceous mass extinction in North America. Proc. Natl Acad. Sci. USA 109, 18857–18861 (2012).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
21.Schoene, B. et al. U-Pb constraints on pulsed eruption of the Deccan Traps across the end-Cretaceous mass extinction. Science 363, 862–866 (2019).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
22.Sprain, C. J. et al. The eruptive tempo of Deccan volcanism in relation to the Cretaceous-Paleogene boundary. Science 363, 866–870 (2019).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
23.Hull, P. M. et al. On impact and volcanism across the Cretaceous-Paleogene boundary. Science 367, 266–272 (2020).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
24.Landman, N. H. et al. Ammonite extinction and nautilid survival at the end of the Cretaceous. Geology 42, 707–710 (2014).ADS 
CAS 
Article 

Google Scholar 
25.Longrich, N. R., Martill, D. M. & Andres, B. Late Maastrichtian pterosaurs from North Africa and mass extinction of Pterosauria at the Cretaceous-Paleogene boundary. PLoS Biol. 16, e2001663 (2018).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
26.Longrich, N. R., Tokaryk, T. & Field, D. J. Mass extinction of birds at the Cretaceous-Paleogene (K-Pg) boundary. Proc. Natl Acad. Sci. USA 108, 15253–15257 (2011).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
27.Longrich, N. R., Bhullar, B.-A. S. & Gauthier, J. A. Mass extinction of lizards and snakes at the Cretaceous-Paleogene boundary. Proc. Natl Acad. Sci. USA 109, 21396–21401 (2012).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
28.Fastovsky, D. E. et al. Shape of Mesozoic dinosaur richness. Geology 32, 877–880 (2004).ADS 
Article 

Google Scholar 
29.Archibald, J. D. in Volcanism, Impacts, and Mass Extinctions: Causes and Effects (eds. Keller, G. & Kerr, A. C.) 213–224 (The Geological Society of America Special Paper 505, 2014).30.Wang, S. C. & Dodson, P. Estimating the diversity of dinosaurs. Proc. Natl Acad. Sci. USA 103, 13601–13605 (2006).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
31.Starrfelt, J. & Liow, L. H. How many dinosaur species were there? Fossil bias and true richness estimated using a Poisson sampling model. Philos. Trans. R. Soc. B Biol. Sci. 371, 20150219 (2016).Article 
CAS 

Google Scholar 
32.Bonsor, J. A., Barrett, P. M., Raven, T. J. & Cooper, N. Dinosaur diversification rates were not in decline prior to the K-Pg boundary. R. Soc. Open Sci. 7, 201195 (2020).ADS 
PubMed 
PubMed Central 
Article 

Google Scholar 
33.Benton, M. J., Wills, M. A. & Hitchin, R. Quality of the fossil record through time. Nature 403, 534–537 (2000).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
34.Alroy, J. et al. Effects of sampling standardization on estimates of Phanerozoic marine diversification. Proc. Natl Acad. Sci. USA 98, 6261–6266 (2001).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
35.Alroy, J. et al. Phanerozoic trends in the global diversity of marine invertebrates. Science 321, 97–100 (2008).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
36.Close, R. A., Evers, S. W., Alroy, J. & Butler, R. J. How should we estimate diversity in the fossil record? Testing richness estimators using sampling-standardised discovery curves. Methods Ecol. Evol. 9, 1386–1400 (2018).Article 

Google Scholar 
37.Silvestro, D., Salamin, N., Antonelli, A. & Meyer, X. Improved estimation of macroevolutionary rates from fossil data using a Bayesian framework. Paleobiology 45, 546–570 (2019).Article 

Google Scholar 
38.Close, R. A., Benson, R. B. J., Saupe, E. E., Clapham, M. E. & Butler, R. J. The spatial structure of Phanerozoic marine animal diversity. Science 368, 420–424 (2020).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
39.Benton, M. J. Scientific methodologies in collision: The history of the study of the extinction of the dinosaurs. Evol. Biol. 24, 371–400 (1990).
Google Scholar 
40.Butler, R. J., Benson, R. B. J., Carrano, M. T., Mannion, P. D. & Upchurch, P. Sea level, dinosaur diversity and sampling biases: Investigating the ‘common cause’ hypothesis in the terrestrial realm. Proc. R. Soc. B Biol. Sci. 278, 1165–1170 (2011).Article 

Google Scholar 
41.Zaffos, A., Finnegan, S. & Peters, S. E. Plate tectonic regulation of global marine animal diversity. Proc. Natl Acad. Sci. USA 114, 5653–5658 (2017).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
42.East, M., Müller, R. D., Williams, S., Zahirovic, S. & Heine, C. Subduction history reveals Cretaceous slab superflux as a possible cause for the mid-Cretaceous plume pulse and superswell events. Gondwana Res. 79, 125–139 (2020).ADS 
Article 

Google Scholar 
43.Grasby, S. E., Them, T. R., Chen, Z., Yin, R. & Ardakani, O. H. Mercury as a proxy for volcanic emissions in the geologic record. Earth Sci. Rev. 196, 102880 (2019).CAS 
Article 

Google Scholar 
44.Miller, K. G. et al. The Phanerozoic record of global sea level change. Science 310, 1293–1298 (2005).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
45.Ray, D. C. et al. The magnitude and cause of short-term eustatic Cretaceous sea-level change: a synthesis. Earth Sci. Rev. 197, 102901 (2019).Article 

Google Scholar 
46.Coiffard, C., Gomez, B., Daviero-Gomez, V. & Dilcher, D. L. Rise to dominance of angiosperm pioneers in European Cretaceous environments. Proc. Natl Acad. Sci. USA 109, 20955–20959 (2012).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
47.Chaboureau, A.-C., Sepulchre, P., Donnadieu, Y. & Franc, A. Tectonic-driven climate change and the diversification of angiosperms. Proc. Natl Acad. Sci. USA 111, 14066–14070 (2014).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
48.Magallón, S., Gómez-Acevedo, S., Sánchez-Reyes, L. L. & Hernández-Hernández, T. A metacalibrated time-tree documents the early rise of flowering plant phylogenetic diversity. N. Phytol. 207, 437–453 (2015).Article 

Google Scholar 
49.Magallón, S., Sánchez-Reyes, L. L. & Gómez-Acevedo, S. L. Thirty clues to the exceptional diversification of flowering plants. Ann. Bot. 123, 491–503 (2019).PubMed 
Article 
PubMed Central 

Google Scholar 
50.Meredith, R. W. et al. Impacts of the Cretaceous terrestrial revolution and KPg extinction on mammal diversification. Science 334, 521–524 (2011).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
51.Grossnickle, D. M. & Newham, E. Therian mammals experience an ecomorphological radiation during the Late Cretaceous and selective extinction at the K–Pg boundary. Proc. R. Soc. B Biol. Sci. 283, 20160256 (2016).Article 

Google Scholar 
52.Liu, L. et al. Genomic evidence reveals a radiation of placental mammals uninterrupted by the KPg boundary. Proc. Natl Acad. Sci. USA 114, E7282–E7290 (2017).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
53.Arbour, V. M., Zanno, L. E. & Gates, T. A. Ankylosaurian dinosaur palaeoenvironmental associations were influenced by extirpation, sea-level fluctuation, and geodispersal. Palaeogeogr. Palaeoclimatol. Palaeoecol. 449, 289–299 (2016).Article 

Google Scholar 
54.Tennant, J. P., Mannion, P. D. & Upchurch, P. Sea level regulated tetrapod diversity dynamics through the Jurassic/Cretaceous interval. Nat. Commun. 7, 12737 (2016).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
55.Silvestro, D., Schnitzler, J., Liow, L. H., Antonelli, A. & Salamin, N. Bayesian estimation of speciation and extinction from incomplete fossil occurrence data. Syst. Biol. 63, 349–367 (2014).PubMed 
PubMed Central 
Article 

Google Scholar 
56.Silvestro, D., Antonelli, A., Salamin, N. & Quental, T. B. The role of clade competition in the diversification of North American canids. Proc. Natl Acad. Sci. USA 112, 8684–8689 (2015).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
57.Lehtonen, S. et al. Environmentally driven extinction and opportunistic origination explain fern diversification patterns. Sci. Rep. 7, 4831 (2017).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
58.Condamine, F. L., Romieu, J. & Guinot, G. Climate cooling and clade competition likely drove the decline of lamniform sharks. Proc. Natl Acad. Sci. USA 116, 20584–20590 (2019).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
59.Signor, P. W. & Lipps, J. H. in Geological Implications of Impacts of Large Asteroids and Comets on The Earth (eds. Silver, L. T. & Schultz, P. H.) vol. 190, 291–296 (Geological Society of America Special Publication, 1982).60.Benson, R. B. J. Dinosaur macroevolution and macroecology. Annu. Rev. Ecol. Evol. Syst. 49, 379–408 (2018).Article 

Google Scholar 
61.Dean, C. D., Chiarenza, A. A. & Maidment, S. C. R. Formation binning: a new method for increased temporal resolution in regional studies, applied to the Late Cretaceous dinosaur fossil record of North America. Palaeontology 63, 881–901 (2020).Article 

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

Google Scholar 
63.Condamine, F. L., Rolland, J. & Morlon, H. Assessing the causes of diversification slowdowns: Temperature-dependent and diversity-dependent models receive equivalent support. Ecol. Lett. 22, 1900–1912 (2019).PubMed 
Article 
PubMed Central 

Google Scholar 
64.Prieto-Márquez, A., Dalla Vecchia, F. M., Gaete, R. & Galobart, À. Diversity, relationships, and biogeography of the lambeosaurine dinosaurs from the European archipelago, with description of the new aralosaurin Canardia garonnensis. PLoS ONE 8, e69835 (2013).65.Prieto-Márquez, A., Fondevilla, V., Sellés, A. G., Wagner, J. R. & Galobart, À. Adynomosaurus arcanus, a new lambeosaurine dinosaur from the Late Cretaceous Ibero-Armorican Island of the European archipelago. Cretac. Res. 96, 19–37 (2019).Article 

Google Scholar 
66.Longrich, N. R., Suberbiola, X. P., Pyron, R. A. & Jalil, N.-E. The first duckbill dinosaur (Hadrosauridae: Lambeosaurinae) from Africa and the role of oceanic dispersal in dinosaur biogeography. Cretac. Res. 120, 104678 (2021).Article 

Google Scholar 
67.Kobayashi, Y., Takasaki, R., Kubota, K. & Fiorillo, A. R. A new basal hadrosaurid (Dinosauria: Ornithischia) from the latest Cretaceous Kita-ama Formation in Japan implies the origin of hadrosaurids. Sci. Rep. 11, 8547 (2021).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
68.Stubbs, T. L., Benton, M. J., Elsler, A. & Prieto-Márquez, A. Morphological innovation and the evolution of hadrosaurid dinosaurs. Paleobiology 45, 347–362 (2019).Article 

Google Scholar 
69.Reest, A. J. van der & Currie, P. J. Troodontids (Theropoda) from the Dinosaur Park Formation, Alberta, with a description of a unique new taxon: Implications for deinonychosaur diversity in North America. Can. J. Earth Sci. 54, 919–935 (2017).70.Hartman, S. et al. A new paravian dinosaur from the Late Jurassic of North America supports a late acquisition of avian flight. PeerJ 7, e7247 (2019).PubMed 
PubMed Central 
Article 

Google Scholar 
71.Horner, J. R., Varricchio, D. J. & Goodwin, M. B. Marine transgressions and the evolution of Cretaceous dinosaurs. Nature 358, 59–61 (1992).ADS 
Article 

Google Scholar 
72.O’Brien, C. L. et al. Cretaceous sea-surface temperature evolution: Constraints from TEX86 and planktonic foraminiferal oxygen isotopes. Earth Sci. Rev. 172, 224–247 (2017).ADS 
Article 
CAS 

Google Scholar 
73.Huber, B. T., MacLeod, K. G., Watkins, D. K. & Coffin, M. F. The rise and fall of the Cretaceous Hot Greenhouse climate. Glob. Planet. Change 167, 1–23 (2018).ADS 
Article 

Google Scholar 
74.Mannion, P. D. et al. A temperate palaeodiversity peak in Mesozoic dinosaurs and evidence for Late Cretaceous geographical partitioning. Glob. Ecol. Biogeogr. 21, 898–908 (2012).Article 

Google Scholar 
75.Forster, A., Schouten, S., Baas, M. & Damsté, J. S. S. Mid-Cretaceous (Albian–Santonian) sea surface temperature record of the tropical Atlantic Ocean. Geology 35, 919–922 (2007).ADS 
Article 

Google Scholar 
76.O’Connor, L. K. et al. Late Cretaceous temperature evolution of the southern high latitudes: a TEX86 perspective. Paleoceanogr. Paleoclimatol. 34, 436–454 (2019).ADS 
Article 

Google Scholar 
77.Linnert, C. et al. Evidence for global cooling in the Late Cretaceous. Nat. Commun. 5, 4194 (2014).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
78.Crane, P. R. & Lidgard, S. Angiosperm diversification and paleolatitudinal gradients in Cretaceous floristic diversity. Science 246, 675–678 (1989).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
79.Condamine, F. L., Silvestro, D., Koppelhus, E. B. & Antonelli, A. The rise of angiosperms pushed conifers to decline during global cooling. Proc. Natl Acad. Sci. USA 117, 28867–28875 (2020).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
80.Condamine, F. L., Rolland, J. & Morlon, H. Macroevolutionary perspectives to environmental change. Ecol. Lett. 16, 72–85 (2013).PubMed 
Article 
PubMed Central 

Google Scholar 
81.Silvestro, D., Cascales-Miñana, B., Bacon, C. D. & Antonelli, A. Revisiting the origin and diversification of vascular plants through a comprehensive Bayesian analysis of the fossil record. N. Phytol. 207, 425–436 (2015).Article 

Google Scholar 
82.Prokoph, A., Shields, G. A. & Veizer, J. Compilation and time-series analysis of a marine carbonate δ18O, δ13C, 87Sr/86Sr and δ34S database through Earth history. Earth Sci. Rev. 87, 113–133 (2008).ADS 
CAS 
Article 

Google Scholar 
83.Miller, K. G. et al. The phanerozoic record of global sea-level change. Science 310, 1293–1298 (2005).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
84.Barrett, P. M. Paleobiology of herbivorous dinosaurs. Annu. Rev. Earth Planet. Sci. 42, 207–230 (2014).ADS 
CAS 
Article 

Google Scholar 
85.Grady, J. M., Enquist, B. J., Dettweiler-Robinson, E., Wright, N. A. & Smith, F. A. Evidence for mesothermy in dinosaurs. Science 344, 1268–1272 (2014).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
86.Eagle, R. A. et al. Isotopic ordering in eggshells reflects body temperatures and suggests differing thermophysiology in two Cretaceous dinosaurs. Nat. Commun. 6, 8296 (2015).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
87.Paladino, F. V., Dodson, P., Hammond, J. K. & Spotila, J. R. Temperature-dependent sex determination in dinosaurs? Implications for population dynamics and extinction. in Paleobiology of the Dinosaurs (ed. Farlow, J. O.) vol. 238, 63–70 (Geological Society of America Special Papers, 1989).88.Vavrek, M. J. & Larsson, H. C. E. Low beta diversity of Maastrichtian dinosaurs of North America. Proc. Natl Acad. Sci. USA 107, 8265–8268 (2010).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
89.Dunne, J. A., Williams, R. J. & Martinez, N. D. Network structure and biodiversity loss in food webs: Robustness increases with connectance. Ecol. Lett. 5, 558–567 (2002).Article 

Google Scholar 
90.Brodie, J. F. et al. Secondary extinctions of biodiversity. Trends Ecol. Evol. 29, 664–672 (2014).PubMed 
Article 
PubMed Central 

Google Scholar 
91.Fraser, D. et al. Investigating biotic interactions in deep time. Trends Ecol. Evol. 36, 61–75 (2021).PubMed 
Article 
PubMed Central 

Google Scholar 
92.Mallon, J. C. Competition structured a Late Cretaceous megaherbivorous dinosaur assemblage. Sci. Rep. 9, 15447 (2019).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
93.Benton, M. J. Progress and competition in macroevolution. Biol. Rev. 62, 305–338 (1987).Article 

Google Scholar 
94.Fricke, H. C. & Pearson, D. A. Stable isotope evidence for changes in dietary niche partitioning among hadrosaurian and ceratopsian dinosaurs of the Hell Creek Formation, North Dakota. Paleobiology 34, 534–552 (2008).Article 

Google Scholar 
95.Mallon, J. C. & Anderson, J. S. Skull ecomorphology of megaherbivorous dinosaurs from the Dinosaur Park Formation (Upper Campanian) of Alberta, Canada. PLoS ONE 8, e67182 (2013).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
96.Nordén, K. K., Stubbs, T. L., Prieto-Márquez, A. & Benton, M. J. Multifaceted disparity approach reveals dinosaur herbivory flourished before the end-Cretaceous mass extinction. Paleobiology 44, 620–637 (2018).Article 

Google Scholar 
97.Lyson, T. R. & Longrich, N. R. Spatial niche partitioning in dinosaurs from the latest Cretaceous (Maastrichtian) of North America. Proc. R. Soc. B Biol. Sci. 278, 1158–1164 (2011).Article 

Google Scholar 
98.Li, Z. et al. Ultramicrostructural reductions in teeth: Implications for dietary transition from non-avian dinosaurs to birds. BMC Evol. Biol. 20, 46 (2020).PubMed 
PubMed Central 
Article 

Google Scholar 
99.Cau, A. et al. Synchrotron scanning reveals amphibious ecomorphology in a new clade of bird-like dinosaurs. Nature 552, 395–399 (2017).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
100.Cau, A. The body plan of Halszkaraptor escuilliei (Dinosauria, Theropoda) is not a transitional form along the evolution of dromaeosaurid hypercarnivory. PeerJ 8, e8672 (2020).PubMed 
PubMed Central 
Article 

Google Scholar 
101.Fowler, D. W., Freedman, E. A., Scannella, J. B. & Kambic, R. E. The predatory ecology of Deinonychus and the origin of flapping in birds. PLoS ONE 6, e28964 (2011).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
102.Frederickson, J. A., Engel, M. H. & Cifelli, R. L. Ontogenetic dietary shifts in Deinonychus antirrhopus (Theropoda; Dromaeosauridae): Insights into the ecology and social behavior of raptorial dinosaurs through stable isotope analysis. Palaeogeogr. Palaeoclimatol. Palaeoecol. 552, 109780 (2020).Article 

Google Scholar 
103.O’Connor, J. et al. Microraptor with ingested lizard suggests non-specialized digestive function. Curr. Biol. 29, 2423–2429 (2019).PubMed 
Article 
CAS 
PubMed Central 

Google Scholar 
104.King, J. L., Sipla, J. S., Georgi, J. A., Balanoff, A. M. & Neenan, J. M. The endocranium and trophic ecology of Velociraptor mongoliensis. J. Anat. 237, 861–869 (2020).PubMed 
PubMed Central 
Article 

Google Scholar 
105.Owocki, K., Kremer, B., Cotte, M. & Bocherens, H. Diet preferences and climate inferred from oxygen and carbon isotopes of tooth enamel of Tarbosaurus bataar (Nemegt Formation, Upper Cretaceous, Mongolia). Palaeogeogr. Palaeoclimatol. Palaeoecol. 537, 109190 (2020).Article 

Google Scholar 
106.Dalman, S. & Lucas, S. New evidence for cannibalism in tyrannosaurid dinosaurs from the Late Cretaceous of New Mexico. N. Mex. Mus. Nat. Hist. Sci. Bull. 82, 39–56 (2021).
Google Scholar 
107.Frederickson, J. A., Engel, M. H. & Cifelli, R. L. Niche partitioning in theropod dinosaurs: Diet and habitat preference in predators from the uppermost Cedar Mountain Formation (Utah, U.S.A.). Sci. Rep. 8, 17872 (2018).108.Hassler, A. et al. Calcium isotopes offer clues on resource partitioning among Cretaceous predatory dinosaurs. Proc. R. Soc. B Biol. Sci. 285, 20180197 (2018).109.Schroeder, K., Lyons, S. K. & Smith, F. A. The influence of juvenile dinosaurs on community structure and diversity. Science 371, 941–944 (2021).110.Currie, P. J., Badamgarav, D., Koppelhus, E. B., Sissons, R. & Vickaryous, M. K. Hands, feet, and behaviour in Pinacosaurus (Dinosauria: Ankylosauridae). Acta Palaeontol. Polon. 56, 489–504 (2011).Article 

Google Scholar 
111.Burns, M. E., Currie, P. J., Sissons, R. L. & Arbour, V. M. Juvenile specimens of Pinacosaurus grangeri Gilmore, 1933 (Ornithischia: Ankylosauria) from the Late Cretaceous of China, with comments on the specific taxonomy of Pinacosaurus. Cretac. Res. 32, 174–186 (2011).Article 

Google Scholar 
112.Burns, M. E., Tumanova, T. A. & Currie, P. J. Postcrania of juvenile Pinacosaurus grangeri (Ornithischia: Ankylosauria) from the Upper Cretaceous Alagteeg Formation, Alag Teeg, Mongolia: Implications for ontogenetic allometry in ankylosaurs. J. Paleontol. 89, 168–182 (2015).113.Botfalvai, G., Prondvai, E. & Ősi, A. Living alone or moving in herds? A holistic approach highlights complexity in the social lifestyle of Cretaceous ankylosaurs. Cretac. Res. 118, 104633 (2021).Article 

Google Scholar 
114.Arbour, V. M. & Zanno, L. E. The evolution of tail weaponization in amniotes. Proc. R. Soc. B Biol. Sci. 285, 20172299 (2018).Article 

Google Scholar 
115.Arbour, V. M. & Zanno, L. E. Tail weaponry in ankylosaurs and glyptodonts: An example of a rare but strongly convergent phenotype. Anat. Rec. 303, 988–998 (2020).Article 

Google Scholar 
116.Van Valen, L. A new evolutionary law. Evol. Theory 1, 1–30 (1973).117.Hagen, O., Andermann, T., Quental, T. B., Antonelli, A. & Silvestro, D. Estimating age-dependent extinction: Contrasting evidence from fossils and phylogenies. Syst. Biol. 67, 458–474 (2018).PubMed 
Article 
PubMed Central 

Google Scholar 
118.Finnegan, S., Payne, J. L. & Wang, S. C. The Red Queen revisited: Reevaluating the age selectivity of Phanerozoic marine genus extinctions. Paleobiology 34, 318–341 (2008).Article 

Google Scholar 
119.Doran, N. A., Arnold, A. J., Parker, W. C. & Huffer, F. W. Is extinction age dependent? PALAIOS 21, 571–579 (2006).ADS 
Article 

Google Scholar 
120.Larson, D. W., Brown, C. M. & Evans, D. C. Dental disparity and ecological stability in bird-like dinosaurs prior to the end-Cretaceous mass extinction. Curr. Biol. 26, 1325–1333 (2016).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
121.Romano, M. Disparity versus diversity in ankylosaurid dinosaurs: Explored morphospace indicates two separate evolutive radiations. Rend. Online Soc. Geol. It. 53, 2–8 (2021).122.Turner, A. H., Montanari, S. & Norell, M. A. A new dromaeosaurid from the Late Cretaceous Khulsan locality of Mongolia. Am. Mus. Novitat. 2020, 1–48 (2021).123.Maryańska, T. & Osmólska, H. Pachycephalosauria, a new suborder of ornithischian dinosaurs. Palaeontol. Polon. 30, 45–102 (1974).
Google Scholar 
124.Sereno, P. C. National Geographic Research: Phylogeny of the bird-hipped dinosaurs (Order Ornithischia). Natl Geogr. Res. 2, 234–256 (1986). https://d3qi0qp55mx5f5.cloudfront.net/paulsereno/i/docs/86-NGRes-PhyloOrnithis_1.pdf?mtime=1591821557.125.Sullivan, R. M. A taxonomic review of the Pachycephalosauridae (Dinosauria: Ornithischia). N. Mex. Mus. Nat. Hist. Sci. Bull. 35, 347–365 (2006).
Google Scholar 
126.Lee, M. S. Y., Cau, A., Naish, D. & Dyke, G. J. Morphological clocks in paleontology, and a mid-cretaceous origin of crown aves. Syst. Biol. 63, 442–449 (2014).PubMed 
Article 
PubMed Central 

Google Scholar 
127.Arbour, V. M. & Evans, D. C. A new ankylosaurine dinosaur from the Judith River Formation of Montana, USA, based on an exceptional skeleton with soft tissue preservation. R. Soc. Open Sci. 4, 161086 (2017).128.McDonald, A. T., Wolfe, D. G. & Dooley, A. C. Jr A new tyrannosaurid (Dinosauria: Theropoda) from the Upper Cretaceous Menefee Formation of New Mexico. PeerJ 6, e5749 (2018).PubMed 
PubMed Central 
Article 

Google Scholar 
129.Longrich, N. R. & Field, D. J. Torosaurus is not Triceratops: Ontogeny in chasmosaurine ceratopsids as a case study in dinosaur taxonomy. PLoS ONE 7, e32623 (2012).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
130.Larson, P. L. in Tyrannosaurid Paleobiology (eds. Parrish, J. M., Molnar, R. A., Currie, P. J. & Koppelhus, E. B.) 15–54 (Indiana University Press, 2013).131.Yun, C. Evidence points out that ‘Nanotyrannus’ is a juvenile Tyrannosaurus rex. PeerJ 3, e1052 (2015).Article 

Google Scholar 
132.Brusatte, S. L. et al. Dentary groove morphology does not distinguish ‘Nanotyrannus’ as a valid taxon of tyrannosauroid dinosaur. Comment on: “Distribution of the dentary groove of theropod dinosaurs: Implications for theropod phylogeny and the validity of the genus Nanotyrannus Bakker et al., 1988. Cretac. Res. 65, 232–237 (2016).Article 

Google Scholar 
133.Schmerge, J. D. & Rothschild, B. M. When a groove is not a groove: Clarification of the appearance of the dentary groove in tyrannosauroid theropods and the distinction between Nanotyrannus and Tyrannosaurus. Reply to Comment on: “Distribution of the dentary groove of theropod dinosaurs: Implications for theropod phylogeny and the validity of the genus Nanotyrannus Bakker et al., 1988. Cretac. Res. 65, 238–243 (2016).Article 

Google Scholar 
134.Xu, X., Zhou, Z., Sullivan, C., Wang, Y. & Ren, D. An updated review of the Middle-Late Jurassic Yanliao biota: Chronology, taphonomy, paleontology and paleoecology. Acta Geol. Sin. 90, 2229–2243 (2016).Article 

Google Scholar 
135.Cau, A., Brougham, T. & Naish, D. The phylogenetic affinities of the bizarre Late Cretaceous Romanian theropod Balaur bondoc (Dinosauria, Maniraptora): Dromaeosaurid or flightless bird? PeerJ 3, e1032 (2015).PubMed 
PubMed Central 
Article 

Google Scholar 
136.Agnolin, F. L. & Motta, M. J. Paravian phylogeny and the dinosaur-bird transition: An overview. Front. Earth Sci. 6, 252 (2019).ADS 
Article 

Google Scholar 
137.Pei, R. et al. Potential for powered flight neared by most close avialan relatives, but few crossed Its thresholds. Curr. Biol. 30, 4033–4046 (2020).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
138.Foth, C. & Rauhut, O. W. M. Re-evaluation of the Haarlem Archaeopteryx and the radiation of maniraptoran theropod dinosaurs. BMC Evol. Biol. 17, 236 (2017).139.Rauhut, O. W., Tischlinger, H. & Foth, C. A non-archaeopterygid avialan theropod from the Late Jurassic of southern Germany. eLife 8, e43789 (2019).PubMed 
PubMed Central 
Article 

Google Scholar 
140.Lefèvre, U. et al. A new Jurassic theropod from China documents a transitional step in the macrostructure of feathers. Sci. Nat. 104, 74 (2017).Article 
CAS 

Google Scholar 
141.Shen, C. et al. A new troodontid dinosaur from the Lower Cretaceous Yixian formation of Liaoning province. China Acta Geol. Sin. 91, 763–780 (2017).Article 

Google Scholar 
142.Arbour, V. M. & Currie, P. J. Euoplocephalus tutus and the diversity of ankylosaurid dinosaurs in the Late Cretaceous of Alberta, Canada, and Montana, USA. PLoS ONE 8, e62421 (2013).143.Arbour, V. M. & Currie, P. J. Systematics, phylogeny and palaeobiogeography of the ankylosaurid dinosaurs. J. Syst. Palaeontol. 14, 385–444 (2016).Article 

Google Scholar 
144.Arbour, V. M., Currie, P. J. & Badamgarav, D. The ankylosaurid dinosaurs of the Upper Cretaceous Baruungoyot and Nemegt formations of Mongolia. Zool. J. Linn. Soc. 172, 631–652 (2014).
Google Scholar 
145.Arbour, V. M. et al. A new ankylosaurid dinosaur from the Upper Cretaceous (Kirtlandian) of New Mexico with implications for ankylosaurid diversity in the Upper Cretaceous of Western North America. PLoS ONE 9, e108804 (2014).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
146.Gradstein, F. M., Ogg, J. G., Schmitz, M. D. & Ogg, G. M. The Geologic Time Scale 2012 (Elsevier B.V., 2012).147.Brown, C. M. & Henderson, D. M. A new horned dinosaur reveals convergent evolution in cranial ornamentation in Ceratopsidae. Curr. Biol. 25, 1641–1648 (2015).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
148.Jerzykiewicz, T., Currie, P. J., Fanti, F. & Lefeld, J. Lithobiotopes of the Nemegt Gobi Basin. Can. J. Earth Sci. https://doi.org/10.1139/cjes-2020-0148 (2021).149.Silvestro, D., Salamin, N. & Schnitzler, J. PyRate: A new program to estimate speciation and extinction rates from incomplete fossil data. Methods Ecol. Evol. 5, 1126–1131 (2014).Article 

Google Scholar 
150.Rambaut, A. R., Drummond, A. J., Xie, D., Baele, G. & Suchard, M. A. Posterior summarization in Bayesian phylogenetics using Tracer 1.7. Syst. Biol. 67, 901–904 (2018).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
151.Brusatte, S. L. et al. Tyrannosaur paleobiology: New research on ancient exemplar organisms. Science 329, 1481–1485 (2010).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
152.Ryan, M. J., Chinnery-Allgeier, B. J. & Eberth, D. A. New Perspectives on Horned Dinosaurs (Indiana University Press, 2010).153.Xu, X., Wang, K., Zhao, X. & Li, D. First ceratopsid dinosaur from China and its biogeographical implications. Chin. Sci. Bull. 55, 1631–1635 (2010).CAS 
Article 

Google Scholar 
154.Hannisdal, B. & Peters, S. E. Phanerozoic Earth system evolution and marine biodiversity. Science 334, 1121–1124 (2011).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
155.Liow, L. H., Reitan, T. & Harnik, P. G. Ecological interactions on macroevolutionary time scales: Clams and brachiopods are more than ships that pass in the night. Ecol. Lett. 18, 1030–1039 (2015).PubMed 
Article 
PubMed Central 

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

Google Scholar 
157.Mayhew, P. J., Bell, M. A., Benton, T. G. & McGowan, A. J. Biodiversity tracks temperature over time. Proc. Natl Acad. Sci. USA 109, 15141–15145 (2012).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
158.Zachos, J., Pagani, M., Sloan, L., Thomas, E. & Billups, K. Trends, rythms, and aberration in global climate 65 Ma to present. Science 292, 686–693 (2001).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
159.Zachos, J. C., Dickens, G. R. & Zeebe, R. E. An early Cenozoic perspective on greenhouse warming and carbon-cycle dynamics. Nature 451, 279–283 (2008).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
160.Cramer, B. S., Toggweiler, J. R., Wright, J. D., Katz, M. E. & Miller, K. G. Ocean overturning since the late cretaceous: Inferences from a new benthic foraminiferal isotope compilation. Paleoceanography 24, 1–14 (2009).Article 

Google Scholar 
161.Barba-Montoya, J., Reis, M., Schneider, H., Donoghue, P. C. J. & Yang, Z. Constraining uncertainty in the timescale of angiosperm evolution and the veracity of a Cretaceous Terrestrial Revolution. N. Phytol. 218, 819–834 (2018).Article 

Google Scholar 
162.Zhang, M., Dai, S., Du, B., Ji, L. & Hu, S. Mid-Cretaceous hothouse climate and the expansion of early angiosperms. Acta Geol. Sin. 92, 2004–2025 (2018).Article 

Google Scholar  More

Filtering threshold for handling spurious sequencesWe first used bacterial communities of known composition (simplified communities) to assess the occurrence of spurious taxa and to determine at which relative abundances they begin to appear. To propose a cutoff that is potentially applicable to different 16S rRNA gene amplicon studies, we included reference data obtained with different variable regions and sequencing pipelines and originating from both in vitro an in vivo communities varying in number and type of species (max. 58) (Tables 1 and 2). To determine a filtering threshold that allowed exclusion of most spurious taxa, we recorded the relative abundance of the first spurious OTU occurring in each of the reference community datasets (Fig. 2a). Median values of approx. 0.12% relative abundance were observed (Fig. 2b). Besides one outlier in the mock communities (0.44% relative abundance), all values were below 0.25% relative abundance.Fig. 2: Determination of filtering thresholds using artificial communities of known composition in vitro (mock; n = 9 different types; 21 replicates in total) and in mice (gnotobiotes; n = 4 different communities; 28 mice in total).a Example of the occurrence of all molecular species detected without filtering in the gut of a gnotobiotic mouse [49]. The arrow indicates the position of the first spurious molecular species, all following taxa being considered as having a high risk of being spurious (light gray bars in the enlarged inset). b Distribution of the relative abundances of first occurring spurious molecular species (as shown in panel a) across all mock communities and samples from gnotobiotes. The orange dashes on the y-axis indicate the consensus threshold of 0.25% relative abundance, above which no spurious taxa occurred with the exception of one outlier in a mock community at a relative abundance of 0.44%. c Comparison of various standard filtering cutoffs (see explanations in the text) in terms of spurious taxa (i.e., those molecular species not matching sequences of the known species contained in the artificial communities). d Corresponding percentages of positive hits retained by the different filtering strategies, with positive hits being defined as the reference sequences found in the respective amplicon datasets. e Percentage of spurious taxa and positive hits in the same reference communities using the DADA2 pipeline for analysis based on amplicon sequence variants (ASVs) [6]. f Effect of filtering thresholds at increments of 0.05% relative abundance on the detection of spurious taxa and positive hits in all mock and gnotobiotic datasets for OTUs (upper panel) and ASVs (lower panel). Lines correspond to mean values; ribbons represent standard deviations.Full size imageWithout any filtering, sequence clustering generated an average of 508 ± 355 OTUs (min. 52; max. 1081) per mock community (10–58 target species in theory) and 105 ± 50 OTUs (min, 55; max. 215) per gnotobiotic community (4–12 target species in theory). Up to 87% of these OTUs were spurious (i.e., they did not match the expected classification of species contained in the corresponding artificial community) (Fig. 2c). On average, the proportion of spurious OTUs in both the mock communities and samples from gnotobiotic mice was slightly lower after removing singletons, although this did not reach statistical significance (50.8 vs. 64.3%, p = 0.227; 57.5% vs. 65.7%; p = 0.70, pairwise comparison by t-test, including Benjamini–Hochberg correction following ANOVA). Interestingly, the proportion of spurious molecular species was higher in gnotobiotic mice independent of filtering (p  0.50) (Fig. 2d). Note that the diversity of reference communities in the gnotobiotic mice was relatively low (4–12 members; Table 2), resulting in a marked drop in the percentage of positive hit (8–25%) when even just one true member is excluded after filtering because of its low relative abundance (which is an expectable event considering a classical, exponentially decreasing distribution of species occurrence in gut environments).We next employed the widely used ASV analysis approach to confirm the aforementioned results. Processing of the same simplified communities generated a total number of 42 ± 25 ASVs (min. 16; max. 98) for mock communities (10–58 target species) and 14 ± 8 ASVs (min. 4; max. 25) for gnotobiotes (4–12 target species). Altogether, a marked decrease in spurious taxa was observed compared with OTU clustering, with an average of 8.6 ± 11.8 and 4.4 ± 6.4% spurious sequences for mock and gnotobiotic communities, respectively (comparison of purple box plots in Fig. 2e, top panels, and Fig. 2c). Of note, the DADA2 pipeline used for the ASV approach does not infer sequence variants that are only supported by a single read (singletons) due to a lack of confidence in their existence relative to sequencing errors. Consequently, data corresponding to “no filtering” with the OTU-based approach were not generated. On average, the first spurious ASV occurred at a relative abundance of 0.10 ± 0.32%. By applying the cutoff of 0.25% relative abundance, spurious sequences were completely removed (except for three outlying samples), albeit with a slight drop in positive hits for both mock and gnotobiotic communities (Fig. 2e).To obtain a more comprehensive view on how filtering thresholds affect the detection of spurious taxa, all datasets (mock and gnotobiotic mice) were processed using a range of relative abundance filtering thresholds (from 0 to 0.5% at increments of 0.05%) after either OTU- or ASV-based processing of raw sequence reads (Fig. 2f). These data indicate that filtering thresholds between 0.1 and 0.3% are appropriate to reduce the occurrence of spurious taxa to 600 of the 678 spurious OTUs occurred in fewer than five of the ten sequencing runs tested, with approximately 450 of them occurring in only one run (Fig. 3c). This observation indicates that the majority of spurious taxa are sporadic cross-contaminations rather than generalist artifacts across sequencing runs, suggesting that fully independent technical replicates would improve data quality. Although most of the spurious taxa were characterized by relative abundances between 0.25 and 2% in the IMNGS-amplicon datasets tested, they represented very dominant populations in a few samples (Fig. 3d).Fig. 3: Origin and occurrence of spurious taxa.a Taxonomic profile and ecological distribution. Inner ring: SILVA-based classification of all non-redundant spurious molecular species at the phylum and family level. Outer colored ring: sample type characterized by the highest prevalence for the given taxon. Outer bars: corresponding highest prevalence values. Only samples with relative abundances >0.25% for any given OTU were counted as positive for prevalence calculation. The total numbers of samples considered were: human, 46,153; soil, 29,864; freshwater, 13,977; mouse, 10,409; marine, 8478. b Distribution of the spurious taxa across sample types. The exclusivity of each OTU for any given sample type was assessed using a Z-test: those assumed to be non-specific for any given sample type appear in red (p 0.25% in at least one replicate were kept). Richness was calculated using ampvis2 [29]. Applying the 0.25% cutoff decreased the number of observed ASVs from 408 ± 71 to 139 ± 5 and, more importantly, the IQR from 101 to 7 (Fig. 6b). Unweighted UniFrac distances within and between runs as calculated using ampvis2 were also compared before and after filtering. Sequences were aligned using MAFFT [30] and phylogeny was inferred using FastTree. Whilst the community makeup in the soil sample varied substantially between sequencing runs without additional filtering, the 0.25% cutoff reduced this variation to the level observed within runs without filtering (Fig. 6c). Replicates within a run were very similar after applying the 0.25% cutoff. Altogether, these data serve as an independent confirmation that stringent filtering delivers more stable values obtained for the exact same sample sequenced in replicates across several sequencing runs. More

2,809 draft MAGs from the rumen ecosystemWe amassed 3.2 terabase pairs (Tbp) of data from 346 publicly available and 66 new rumen metagenome datasets (Supplementary Table 1). The metagenomes were from cattle (312 samples, 2.1 Tbp), sheep (75 samples, 888.4 gigabase pairs (Gbp)), moose (9 samples, 108.8 Gbp), deer (8 samples, 62.9 Gbp), and bison (8 samples, 52.3 Gbp). Metagenomes were assembled independently to reduce the influence of strain variation and improve the recovery of closely related genomes18,19. Following refinement, dereplication, and filtering of resulting population genomes, we identified 2,809 nonredundant MAGs satisfying the following criteria: dRep20 genome quality score ≥60, ≥75% complete, ≤10% contamination, N50 ≥5 kbp, and ≤500 contigs.The median estimated completeness and contamination of the MAGs were 89.7% and 0.9%, respectively (Fig. 1a and Supplementary Data 1). Further, recovered MAGs had a median genome size of 2.2 Mbp, a median of 131 contigs, and a median N50 of 28.3 kbp (Fig. 1b). The proposed minimum information about a MAG (MIMAG) specifies high-quality draft genomes to have an estimated ≥90% completeness, ≤5% contamination, at least 18 tRNAs, and contain 23S, 16S, and 5S rRNA genes21. It remains challenging to reconstruct rRNA genes from short metagenomic reads due to the high sequence similarity of rRNA genes in closely related species. As a result, despite high estimated completeness and low contamination rates, only 20 MAGs meet the MIMAG standards for a high-quality draft genome. We identified a 16S rRNA gene in 197 of the MAGs. The remaining MAGs are characterized as medium-quality MAGs under the MIMAG standards.Fig. 1: Genomic properties of 2,809 rumen MAGs.a CheckM completeness and contamination estimates for the 2,809 population genomes recovered from rumen metagenomes. The size of the point on the scatter plot corresponds to the dRep genome quality score, where Quality = Completeness − (5 ⋅ Contamination) + (Contamination ⋅ (Strain Heterogeneity/100)) + 0.5 ⋅ (({mathrm{log}},)(N50). The reported MAGs meet the following minimum criteria: genome quality score ≥60, ≥75% complete, ≤10% contamination, N50 ≥5 kbp, and ≥500 contigs. b The frequency distribution of the number of contigs and genome sizes of reconstructed MAGs.Full size imageThe majority of bacterial MAGs belonged to phyla Firmicutes or Bacteroidota (2,326; Fig. 2a and Supplementary Data 1). Additionally, we assembled 12 bacterial genomes from the superphylum Patescibacteria. At lower taxonomic ranks, Lachnospiraceae (415) and Prevotella (398) were the dominant family and genus identified among the assembled bacterial genomes. The most prevalent archaeal family and genus were Methanobacteriaceae (45) and Methanobrevibacter (35), respectively (Fig. 2b). The recovered MAGs represent several new taxonomic lineages, as four genomes could not be classified at the rank of order, 16 at the rank of family, and 243 at the genus rank.Fig. 2: Phylogenetic relationships and coverage patterns of near-complete MAGs.a Phylogenomic analysis of 1,163 near-complete (≥90% complete, ≤5% contamination, and N50 ≥15 kbp) bacterial MAGs and (b) 20 near-complete archaeal MAGs inferred from the concatenation of phylogenetically informative proteins. Layers below the genomic trees designate bacterial phylum or archaeal genus based on GTDB taxonomic assignments, genomic size (0–5 Mbp), and the mean number of bases with ≥1× coverage in a rumen metagenomic dataset (layer color indicates the ruminant the data was collected from). The mean number of bases with ≥1× coverage was used as input for hierarchical clustering of rumen metagenomic datasets based on Euclidean distance and Ward linkage. The bacterial and archaeal phylogenetic trees are provided as Supplementary Data 6 and Supplementary Data 7, respectively.Full size imageSpecies-level overlap between reference genomes, the Hungate1000 Collection, and rumen MAGsTo further characterize the assembled genomes, we compared the MAGs to other rumen-specific genome collections, specifically genomes generated from the Hungate1000 project3 and MAGs identified from the Stewart et al. studies4,5. We clustered genomes based on approximate species-level thresholds (≥95% ANI) and calculated the intersection between MAGs in the current study and the Hungate1000 Collection (410 genomes)3, MAGs from Stewart et al. (4,941 genomes)4,5, and a dereplicated genome collection from the GTDB (22,441 genomes, see Methods)22, which includes reference isolate genomes and some environmental MAGs23. It should be noted that we used the raw data from the first of the Stewart et al. studies4 (Supplementary Table 1), but with different assembly and binning approaches. Approximately one-third of the MAGs (1,007) did not exhibit ≥95% ANI with a genome in the GTDB, Stewart et al. MAGs, or the Hungate1000 isolates (Fig. 3a). When considering the pairwise intersections between the datasets, 98 (3.5%), 933 (33.2%), and 1,438 (51.2%) of the MAGs in the current study had ≥95% ANI with a genome in the Hungate1000 Collection3, GTDB22, and Stewart et al.4,5, respectively. One hundred twenty-one (29.5%), 552 (2.5%), and 3,125 (63.2%) of the genomes from the Hungate1000 Collection3, GTDB22, and Stewart et al.4,5 displayed ≥95% ANI with a MAG from the current study. Together, these results indicate that we recovered a majority of previous rumen genomic diversity with additional lineages not previously identified in other major rumen genomic collections.Fig. 3: Genomes sharing ≥95% ANI between databases and the characterization of rumen-specific 95% ANI clusters.a The approximate number of species overlapping amongst rumen-specific and reference genomic datasets. Genomes demonstrating ≥95% ANI were considered to be shared between two datasets. Presented are a subset of intersections in which a MAG from the current study was the query genome. b The number of genomes comprising each of the 3,541 95% ANI clusters generated from 8,160 rumen microbial genomes in the current study, the Hungate1000 Collection3, and Stewart et al. studies4, 5. c Rarefaction analysis based on subsampling 95% ANI clusters at steps of 500 genomes indicates the 8,160 genomes from recently published rumen genomic collections still only represent a fraction of expected microbial species diversity in the rumen ecosystem. Phylogenomic relationships of the 1,781 near-complete bacterial (d) and 35 near-complete archaeal (e) representative genomes with the highest dRep genome quality score from the 3,541 95% ANI clusters generated from 8,160 rumen-specific genomes. Near-complete genomes were defined as being ≥90% complete, having ≤5% contamination, and contig N50 ≥15 kbp. Layers surrounding the genomic trees indicate the bacterial phyla or archaeal genera and the log normalized number of genomes from each rumen genomic collection belonging to the same 95% ANI cluster. The bacterial and archaeal phylogenetic trees are provided as Supplementary Data 8 and Supplementary Data 9, respectively.Full size imageWe applied an additional clustering approach to identify the approximate number of species represented by the rumen-specific genomes assembled in this study, in the Hungate1000 Collection3, and Stewart et al.4,5. A 95% ANI threshold yielded 3,541 clusters from the combination of the datasets (Supplementary Data 2). Of the 3,541 clusters, 2,024 contained a MAG from the current study, and 1,135 were composed exclusively of MAGs from the current study. In comparison, 2,175 and 286 clusters were comprised of genomes from Stewart et al.4,5 and the Hungate1000 Collection3, respectively. The majority of 95% ANI clusters (2,166) are only comprised of a single genome (Fig. 3b). Furthermore, a rarefaction curve suggests the 8,160 genomes from the genomic collections analyzed here only represent a fraction of the estimated microbial species diversity in the rumen (Fig. 3c). The genome with the best dRep score from each cluster was used to generate a phylogenetic tree highlighting the species diversity within each rumen genomic collection and represents the vast diversity of rumen bacterial (Fig. 3d) and archaeal (Fig. 3e) genomes published to date.As stated previously, the median genome size of reconstructed MAGs was 2.2 Mbp, smaller than the median size of genomes from the Hungate1000 project (3.1 Mbp)3. To provide an assessment at a finer resolution, genome sizes of MAGs and Hungate1000 genomes3 belonging to the same 95% ANI cluster were compared (Supplementary Fig. 1). Adjusted sizes of MAGs and Hungate1000 genomes that are ≥95% complete displayed a regression coefficient of 0.96 with a slope of 0.86, indicating the binning process likely did not lead to extensive losses and systematic biases in the reconstructed genomes. Instead, it further highlights that current culturing approaches have not brought large portions of rumen microbial diversity into culture and putatively supports previous findings from the human gut that revealed genome-reduction in uncultured bacteria24.Rumen metagenome classification rates using reference and rumen-specific genomesUtilizing an approach similar to Stewart et al.4,5, we investigated the influence of MAGs on rates of metagenomic read classification. The baseline for read classification was the standard Kraken database containing bacterial, archaeal, fungal, and protozoal RefSeq genomes25. Each rumen-specific dataset was incrementally added to the Kraken RefSeq genomic database in the following order to build new databases: the Hungate1000 Collection3, MAGs from Stewart et al.4,5, and MAGs from the current study. Each individual and collective database was used for classification of sample reads that underpinned metagenomic binning and from a rumen metagenomic dataset not used in the reconstruction of MAGs26. MAGs from the current work classified more reads from deer, moose, and sheep metagenomes, while the more numerous MAGs from Stewart et al.4,5 classified more reads from bison and cattle metagenomes (Supplementary Fig. 2a). The addition of MAGs improves classification relative to databases primarily based on cultured isolates, like the Hungate1000 Collection3 (Supplementary Fig. 2b). Using the combination of all reference and rumen-specific genomes, the median classification rate on an independent set of cattle metagenomes was 62.6%.Phylogenetic characterization of biosynthetic gene clustersMicrobial genome mining is a powerful tool for natural product discovery. We sought to explore the extent of secondary metabolite diversity coded by the MAGs in the current study, the Hungate1000 Collection3, and Stewart et al. MAGs4,5. We identified 14,814 BGCs encoded by the 8,160 rumen-specific genomes using antiSMASH27 (Fig. 4a and Supplementary Data 3). The majority of BGCs were NRPS (5,346), followed by aryl polyenes (2,800), sactipeptides (2,126), and bacteriocins (1,943). Only a few PKS were identified (75). Firmicutes harbored the vast majority of clusters for NRPS, sactipeptide, lantipeptide, lassopeptide, and bacteriocin synthesis (Fig. 4b). At lower taxonomic ranks, DTU089 (979), Bacteroidaceae (934), and Lachnospiraceae (923) coded for the bulk of NRPS gene clusters. Moreover, Acidaminococcaceae genomes contained 21.2% of identified bacteriocins and Ruminococcus spp. possessed the bulk of sactipeptides and lantipeptides. Archaea were predicted to code 737 BGCs, including an average of 3.8 NRPS gene clusters per genome (Fig. 4a).Fig. 4: Characterization of BGCs from 8,160 rumen genomes and MAGs.a Number and types of BGCs identified from select phyla in genomes from the Hungate1000 Collection3, Stewart et al. studies4, 5, and the current study. b Phylogenomic analysis of 1,766 near-complete Firmicutes genomes inferred from the concatenation of phylogenetically informative proteins. The inner layer surrounding the genomic tree designates taxonomic annotations, while the remaining layers depict the log normalized number of BGCs in the genome with the ascribed function. Bacterial class and order labels are displayed for those lineages in which more than 50 genomes were identified. Near-complete genomes were defined as being ≥90% complete, having ≤5% contamination, and contig N50 ≥15 kbp. The phylogenetic tree is provided as Supplementary Data 10. c A relational network of NRPS gene clusters in Firmicutes, Bacteroidota, and Euryarchaeota highlights the similarity of NRPS BGCs from Euryarchaeota and Firmicutes. Edge weight represents the similarity of two BGCs, as determined by BiG-SCAPE (i.e. darker edges demonstrate more similarity between two BGCs). Edges are only shown for BGCs with ≥0.3 BiG-SCAPE similarity. Nodes from each phylum are duplicated to illustrate intra-phylum relationships and nodes along a given axis are ordered alphabetically by taxonomic family. d The association between genome phylogeny and the similarity of NRPS gene clusters coded by near-complete Euryarchaeota genomes. BGCs designated as NRPS were clustered with BiG-SCAPE. The relationship between NRPS clusters was portrayed through the hierarchical clustering of pairwise inter-cluster similarities. The number of NRPS clusters coded by each genome (range of 0–3) is presented alongside the assigned genus. A group of Methanobrevibacter genomes, likely of the same species (≥95% ANI), possessed very similar NRPS clusters (highlighted in blue). Yet, phylogenetically closely related genomes, belonging to two different 95% ANI clusters, did not code for any identified NRPS gene clusters (highlighted in red). The phylogenetic tree is based on the concatenation of 122 phylogenetically informative archaeal proteins and is available as Supplementary Data 11.Full size imageNRPS exhibit high molecular and structural diversity resulting in a wide array of biological activities. The diversity of NRPS, combined with their proteolytic stability and selective bioactivity, has resulted in the development of many NRPS as antimicrobials and other therapeutic agents28. Given the prevalence of NRPS among the recovered MAGs (Fig. 4a), the peptides appear to be important bioactive metabolites in the rumen. To gain fundamental insight into the phylogenetic diversity of rumen NRPS, we built a network based on BGC similarity using BiG-SCAPE29. BiG-SCAPE uses protein domain content, order, copy number, and sequence identity to calculate a distance metric. We assessed the similarity of NRPS gene clusters identified in Firmicutes, Bacteroidota, and Euryarchaeota, as these three phyla coded for 96.4% of assembled NRPS gene clusters from rumen genomes. With a BiG-SCAPE similarity threshold of 0.3, the resulting network consisted of 3,436 nodes (NRPS BGCs on contigs ≥10 kbp) and 79,112 edges (Fig. 4c and Supplementary Data 4). As expected, the network analysis depicted high inter- and intra-phylum genetic diversity among the NRPS gene clusters. The median intra-phylum, -family, and -genus similarity was 0.40, 0.44, and 0.46, respectively, while the median inter-phylum, -family, and -genus similarity was 0.32, 0.34, and 0.34, respectively. Further, only 2.6% of edges were inter-phylum and 69.0% were intra-family. Of the 6,594 Euryarchaeota edges, 8.1% were Euryarchaeota-Firmicutes (median similarity of 0.32) and 2.0% of edges were Euryarchaeota-Bacteroidota (median similarity of 0.31). To further examine the phylogenetic relationships of rumen Euryarchaeota NRPS, we clustered 265 NRPS gene clusters (≥10 kbp) from 85 near-complete Euryarchaeota genomes at a higher similarity threshold of 0.75, yielding 57 NRPS clusters (Fig. 4d). The distribution of NRPS clusters amongst the genomes suggests there exists a strong relationship between methanogen phylogeny and NRPS similarity. Only Methanobrevibacter genomes contain NRPS gene clusters, and genomes of the same species often possessed many of the same NRPS clusters (see genomes highlighted in blue in Fig. 4d). However, there are instances in which closely related methanogens code for a contrasting pattern of NRPS clusters or no NRPS clusters at all (see genomes highlighted in red in Fig. 4d).Bacteriocins likely serve as regulatory elements in complex microbial communities such as the rumen. Consequently, bacteriocins have been studied and characterized for their bactericidal activity and as agents that modulate the microbiota structure and function30. In particular, lanthipeptides, a class of ribosomally synthesized and post-translationally modified peptides (RiPPs) with thioether cross-linked amino acids31, are of pharmaceutical, preservative, and agricultural interest due to their strong antimicrobial properties against gram-positive pathogens31,32,33, low levels of antimicrobial resistance34, and stability35. We identified 195 rumen lanthipeptide BGCs from the Hungate1000 genomes and MAGs from Stewart et al. and the current study. Rumen lanthipeptide BGCs were clustered with 22,870 lanthipeptide BGCs from RefSeq genomes36,37 into gene cluster families (GCFs; groups of BGCs that may generate highly similar products). Clustering with BiG-SCAPE29 yielded 4,565 GCFs, 120 of which contained a rumen lanthipeptide. The 120 GCFs were composed of 519 lanthipeptide BGCs, where 324 were from RefSeq isolates and 195 from rumen genomes (Fig. 5a). The 324 RefSeq BGCs fell into only 18 GCFs. Lanthipeptides from the Hungate1000 isolates clustered into 36 GCFs, while rumen MAG lanthipeptides belonged to 92 GCFs, 82 of which were exclusively composed of MAG lanthipeptides. Together, this evidence suggests rumen MAGs code for diverse and novel lanthipeptides not represented in cultured isolates, including the Hungate Collection.Fig. 5: Phylogenetic diversity of 195 lanthipeptide BGCs coded by rumen genomes.a Network depicting the similarity between lanthipeptide BGCs identified from complete and draft isolate genomes in RefSeq and rumen genomes of the Hungate1000 collection, Stewart et al. MAGs, and MAGs from the current study. The BGCs were clustered into gene cluster families (GCFs) with BiG-SCAPE29. Only the GCFs containing a rumen genome and at least two BGCs were visualized. Nodes in the network represent BGCs and edges connect BGCs with BiG-SCAPE defined similarity ≥0.3. b Phylogenetic relationships of 120 near-complete rumen bacterial genomes coding for lanthipeptide BGCs. Near-complete genomes were defined as being ≥90% complete, having ≤5% contamination, and contig N50 ≥15 kbp. Layers surrounding the genomic trees indicate the bacterial phyla and family, if the genome is a MAG or Hungate Collection isolate, and the class of lanthipeptide, as predicted by antiSMASH27. Genomes without an indicated lanthipeptide class were not classified by antiSMASH. The phylogenetic tree is based on the concatenation of 120 phylogenetically informative bacterial proteins and is available as Supplementary Data 12.Full size imageWe sought to further examine the differences in rumen MAG lanthipeptides relative to isolates and the taxonomic diversity of rumen microbes coding for lanthipeptides. The 195 rumen lanthipeptides were mainly found in Firmicutes genomes, with a subset from Bacteroidota and Actinobacteriota (Fig. 5b). Fifty-two of the 55 lanthipeptides from the Hungate Collection isolates were from Firmicutes (94.5%). At the family-level, these 52 Firmicutes BGCs were distributed evenly between Lachnospiraceae and Streptococcaceae. In contrast, 19.2% and 8.6% of lanthipeptides from rumen MAGs belonged to Bacteroidota and Actinobacteriota, respectively. Lanthipeptides from MAGs were also found in Muribaculaceae and Oscillospiraceae. Moreover, 26.4% of rumen MAG lanthipeptides, compared to 3.6% of Hungate Collection isolates, were found in Eubacterium genomes. The majority of Eubacterium MAG lanthipeptides (62.1%) belonged to a single GCF, suggesting they code for very similar products. Lastly, antiSMASH predicted the bulk of the rumen lanthipeptides were Class II lanthipeptides, with fewer Class I and Class III types (Fig. 5b). Nearly all of the Class I lanthipeptides were from Hungate isolates. The above analysis of lanthipeptide diversity further supports that rumen MAGs code for novel secondary metabolites not represented in cultured isolates.We aligned previously published rumen metatranscriptome data from steers characterized as having high and low feed efficiency to the BGCs to demonstrate if the identified BGCs are active and to explore potential ecological roles of secondary metabolites. Despite data from the metatranscriptome study not being applied to reconstruct genomes in the current study, we identified the expression of 554 gene clusters from rumen-specific genomes in the 20 metatranscriptomes (≥100 aligned reads). Metatranscriptome read count data were normalized independently for each genome to better account for the variation in taxonomic composition across samples38. Genome-specific normalization resulted in the identification of 17 differentially expressed gene clusters between steers with high and low feed efficiency (DESeq239 false discovery rate adjusted P  More

1.Bálint M, Domisch S, Engelhardt CHM, Haase P, Lehrian S, Sauer J, et al. Cryptic biodiversity loss linked to global climate change. Nat Clim Chang. 2011;1:313–8.Article 

Google Scholar 
2.Parmesan C, Yohe G. A globally coherent fingerprint of climate change impacts across natural systems. Nature. 2003;421:37–42.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
3.Blois JL, Zarnetske PL, Fitzpatrick MC, Finnegan S. Climate change and the past, present, and future of biotic interactions. Science. 2013;341:499–504.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
4.Haines A, Ebi K. The imperative for climate action to protect health. N. Engl J Med. 2019;380:263–73.PubMed 
Article 
PubMed Central 

Google Scholar 
5.Deutsch CA, Tewksbury JJ, Tigchelaar M, Battisti DS, Merrill SC, Huey RB, et al. Increase in crop losses to insect pests in a warming climate. Science. 2018;361:916–9.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
6.Kattwinkel M, Jan-Valentin K, Foit K, Liess M. Climate change, agricultural insecticide exposure, and risk for freshwater communities. Ecol Appl. 2011;21:2068–81.PubMed 
Article 
PubMed Central 

Google Scholar 
7.Moe SJ, De Schamphelaere K, Clements WH, Sorensen MT, Van den Brink PJ, Liess M. Combined and interactive effects of global climate change and toxicants on populations and communities. Environ Toxicol Chem. 2013;32:49–61.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
8.Potts SG, Biesmeijer JC, Kremen C, Neumann P, Schweiger O, Kunin WE. Global pollinator declines: trends, impacts and drivers. Trends Ecol Evol. 2010;25:345–53.PubMed 
Article 
PubMed Central 

Google Scholar 
9.Moran EV, Alexander JM. Evolutionary responses to global change: Lessons from invasive species. Ecol Lett. 2014;17:637–49.PubMed 
Article 
PubMed Central 

Google Scholar 
10.Harwood AD, You J, Lydy MJ. Temperature as a toxicity identification evaluation tool for pyrethroid insecticides: toxicokinetic confirmation. Environ Toxicol Chem. 2009;28:1051–8.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
11.Guo L, Su M, Liang P, Li S, Chu D. Effects of high temperature on insecticide tolerance in whitefly Bemisia tabaci (Gennadius) Q biotype. Pestic Biochem Physiol. 2018;150:97–104.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
12.Mao K, Jin R, Li W, Ren Z, Qin X, He S, et al. The influence of temperature on the toxicity of insecticides to Nilaparvata lugens (Stål). Pestic Biochem Physiol. 2019;156:80–86.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
13.Verheyen J, Delnat V, Stoks R. Increased daily temperature fluctuations overrule the ability of gradual thermal evolution to offset the increased pesticide toxicity under global warming. Environ Sci Technol. 2019;53:4600–8.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
14.Moran NA. Symbiosis as an adaptive process and source of phenotypic complexity. Proc Natl Acad Sci USA. 2007;104:8627–3863.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
15.Kikuchi Y, Hayatsu M, Hosokawa T, Nagayama A, Tago K, Fukatsu T. Symbiont-mediated insecticide resistance. Proc Natl Acad Sci USA. 2012;109:8618–22.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
16.Jones RM, Desai C, Darby TM, Luo L, Wolfarth AA, Scharer CD, et al. Lactobacilli modulate epithelial cytoprotection through the Nrf2 pathway. Cell Rep. 2015;12:1217–25.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
17.Cheng D, Guo Z, Riegler M, Xi Z, Liang G, Xu Y. Gut symbiont enhances insecticide resistance in a significant pest, the oriental fruit fly Bactrocera dorsalis (Hendel). Microbiome. 2017;5:13.PubMed 
PubMed Central 
Article 

Google Scholar 
18.Pang R, Chen M, Yue L, Xing K, Li T, Kang K, et al. A distinct strain of Arsenophonus symbiont decreases insecticide resistance in its insect host. PLoS Genet. 2018;14:e1007725.PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
19.Kikuchi Y, Tada A, Musolin DL, Hari N, Hosokawa T, Fujisaki K, et al. Collapse of insect gut symbiosis under simulated climate change. mBio. 2016;7:e01578–16.PubMed 
PubMed Central 
Article 

Google Scholar 
20.Corbin C, Heyworth ER, Ferrari J, Hurst GDD. Heritable symbionts in a world of varying temperature. Heredity. 2017;118:10–20.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
21.Jia FX, Yang MS, Yang WJ, Wang JJ. Influence of continuous high temperature conditions on Wolbachia infection frequency and the fitness of Liposcelis tricolor (Psocoptera: Liposcelididae). Environ Entomol. 2009;38:1365–72.PubMed 
Article 
PubMed Central 

Google Scholar 
22.Burke G, Fiehn O, Moran N. Effects of facultative symbionts and heat stress on the metabolome of pea aphids. ISME J. 2010;4:242–52.PubMed 
Article 
PubMed Central 

Google Scholar 
23.Fan Y, Wernegreen JJ. Can’t take the heat: high temperature depletes bacterial endosymbionts of ants. Micro Ecol. 2013;66:727–33.Article 

Google Scholar 
24.Hussain M, Akutse KS, Ravindran K, Lin Y, Bamisile BS, Qasim M, et al. Effects of different temperature regimes on survival of Diaphorina citri and its endosymbiotic bacterial communities. Environ Microbiol. 2017;19:3439–49.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
25.Engl T, Eberl N, Gorse C, Krüger T, Schmidt THP, Plarre R, et al. Ancient symbiosis confers desiccation resistance to stored grain pest beetles. Mol Ecol. 2018;27:2095–108.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
26.Zhang XJ, Yu XP, Chen JM. High Temperature effects on yeast-like endosymbiotes and pesticide resistance of the small brown planthopper, Laodelphax striatellus. Rice Sci. 2008;15:326–30.CAS 
Article 

Google Scholar 
27.Zhang B, Zuo TQ, Li HG, Sun LJ, Wang SF, Zhang CY, et al. Effect of heat shock on the susceptibility of Frankliniella occidentalis (Thysanoptera: Thripidae) to insecticides. J Integr Agric. 2016;15:2309–18.CAS 
Article 

Google Scholar 
28.Karimzadeh R, Javanshir M, Hejazi MJ. Individual and combined effects of insecticides, inert dusts and high temperatures on Callosobruchus maculatus (Coleoptera: Chrysomelidae). J Stored Prod Res. 2020;89:10693.Article 

Google Scholar 
29.Michigan State University. Arthropod Pesticide Resistance Database (APRD). East Lansing: Michigan State University; 2020. http://www.pesticideresistance.com/.30.Ju JF, Bing XL, Zhao DS, Guo Y, Xi Z, Hoffmann AA, et al. Wolbachia supplement biotin and riboflavin to enhance reproduction in planthoppers. ISME J. 2019;14:676–87.PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
31.Zhang Y, Tang T, Li W, Cai T, Li J, Wan H. Functional profiling of the gut microbiomes in two different populations of the brown planthopper. Nilaparvata lugens J Asia Pac Entomol. 2018;21:1309–14.Article 

Google Scholar 
32.Ye YH, Seleznev A, Flores HA, Woolfit M, McGraw EA. Gut microbiota in Drosophila melanogaster interacts with Wolbachia but does not contribute to Wolbachia-mediated antiviral protection. J Invertebr Pathol. 2017;143:18–25.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
33.Yamada R, Floate KD, Riegler M, O’Neill SL. Male development time influences the strength of Wolbachia-induced cytoplasmic incompatibility expression in Drosophila melanogaster. Genetics. 2007;177:801–8.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
34.Wari D, Kabir MA, Mujiono K, Hojo Y, Shinya T, Tani A, et al. Honeydew-associated microbes elicit defense responses against brown planthopper in rice. J Exp Bot. 2019;70:1683–96.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
35.Miller ALE, Tindall K, Leonard BR. Bioassays for monitoring insecticide resistance. J Vis Exp. 2010;46:2129.
Google Scholar 
36.Zhang J, Zhang Y, Wang Y, Yang Y, Cang X, Liu Z. Expression induction of P450 genes by imidacloprid in Nilaparvata lugens: a genome-scale analysis. Pestic Biochem Physiol. 2016;132:59–64.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
37.Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2-ΔΔCT method. Methods. 2001;25:402–8.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
38.Noda H, Koizumi Y, Zhang Q, Deng K. Infection density of Wolbachia and incompatibility level in two planthopper species, Laodelphax striatellus and Sogatella furcifera. Insect Biochem Mol Biol. 2001;31:727–37.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
39.Martin M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet J. 2011. https://doi.org/10.14806/ej.17.1.20040.Callahan BJ, McMurdie PJ, Rosen MJ, Han AW, Johnson AJA, Holmes SP. DADA2: high-resolution sample inference from illumina amplicon data. Nat Methods. 2016;13:581–3.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
41.Katoh K, Misawa K, Kuma KI, Miyata T. MAFFT: a novel method for rapid multiple sequence alignment based on fast Fourier transform. Nucleic Acids Res. 2002;30:3059–66.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
42.Bokulich NA, Kaehler BD, Rideout JR, Dillon M, Bolyen E, Knight R, et al. Optimizing taxonomic classification of marker-gene amplicon sequences with QIIME 2’s q2-feature-classifier plugin. Microbiome. 2018;6:90.PubMed 
PubMed Central 
Article 

Google Scholar 
43.Liu S, Ding Z, Zhang C, Yang B, Liu Z. Gene knockdown by intro-thoracic injection of double-stranded RNA in the brown planthopper, Nilaparvata lugens. Insect Biochem Mol Biol. 2010;40:666–71.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
44.Tai V, James ER, Nalep CA, Scheffrahn RH, Perlman SJ, Keelinga PJ. The role of host phylogeny varies in shaping microbial diversity in the hindguts of lower termites. Appl Environ Microbiol. 2015;81:1059–70.PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
45.Bale JS, Hayward SAL. Insect overwintering in a changing climate. J Exp Biol. 2010;213:980–94.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
46.Rahmstorf S, Cazenave A, Church JA, Hansen JE, Keeling RF, Parker DE, et al. Recent climate observations compared to projections. Science. 2007;316:709.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
47.Radchuk V, Reed T, Teplitsky C, van de Pol M, Charmantier A, Hassall C, et al. Adaptive responses of animals to climate change are most likely insufficient. Nat Commun. 2019;10:3019.Article 
CAS 

Google Scholar 
48.Iwamura T, Guzman-Holst A, Murray KA. Accelerating invasion potential of disease vector Aedes aegypti under climate change. Nat Commun. 2020;11:2130.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
49.Li J, Mao T, Wang H, Lu Z, Qu J, Fang Y, et al. The CncC/keap1 pathway is activated in high temperature-induced metamorphosis and mediates the expression of Cyp450 genes in silkworm, Bombyx mori. Biochem Biophys Res Commun. 2019;541:1045–50.Article 
CAS 

Google Scholar 
50.Kalsi M, Palli SR. Transcription factor cap n collar C regulates multiple cytochrome P450 genes conferring adaptation to potato plant allelochemicals and resistance to imidacloprid in Leptinotarsa decemlineata (Say). Insect Biochem Mol Biol. 2017;83:1–12.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
51.Kalsi M, Palli SR. Transcription factors, CncC and Maf, regulate expression of CYP6BQ genes responsible for deltamethrin resistance in Tribolium castaneum. Insect Biochem Mol Biol. 2015;65:47–56.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
52.Misra JR, Lam G, Thummel CS. Constitutive activation of the Nrf2/Keap1 pathway in insecticide-resistant strains of Drosophila. Insect Biochem Mol Biol. 2013;43:1116–24.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
53.Tang B, Cheng Y, Li Y, Li W, Ma Y, Zhou Q, et al. Adipokinetic hormone regulates cytochrome P450-mediated imidacloprid resistance in the brown planthopper, Nilaparvata lugens. Chemosphere. 2020;259:127490.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
54.Cheng Y, Li Y, Li W, Song Y, Zeng R, Lu K. Inhibition of hepatocyte nuclear factor 4 confers imidacloprid resistance in Nilaparvata lugens via the activation of cytochrome P450 and UDP-glycosyltransferase genes. Chemosphere. 2021;263:128269.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
55.Li Y, Liu X, Wang N, Zhang Y, Hoffmann AA, Guo H. Background-dependent Wolbachia-mediated insecticide resistance in Laodelphax striatellus. Environ Microbiol. 2020;22:2653–63.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
56.Berticat C, Rousset F, Raymond M, Berthomieu A, Weill M. High Wolbachia density in insecticide-resistant mosquitoes. Proc R Soc B Biol Sci. 2002;269:1413–6.Article 

Google Scholar 
57.Zhang G, Hussain M, O’Neill SL, Asgari S. Wolbachia uses a host microRNA to regulate transcripts of a methyltransferase, contributing to dengue virus inhibition in Aedes aegypti. Proc Natl Acad Sci USA. 2013;110:10276–81.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
58.Bi J, Sehgal A, Williams JA, Wang YF. Wolbachia affects sleep behavior in Drosophila melanogaster. J Insect Physiol. 2018;107:81–88.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
59.Roughgarden J, Gilbert SF, Rosenberg E, Zilber-Rosenberg I, Lloyd EA. Holobionts as units of selection and a model of their population dynamics and evolution. Biol Theory. 2018;13:44–65.Article 

Google Scholar 
60.Pan X, Zhou G, Wu J, Bian G, Lu P, Raikhel AS, et al. Wolbachia induces reactive oxygen species (ROS)-dependent activation of the Toll pathway to control dengue virus in the mosquito Aedes aegypti. Proc Natl Acad Sci USA. 2012;109:E23–31.PubMed 
Article 
PubMed Central 

Google Scholar 
61.Gong JT, Li Y, Li TP, Liang Y, Hu L, Zhang D, et al. Stable introduction of plant-virus-inhibiting Wolbachia into planthoppers for rice protection. Curr Biol. 2020;30:4837–45.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
62.Elzaki MEA, Li ZF, Wang J, Xu L, Liu N, Zeng RS, et al. Activiation of the nitric oxide cycle by citrulline and arginine restores susceptibility of resistant brown planthoppers to the insecticide imidacloprid. J Hazard Mater. 2020;396:122755.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
63.Werren JH. Biology of Wolbachia. Annu Rev Entomol. 1997;42:587–609.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
64.Kokou F, Sasson G, Nitzan T, Doron-Faigenboim A, Harpaz S, Cnaani A, et al. Host genetic selection for cold tolerance shapes microbiome composition and modulates its response to temperature. Elife. 2018;77:e36398.Article 

Google Scholar  More