Ecology

Proteobacteria are major contributors to marine microbial phosphonate cyclingDatabases for all putative sequences of genes for phosphonate production (pepM, aepY, phpC, mpnS, hepD), substrate-specific catabolism (phnAWXYZ, palA), and broad-specificity catabolism (phnIJM) were created using available public genomes from JGI IMG/MER and GORG-Tropics. Gene identity was verified by the presence of catalytically essential residues (Supplementary Table S2). Phosphonate genes were identified in 10,337 genomes of bacteria and archaea spanning over 100 unique classes, suggesting a wide variety of microorganisms mediate phosphonate production and catabolism (Fig. 2, Supplementary Dataset S1). A high proportion of all collected sequences affiliated with Proteobacteria (Gamma, Alpha, and Beta classes), averaging 52% of the production genes, 78% of substrate-specific catabolism genes, and 88% of broad-specificity catabolism genes before dereplication (Fig. 2).Fig. 2: Phosphonate gene and genome count with taxonomic distribution.Number of sequences and genomes collected for study (A, D, G) with distribution of class-level taxa for all redundant sequences (B, E, H) and marine redundant (C, F, I) sequences. Results are shown for selected genes representing phosphonate (A–C) production, (D–F) substrate-specific catabolism, and (G–I) broad-specificity catabolism. The taxa shown are the 15 classes with the highest representation across all databases.Full size imageOf the 10,337 genomes, 1556 (15%) were confirmed to be marine organisms from 35 different classes (Fig. 2, Supplementary Dataset S1). Proteobacteria had even greater representation in the subset of marine genomes, averaging 65% of marine production genes, 88% of marine substrate-specific catabolism genes, and 96% of marine broad-specificity catabolism genes from the redundant databases (Fig. 2). The dominance of Alphaproteobacteria in the marine subset may be attributed to the wide variety of Pelagibacterales bacterium captured in the database, making up 426 (27%) of the 1556 genomes involved in all three categories of phosphonate cycling. Rhodobacterales (Ascidiaceihabitans sp., Roseovarius sp., Sulfitobacter sp., Labrenzia sp., and Phaeobacter sp.) alongside Rhodospirillales (Thalassobaculum sp., Thalassospira sp., Roseospira sp., Varunaivibrio sp., and Oceanibaculum sp.) were also highly represented among the marine subset with 214 (14%) and 251 (16%) genomes, respectively (Supplemental Dataset S1), though these taxa primarily show potential for phosphonate catabolism rather than production. Vibrionales were well represented in the JGI IMG/MER marine genome subset with 107 (7%) genomes spanning 59 different species including Vibrio lentus, Vibrio breoganii, and Vibrio splendidus.Diverse taxa encode the capacity to produce phosphonate derivativesPhosphonate production is widespread and distributed throughout many different bacteria and archaea. Genes responsible for the first two steps in phosphonate production, pepM and aepY, had the broadest taxonomic distribution within the redundant databases (Shannon indices of 2.66 and 2.76) for all genes in this study, distributed with 0.59 and 0.61 evenness from 70 and 72 unique, verified classes, respectively. Their broad distribution further highlights the ubiquity and necessity of phosphonate compounds to microbial life and function across all environments. Within the marine setting, both pepM and aepY have reference sequences from 22 unique, verified classes which is the second highest class representation in the marine genome subset (Fig. 2). The marine subset of pepM and aepY also have the highest Shannon indices (1.76 and 1.92) distributed with 0.53 and 0.58 evenness, respectively. A majority (87%) of the Alphaproteobacteria phosphonate producers are Pelagibacterales bacterium with other notable taxa including Bacteria: Candidatus Actinomarinaceae, Prochlorococcus sp., Synechococcus sp., Nitrosococcus sp., and MG-I Archaea: Candidatus Nitrosomarinus catalina, Nitrosopumilus maritimus, alongside other unidentified Crenarchaeota and Thaumarchaeota genomes.The gene phpC was found in less than half the number of genomes than pepM and aepY, and encoded by fewer classes in both the general database (47) and marine subset (10). In the full databases, the distribution of retrieved phpC sequences are similar to pepM and aepY with respect to taxonomic ranking, Shannon index (2.49), and evenness (0.60) (Fig. 2A–C, Supplementary Table 5). Within the marine subset, phpC has less Shannon index (1.61) but greater evenness (0.61) than the marine subset of pepM and aepY. All three upstream phosphonate production genes (pepM, aepY, phpC) are found together within Pelagibacterales bacterium, Prochlorococcus sp., Thaumarchaeota, and Crenarchaoeta alongside other taxa such as Oceanospirillales sp., Arenimonas donghaensis, Desulfuromusa kysingii, and Cellulosilyticum lentocellum.We further investigated the relationship between pepM, aepY, and phpC by examining co-occurrence in genomes and synteny with the general, redundant databases. The first two steps in phosphonate biosynthesis are intimately linked (Fig. 3). Out of all genomes with pepM, 86% have aepY, and out of all genomes with aepY, 90% have pepM. By contrast, phpC is not as closely tied to pepM and phosphonate production. We found phpC in just over 20% of genomes with the capability of phosphonate production (Fig. 3), implying that a majority of bacterial and archaeal phosphonate production stops at the production of phosphonoacetaldehyde or 2-AEP (Fig. 1A). Furthermore, half of the phpC genes were not associated with phosphonate production, given 53% of genomes with phpC did not have pepM and 54% did not have aepY (Fig. 3). In these instances, microbes may use phpC within a 2-AEP substrate-specific catabolism operon (Fig. 3) that allows phosphonate compounds to be synthesized by transforming 2-AEP with phnW and phpC into 2-HEP (Figs. 1A and 3). By repurposing 2-AEP, individuals can still create the specific compound needed while bypassing the energetically unfavourable first step of phosphonate production.Fig. 3: Co-occurrence of phosphonate cycling genes within the same genome and examples of genetic organization of phosphonate cycling genes.The heatmap displays co-occurrence of phosphonate cycling genes. Each column represents the subset of all genomes which contain the source gene and the heatmap value represents the fraction of the source genomes which also contain the co-occurring gene. Heatmap values are not symmetrical due to differing number of genomes represented in each column, database size listed above each column. Examples for phosphonate cycling genomic neighbourhoods were chosen to maximize diversity in synteny with examples from both Bacteria and Archaea where applicable. Several phosphonate-specific ABC transport system clusters are labelled as follows: phnC = phosphonate transport system ATPase; phnD = phosphonate transport system substrate-binding; phnE = phosphonate transport system permease; phnS = 2-AEP transport system substrate-binding; phnT = 2-AEP transport system ATP-binding; phnV = 2-AEP transport system permease; palC = transport system permease; palD = transport system ATP-binding; palE = transport system permease. Genes are colour coded by: red = lyase; orange = transcriptional regulator; yellow = hydrolase; green = transferase; light blue = oxidoreductase; dark blue = transaminase; purple = kinase; pink = isomerase; brown = transport; white = synthase; black = uncharacterized protein; grey = unknown.Full size imageA narrow but diverse selection of taxa encoded MpnS, the marker gene for Mpn production and a key determinant in marine methane production. We observed distinct clades of this enzyme in autotrophic archaea and heterotrophic bacteria (Fig. 2B, C). Within the marine ecosystems, Pelagibacterales, Rhodospirillales, Rickettsiales, Oceanospirillales, Flavobacteriales, and Synechococcales are bacterial candidates for MPn production alongside Thaumarchaeota and Crenarchaeota archaeon (Fig. 2B, C). While six of the bacterial genomes with MpnS also encoded genes for phosphonate catabolism, none of the archaeal MPn producers showed capacity for catabolism (Supplementary Dataset S1). The genomic neighbourhoods for general phosphonate production (pepM, aepY, phpC) and MPn production (mpnS) in both bacteria and archaea include genes such as glycosyltransferase, lipopolysaccharide choline phosphotransferase, choline kinase, adenylyltransferase, and arylsulfatase A (Fig. 3) suggesting the potential for synthesis of (methyl)phosphonate esters [93]. This is consistent with previous analysis [29] of the Nitrosopumilus maritimus SCM1 MPn production genomic neighbourhood and biophysical evidence that MPn producing archaea synthesize an exopolysaccharide modified with MPn similar to 2-AEP modified polymers.Contrary to the diversity of the other phosphonate production databases, the hepD database has low Shannon index (0.62) and evenness (0.45) with 79% of sequences mapping to Actinomycetia including Streptomycetales and Corynebacteriales (Fig. 2B, C). The marine subset has lower Shannon index (0.28) and evenness (0.41) where all sequences derive from Pelagibacterales except one from Prochlorococcus sp. The genomic neighbourhood of HMP production may contain genes for cell surface modification such as acetyltransferase, peptidoglycan biosynthesis, and adenylylsulfate transferase, suggesting that some organisms may use HMP as a conjugate for membrane-associated or exported macromolecules similar to theories on MPn utilization. Other examples of hepD synteny contain more specific genes such as the HMP dehydrogenase or other enzymes for downstream modification (Fig. 3).Marine proteobacteria encode genes for substrate-specific and broad-specificity phosphonate catabolismGenes for marine substrate-specific phosphonate catabolism were widespread among Proteobacterial classes, and to a lesser extent amongst other classes including Bacilli, Planctomycetes, and Synechococcus (Fig. 2E,F). Marine substrate-specific catabolism has lower average Shannon index (1.00) and evenness (0.43) than the three general production genes (pepM, aepY, phpC). The most widespread of these genes was phnW, likely due to its pivotal role in 2-AEP transformations as a precursor reaction to phnAY or phnX (Fig. 1, Supplementary Table 5). Marine hydrolases for 2-AEP catabolism, phnA, phnX, and phnZ, have similar Shannon indices (mean: 1.11 ± 0.05) and evenness (mean: 0.41 ± 0.03) (Fig. 2E, F, Supplementary Table 5).While not exclusive, sequenced references demonstrate a strong taxonomic partition between Proteobacterial classes for 2-AEP catabolism pathways phnAWY and phnWX. Over 74% of marine genomes with phnAWY are Alphaproteobacteria, in particular Rhodobacterales species such as Roseovarius nubinhibens, Marivita geojedonensis, and Pelagicola litoralis. On the contrary, ~80% of marine genomes with phnWX are Gammaproteobacteria, specifically of Vibrionales, Oceanospirillales, and Alteromonadales including a wide range of species from Vibrio, Photobacterium, Marinobacterium, Halomonas, and Pseudoalteromonas.Taxonomic distribution for marine phnZ was 72% Alphaproteobacteria with Pelagibacterales making up 45% of marine phnZ sequences. Note that phnZ has the most (17) reference sequences from marine Cyanobacteriia, specifically Prochlorococcus sp., than any other phosphonate catabolizing gene. Lack of marine sequence representatives for catabolism of phosphonopyruvate by palA suggests that either the substrate is uncommon, therefore the function unnecessary, or marine microbes have other methods of catabolizing phosphonopyruvate, perhaps by the C-P lyase. Overall taxonomic distribution of phosphonate substrate-specific catabolism, specifically targeting 2-AEP, suggests said function is essential to many marine heterotrophs within Alphaproteobacteria and Gammaproteobacteria. However 2-AEP catabolism appears to be less universally important than phosphonate production to marine microbial life since the required genes are found in a less diverse selection of taxa.Genetic organization for substrate-specific catabolism genes, particularly those targeting 2-AEP, varied widely in line with the numerous options for 2-AEP catabolism (Fig. 3). Though some bacteria specialize in a single 2-AEP degradation pathway such as only containing phnWAY, others contained multiple hydrolases for 2-AEP catabolism with some incorporating phpC into a 2-AEP specific catabolism operon (Fig. 3). When a genome has two hydrolases for phosphonate catabolism, often phnZ was paired with either phnA or phnX. Co-occurrence between phnZ and either phnA or phnX ranged between 30-50%, whereas co-occurrence between phnA and phnX was between 6-12% (Fig. 3). This discrepancy in co-occurrence may be due to the metabolic similarity between phnA and phnX, where having both may be redundant. Both of these enzymes rely on phnW for 2-AEP catabolism and produce carbon metabolites, whereas phnZ does not need phnW and produces the amino acid glycine (Fig. 1B).C-P lyase genes representing substrate non-specific catabolism were overwhelmingly attributed to Alphaproteobacteria which consisted over 75% of all collected marine sequences for phnIJM (Fig. 2H, I). A wide variety of Rhodobacterales, spanning 55 different genus are the most numerous representatives, followed by Pelagibacterales and Rhodospirillales. The genes in all three databases have very high genome co-occurrence, 89–99%, as expected given all three operate within the same enzyme complex (Fig. 3). Gene co-occurrence, Shannon index, and evenness is lower for phnM than the other two C-P lyase components, phnI and phnJ, likely due to instances of organisms containing two copies of phnM where one copy lies outside the C-P lyase operon [94]. C-P lyase gene databases have lower Shannon index (mean 0.81 ± 0.11) than phosphonate production and 2-AEP substrate-specific catabolism genes (phnAWXZ) (Fig. 3D, G), suggesting broad-specificity phosphonate catabolism by the C-P lyase is a narrowly distributed function (Supplementary Table 5). Organization of C-P lyase operons held the most consistency between example genomes, likely due to the high number of genes simultaneously utilized for lyase construction. These operons encoded a consecutive string of lyase subunits, including a generic phosphonate transporter (phnCDE) and GntR transcriptional regulator (Fig. 3). C-P lyase genes had low genomic co-occurrence with all other phosphonate cycling genes with notable co-occurrence between phnW at 26%, phnZ at 21%, and phnX and 18% (Fig. 3). The low rate of co-occurrence may be due to redundancy in function for P harvesting between the C-P lyase and substrate-specific catabolism. In some cases there are instances of a substrate-specific hydrolase gene located within the C-P lyase operon (Fig. 3).Phosphonate biosynthesis genes are globally prevalent in oceans and increase in mesopelagic watersFollowing curation of phn-gene databases, we analysed 121 metagenomes and 91 metatranscriptomes from the publicly available TARA Oceans expedition (spanning samples from the Atlantic, Indian, Pacific, and Southern Ocean and Red Sea) to investigate the global potential for marine phosphonate cycling. Measuring the proportion of the community capable of performing specific tasks through metagenomics indicates the long-term selective pressures that shape P-cycling and microbial communities.Potential for phosphonate production (pepM) was globally ubiquitous across all depths, with 14–17% of the community encoding in the surface waters and deep chlorophyll maximum (DCM), increasing to 45% in mesopalagic waters (Fig. 4A, Supplementary Tables S6 and S7), highlighting the importance of phosphonate compounds to marine microbial communities. Relative abundance of phpC was 64–76% that of pepM and aepY across all depths (Fig. 4A). We observed significant increase in relative abundance between the surface and mesopelagic for pepM (ANOVA: F = 1262, p  More

DEGs under low-temperature stressThe results from the field experiments indicated that the daily mean values of A, Fvʹ/Fmʹ, ETR and Fv/Fm decreased in the LT group, the PSII function was impaired, and the photosynthetic capacity was weakened. Through the specific analysis of the “Photosynthesis” pathway (pathway ID ko00195) in the LT group, it was found that PSII, the cytochrome b6f. complex (Cyt b6f.), PSI and ATPase exhibited differential gene expressions. Figure 9 shows the structural pattern diagram for photosynthesis. The parts marked by white boxes indicate that the structure has DEGs. The gene expressions of CP43, CP47, D1 protein and Cytb559 of PSII changed. The inner peripheral antenna pigment proteins, CP43 and CP47, of PSII bind to chlorophyll. They accept the excitation energy transferred from the surrounding antenna complex and transfer this energy to the reaction centre complex. Changes in CP43 and CP47 affect the absorption and transmission of light energy. In the PSII reaction centre, light energy is converted into chemical energy. P680 absorbs light and is excited to become P680*, and then transfers electrons to pheophytin (Pheo). At the same time, the PSII oxygen-evolving complex obtains electrons from water molecules, the water molecules are split and releases oxygen and protons. As one of the two core proteins that compose the reaction centre complex, the D1 protein combines with various cofactors that are related to the original charge separation and electron transfer. The D1 protein plays an important role in the process of photosynthetic electron transfer. Studies have found that low temperatures can induce allosteric inactivation of the D1 protein, which results in changes in the structure of thylakoid membranes and hinders electron transfer8. As part of the reaction centre, Cytb559 can adjust the photoinhibition sensitivity of PSII through redox changes so that the PSII reaction centre is protected from damage9. The light energy absorption, energy conversion and electron transfer functions of PSII are impaired, which result in significant decreases in Fv/Fm to levels far below the normal value. The results of Xiangchun Song are similar to those presented in this paper: the PS II reaction centre of A. mongolicus seedlings is irreversibly inactivated or the thylakoid membrane is damaged under subzero low temperature stress, which may produce serious photoinhibition. However, Song believes that the peripheral antenna component of the optical system is more affected than the core complex at low temperatures, which was not observed in the corresponding results in this study10.Figure 9Photosynthesis of A. mongolicus under low-temperature stress. The areas outlined by white boxes indicate the differentially expressed genes in these structures.Full size imageThe gene expressions of Cyt b6, PrtD and Cyt f in Cyt b6f. changed. Cyt b6f. changes not only affect the electron transport function of photosynthesis but also affect ATP synthesis. Pheo transfers the received electrons to plastid quinone (PQ). PQ receives electrons and protons to form plastid hydroquinone (PQH2). Then, the electrons of PQH2 are transferred to plastid cyanin (PC) on PSI through Cyt b6f., and hydrogen protons are released into the cavity of the thylakoid to form a transmembrane proton gradient. The transmembrane proton gradient is the driving force for ATP synthesis.The function of PSI is to transfer electrons from PC to ferredoxin for the reduction of NADP+. Recent studies have found that PSI is more sensitive to light and more prone to selective photoinhibition than PS II under low temperature and weak light conditions11,12. The KEGG analysis results indicated that the LHCI complex, PsaF and PsaE subunits of PSI showed differential gene expressions. The main function of the LHCI light-harvesting pigment protein complex is to capture light energy. PsaF is a low-molecular-weight protein that is distributed in the membrane. Some studies have suggested that the N-terminal amino acid sequence of eukaryotic PsaF is involved in the binding of PSI and PC13. PsaE, PsaD and PsaC together form the docking site of ferredoxin on the PSI receptor side14,15. Ferredoxin and ferredoxin-NADP+ reductase in the photosynthetic electron transport chain are also affected, which results in hindrance of NADPH synthesis. The F-type H+/Na+ transport ATPase subunits also show differential gene expressions, which lead to impaired ATP synthesis. Low temperatures affect the ability to absorb light energy, transfer electrons, convert light energy into electric energy, and synthesize NADPH as well as ATP, which ultimately lead to declines in Fv’/Fm’ and ETR and impair the photosynthesis capacity of A. mongolicus.Compared with the light reaction, low temperatures have a greater impact on the dark reaction. Because the dark reaction process is composed of many complex enzymatic reactions, the enzyme activity is very susceptible to temperature. The KEGG results show that 13 related enzymes were differentially expressed in the “carbon sequestration of photosynthesis” (ko00710). The Rubisco enzyme is a key enzyme that determines the direction and efficiency of photosynthetic carbon metabolism in C3 plants and is sensitive to temperature16. The results also show that the expression levels of 10 differentially expressed genes of Rubisco enzymes all declined. In the Calvin cycle, the gene expressions of only transketolase and glyceraldehyde-3-phosphate dehydrogenase are not sensitive to temperature. In addition, the reduction phase of the dark reaction requires the use of NADPH and ATP that are produced by the light reaction. The inhibition of NADPH and ATP synthesis will inevitably affect the normal progression of the Calvin cycle.Chloroplast respiration is an O2-dependent electron transport pathway in chloroplasts. Chloroplast respiration includes the nonphotochemical reduction of PQ by NAD(P) H and the reoxidation of PQ by terminal oxidase, which can consume excess electrons to protect plants from damage due to photooxidation.Figure 10 shows the partial KEGG enrichment metabolic pathway in the LT group. There were three significant enrichment pathways related to carbohydrate metabolism: fructose and mannose metabolism (ko00051), butanoate metabolism (ko00650) and C5-branched dibasic acid metabolism (ko00660). The metabolism of fructose and mannose includes the ascorbic acid biosynthetic pathway. Ascorbic acid (ASA), also known as vitamin C, can be used as a cofactor of violaxanthin de-epoxidase to participate in the lutein cycle and consume excess light energy and protect plants from harm.Figure 10The regulatory mechanism of A. mongolicus under low-temperature stress. The white ovals represent the enriched metabolic pathways. The blue rectangles represent significantly enriched KEGG metabolic pathways. The pathways are followed by the physiological structures and substances or physiological processes in which the expressions of related genes change.Full size imageLow temperatures damage cell membranes first. Increasing the mass fraction of unsaturated fatty acids in the membrane is beneficial to improve the stability and fluidity of the membrane. Some studies have shown that the degree of unsaturation of fatty acids in adult leaves of A. mongolicus that grow naturally in the field is lower in summer and higher in autumn and winter17. The significantly enriched pathways related to unsaturated fatty acid metabolism were alpha-linolenic acid metabolism (ko00592), linoleic acid metabolism (ko00591) and arachidonic acid metabolism (ko00590). Various proteins, such as linoleate 13S-lipoxygenase and cytochrome P450 family 2 subfamily J (CYP2J), which are involved in the metabolism of linoleic acid, showed differences in their gene expressions. Linoleate 13S-lipoxygenase is a common lipoxygenase in plants that can catalyse the production of precursors of several important compounds, including jasmonic acid. CYP2J is a group of P450 haem thiolate proteins, which are mainly distributed on the endoplasmic reticulum and inner mitochondrial membrane and are involved in the synthesis of sterol hormones, including brassinosteroids. Because light systems are distributed on the thylakoid membrane, damage to this membrane will affect the progress of plant photosynthesis.Plant hormone signal transduction (ko04075) plays an important role in plant resistance to stress. Studies have shown that JAs have physiological functions, such as inducing stomatal closure, inhibiting photosynthesis, promoting respiration and promoting leaf senescence18,19. Treating plants with exogenous methyl jasmonate can induce the transcription of the heat shock protein family, increase the synthesis of antioxidants, reduce lipoxygenase activity and enhance the ability of plants to resist cold damage20.Figure 11 shows the regulatory mechanism of A. mongolicus in the HL group. The MapMan analysis results show that the DEGs of the LHCII complex and those for the assembly and maintenance of PSII are significantly changed. LHCII contains chlorophyll and carotenoids, which can capture and transmit light energy. Chlorophyll is an important photosynthetic pigment that captures light energy and drives electrons to the reaction centre. The chlorophyll molecule in the reaction centre is related to photochemical quenching. The entire chlorophyll biosynthesis process (e.g., L-glutamyl-tRNA → chlorophyll a → chlorophyll b) involves 15 enzymes. The analysis found that 4/5 of the enzymes’ expression genes were changed. Carotenoids include carotene and lutein, and their synthesis is affected by high temperatures. Lutein participates in the lutein cycle, which can dissipate excess light energy and prevent membrane lipids from being peroxidized and thus maintain the stability of the thylakoid membrane structure and protect A. mongolicus. from high temperature stress and strong light stress.Figure 11The regulatory mechanism of A. mongolicus. under high-temperature stress. The white ovals represent enriched metabolic pathways. The red rectangles represent significantly enriched KEGG metabolic pathways. The pathways are followed by the physiological structures and substances or physiological processes in which the expressions of related genes change.Full size imageThe D1 protein in the PSII reaction centre is rapidly degraded under strong light conditions. To maintain the normal physiological needs of plants, the degraded D1 protein will be replaced by the new D1 protein that is produced by the repair mechanism. The reversible inactivation of the PSII reaction centre can protect the photosynthetic system and avoid destruction. This may be the reason for the significant changes in the DEGs that are involved in the assembly and maintenance of PSII.Rubisco is the main site for high-temperature inhibition of the Calvin cycle16. The KEGG analysis found that there were 7 (4↑, 3↓) DEGs of Rubisco. SBPase catalyses the conversion of sedum heptulose-1,7-diphosphate (SBP) into sedum heptulose-7-phosphate (S7P) in the renewal phase. Under low-temperature stress, only transketolase and glyceraldehyde-3-phosphate dehydrogenase remained unchanged in the Calvin cycle. In addition, NDH-mediated cyclic electron transfer may decreased the photooxidation damage that is caused by high-temperature stress by shunting the excess electrons that were generated by the inhibition of CO2 assimilation to the chloroplast respiratory pathway21.In the HT group, the net photosynthetic rates of the leaves showed two peaks on the diurnal change curves, and there was an obvious phenomenon of midday photosynthesis depression. The daily average A values were greater than those of the CK group. These results show that A. mongolicus has a complete photosynthetic structure protection mechanism and can adapt to high-temperature environments. The pathway of significant enrichment related to carbohydrate metabolism in the HT group was the same as that in the LT group. The enrichment degrees of the fructose and mannose metabolic pathways were higher only in the HT group, and C5-branched dibasic acid metabolism and butanoate metabolism were higher in the LT group.Under high temperature and strong light conditions, the balance between production and removal of reactive oxygen species (ROS) in plant cells was broken, and large amounts of reactive oxygen species accumulated in the cells. Active oxygen can cause lipid peroxidation of the biomembrane, enlarge membrane pores, increase the permeability, and affect the spatial structures of enzymes on the membrane, which thus leads to chloroplast destruction. In severe cases, ROS will cause serious injury or even death to plants22. The gene expressions of FabH and acetyl-CoA carboxylase (ACCase) changed during the synthesis of unsaturated fatty acids in the HT group.There are two types of active oxygen scavenging mechanisms in plants. (1) The enzymatic detoxification system: superoxide dismutase (SOD), ascorbate peroxidase (APX), and catalase (CAT). (2) Nonenzymatic antioxidants: ASA, carotenoids, glutathione, mannitol, and flavonoids23.Secondary metabolites result from long-term adaptation of plants to their environments. They can improve the ability of plants to protect themselves, compete for survival, and coordinate the relationship between plants and the environment. The significant enrichment pathways related to the biosynthesis of secondary metabolites in the HT group consisted of phenylpropane biosynthesis (ko00940), flavonoid biosynthesis (ko00941) and isoflavone biosynthesis (ko00943). The phenylpropanoid biosynthesis pathway is one of the three main secondary metabolic pathways in plants. It starts from phenylalanine and generates different phenylpropane metabolites through multistep reactions, such as flavonoids, isoflavones, anthocyanins and lignin24,25. Anthocyanins can protect plants from light damage by quenching free oxygen radicals and reducing the absorption of light energy. Hughes studied 10 species of evergreen broad-leaved trees and found that red leaves containing anthocyanins always maintained higher Fv/Fm levels than green leaves. Fv’/Fm’ is related to nonphotochemical quenching. This means that trees with red leaves rely more on the light-damage defence function of anthocyanins than on the light-damage defence mediated by lutein26.Riboflavin metabolism (ko00740) and biotin metabolism (ko00780) are two significantly enriched cofactors and vitamin metabolic pathways. Riboflavin is the precursor of flavin mononucleotide (FMN) and flavin adenine dinucleotide (FAD). As a prosthetic group of flavinases, FAD participates in multiple biochemical processes, such as mitochondrial electron transport, photosynthesis, fatty acid oxidation and folate metabolism, in plants27. Riboflavin can induce antioxidant accumulations in plant cells and can also promote plant growth by affecting the ethylene signalling pathway28. Biotin (e.g., VH or VB7), as an essential cofactor for biotin-dependent carboxylase, plays an important role in the life activities of plants. Common biotin-dependent carboxylase enzymes are pyruvate carboxylase (PC) and ACCase. PC is present in the mitochondria and participates in the replenishment mechanism of the tricarboxylic acid cycle. ACCase plays a pivotal role in the feedback regulation of fatty acid synthesis and is the site of action for the feedback regulation of fatty acid synthesis29.The four pathways related to amino acid metabolism showed differences in the HT group. The enrichment degrees of each pathway were as follows: valine, leucine and isoleucine biosynthesis (ko00290)  > biosynthesis of amino acids (ko01230)  > lysine biosynthesis (ko00300)  > glycine, serine and threonine metabolism (ko00260). The branched chain amino acids, valine, leucine and isoleucine and their derivatives, are beneficial to plant growth and plant responses to stress30. As an essential amino acid, lysine metabolism affects many physiological reactions, such as the tricarboxylic acid cycle, abiotic and biotic stress responses, and starch metabolism31. The glycine, serine and threonine metabolic pathways combined with the GO enrichment results showed that the genes related to glycine catabolism and glycine dehydrogenation/decarboxylase activity changed greatly. It is known that when the activity of mitochondrial glycine decarboxylase increases, both photorespiration and photosynthesis will increase32.In terms of hormones, salicylic acid, cytokinin, and abscisic acid (ABA) can improve plant active oxygen scavenging ability. Salicylic acid can decrease the damage to seedlings due to high temperatures by improving the ability of plants to resist oxidative stress and increasing the contents of osmotic adjustment substances in cells33. Salicylic acid also has the function of delaying the degradation of D1 protein and speeding up the recovery of D1 protein when high temperatures are no longer present34. ABA can improve the heat tolerance of plants by regulating the expressions of heat stress-induced genes at the transcriptional level35.In conclusion, A. mongolicus has weak resistance to low temperatures and good adaptation to high temperatures. At the physiological level, under low-temperature stress, the proportion of Y (NO) increased, the function of PSII was damaged, and photosynthesis was inhibited. A. mongolica maintains normal physiological activities by regulating the circadian rhythm, increasing the synthesis of unsaturated fatty acids and changing the effects of plant hormones. Under high-temperature stress, A. mongolicus maintains normal photosynthesis by adjusting gsw as well as water utilization and by increasing the proportion of Y (NPQ). At the same time, A. mongolicus uses LHCII to consume excess energy, continuously assembles and maintains the normal function of PSII, and changes the types of antioxidants, such as by synthesizing anthocyanins, flavonoids, and isoflavones, to protect itself from injury. In addition, the porphyrin and chlorophyll metabolisms, carotenoid metabolism, plant hormones, amino acid metabolism, unsaturated fatty acid synthesis and other metabolic pathways that are related to the differentially expressed genes changed greatly. More

The avian plumage color gamut is much more diverse than previously estimated2. We demonstrate that hummingbird barbule structural colors contribute substantially to the total color diversity of living birds, occurring in areas of the avian color space that were sparsely occupied in Stoddard and Prum2, which most notably included saturated blues, greens, and true purples (blue + red). Such regions of the avian color space were suggested to be unoccupied because these colors are challenging to create, rather than because they might function poorly for communication2. Our results support this hypothesis because hummingbird coloration densely occupies these regions of the avian color gamut (Fig. 2d), using plumage patches that generally play particularly important roles in hummingbird communication, such as throat and crown plumage patches (Supplementary Fig. 5)16,17. The greater color diversity uncovered by our study suggests that barbule structural coloration is the most versatile class of all plumage coloration mechanisms and poses the least constraints on the evolvability of plumage color diversity. Barbule structural colors evolve through changes in the size, shape, spacing, and refractive index of barbule melanosome nanostructures, but little is known about how changes in these parameters themselves evolve18.The UV/V + green region of avian color space remains mostly unoccupied (Fig. 2c, d). It is challenging to create colors with separate reflectance peaks within the wavelength sensitivities of non-adjacent color cones because the peaks must be highly saturated to avoid stimulating neighboring cones2. However, this idea does not explain why there are far more true purple (blue + red) than UV/V + green plumage colors. Notably, birds particularly fail to fill the more UV/V regions (those closer to the UV/V vertex) of UV/V + green color space, which might indicate that it is difficult to create spectra with uv/v wavelength peaks higher than those in the m wavelengths.The differences between our methods and those of Stoddard and Prum2 likely contribute in part to the larger gamut size when comparing species data but not overall data. While the number of species included in our study was comparable to that of Stoddard and Prum2 (114 vs 111 species, respectively), we measured almost twice as many plumage patches as they did (+1600 vs. 965 patches). To prevent erroneous distortion to iridescent colors we did not average the three measurements per patch. Both studies measured six standard patches for all species and additional patches if necessary to capture other plumage color variation. The larger number of plumage patches we measured reflects how color diverse hummingbird plumages are. Our methods preserved the natural variation in hue due to iridescence and avoided the distorted flattening caused by averaging highly saturated peaks with slightly different peak hues. Although our methods are biased toward increasing variation, they are necessary to accurately capture the phenomenon of iridescent hummingbird coloration.There are multiple reasons why the hummingbird color gamut is so diverse. The size of the hummingbird color gamut, like the achieved color gamut of any clade, constitutes a combination of the history of selection on color function, the clade’s evolved capacities for color production, the age of the clade, and the number of species. Hummingbirds excel at all these criteria. The 336 species of extant hummingbirds have radiated rapidly over the last 22 million years19. Hummingbird plumage color diversity has evolved through a long history of persistent sexual and social selection on plumage coloration. Hummingbirds have polygynous breeding systems characterized by female only parental care, female mate choice, and often elaborate male courtship displays. Intersexual selection in hummingbirds has contributed to elaborate radiation in brilliant plumage coloration as well as vocalizations and non-vocal feather sounds14,16,20. Hummingbird plumage color evolution rates have even been shown to positively correlate with hummingbird speciation rates14. Furthermore, in some species, brilliant monomorphic plumage ornaments apparently function in aggressive, intra- and interspecific defense of floral resources21 and appear to be associated with socioecological features related to resource competition19. Our finding that crown and throat patches, which flash brilliantly when the head of the bird is oriented toward the observer, are more diverse in coloration than other plumage regions highlights the role of plumage coloration in direct inter-individual communication and social interactions.The mechanistic properties of hummingbird barbule structural color further explain the exceptional diversity of hummingbird plumage coloration. Hummingbird barbule structural coloration is among the most complex plumage coloration mechanisms, comprised of stacks of hollow, air-filled melanosomes, surrounded by a thin superficial, solid keratin cortex as well as sometimes superficial, miniature melanin platelets which lie just beneath this cortex9,10,11,12,13. Complex nanostructures allow for independent tuning of multiple components, and, hence, greater achievable color diversity12,18,22. Barbule structural color permits the production of any peak-reflected wavelength by varying the thickness of melanosome arrays, which can produce a diversity of single-peak spectra-hues, such as the unusual diversity of greens, blues, and blue + greens seen in hummingbirds (Fig. 2b). Hummingbird melanosomes are among the most unusual in birds in being both disc-shaped and air-filled9,10,11,12,13,23. The air in the center of hummingbird melanosomes approaches the maximum possible biological difference in refractive index (air = 1.0, melanin = ~1.7), which results in the efficient production of brilliant colors with the fewest layers of melanosomes, such that resulting spectra are narrow and near saturation13,24. Such spectra can thereby create colors that extend further in color space (Fig. 2a–c).Barbule structural color also allows for the production of plumage spectra with multiple saturated peaks, creating saturated color combinations that are not as commonly produced via other plumage coloration mechanisms. However, researchers have yet to identify exactly how hummingbird multipeak spectra are produced12,13, emphasizing the need for further analyses of the optics of hummingbird feathers. Many hummingbird melanosome arrays are non-ideal– i.e., the products of the thicknesses and refractive indices of the melanin and air cavity layers are not equal25. Non-ideal thin films can create more highly saturated, pure tone colors of the primary peak while also introducing additional, harmonic spectral peaks at shorter wavelengths25, which allows for complex reflectance spectra with multiple bright peaks within the avian visible spectrum. Also, melanosome arrays with a large average layer thickness ( >~300 nm) can create colors with fundamental interference peaks in the infrared and multiple, harmonic peaks in the avian visible range (300–700 nm). The presence of minute, superficial melanin platelets below the cortex in hummingbird barbules is also correlated with secondary, lower wavelength reflectance peaks, but the precise optical mechanism remains to be established12. These different nanostructural elements all contribute to distinctive multipeak reflectance spectra that can stimulate non-adjacent color cone combinations, which Stoddard and Prum2 identified as particularly difficult to accomplish: UV/V-purple (uv/v + s + l wavelengths; Schistes geoffroyi cheek, Fig. 4g); true purple (s + l wavelengths; Atthis ellioti gorget, Fig. 4h); UV/V-green (uv/v + m; Schistes geoffroyi crown, Fig. 4a); and UV/V-red (uv/v + l; Heliangelus viola, Fig. 4b). With multipeak spectra the potential for creating new and different colors is greatly expanded, allowing for a more versatile evolution of novel colors.Unexpectedly, the hummingbird plumage color gamut is larger in volume when modeled with the VS-type (34.2%) than with the UVS-type (29.6%) visual system. This apparently unique result contrasts notably with both Stoddard and Prum’s2 and our revised estimate of the color gamut of all birds combined– VS gamut = 40.5%; UVS gamut = 47.3%. Multiple previous analyses have shown that the UVS cone-type visual system does a more efficient job of discriminating the colors of natural objects because of the broader separation between the peak spectral sensitivities of the uv and s (blue) cone types2,26,27. Because the UVS-type visual system produces an even greater increase in color volume for a diverse plant color data set over the VS-type visual system, Stoddard and Prum2 rejected the hypothesis that the UVS-type visual system had specifically evolved to expand the diversity of avian color stimuli.However, our observations that the hummingbird plumage gamut is substantially greater in volume with the VS-visual system than with the more efficient UVS-visual system strongly suggests another hypothesis: Hummingbird plumage may have specifically evolved to be more diverse within the hummingbird VS-type color visual system via selection for highly saturated plumage colors. Given diversity in hue, the way to achieve greater color gamut volume, i.e., greater plumage color diversity, is through highly chromatic color vectors that extend toward the limits of the color space. The two visual systems map variation in wavelength to different maximum potential chroma—i.e., wavelengths with color vectors that extend toward the edges, faces, and vertices of the tetrahedron6. Color vectors that extend towards the vertices, i.e., plumage that best corresponds to a singular cone type’s peak sensitivity, have the highest maximum potential chroma because vertices are the regions furthest away from the tetrahedron’s center. Thus, hummingbird plumages may have specifically evolved to have maximum chroma within their own VS-visual system via peaks that correspond most closely to the peak sensitivities of the VS- rather than the UVS-visual system. For example, when comparing the UVS and VS plumage color gamuts for hummingbirds, it is notable that hummingbird coloration extends much further into the UV/V regions of color space for the VS-visual system (Supplementary Fig. 2). While in the VS system these color points map toward the v vertex, in the UVS-visual system they map towards the uv-s edge and the uv-s-l face. Such color vectors that contribute to expanded color volume of the VS gamut could have evolved by sexual or social selection for highly saturated plumage colors that are near in hue to the specific sensitivity peaks of hummingbird receptor cone types. Such selection could note preferences within some hummingbird species for hues with maximally possible chroma, not merely for maximal chroma of a given hue.Hummingbirds have tetrachromatic color vision with substantial sensitivity in the near ultraviolet28,29. Recently, Stoddard et al.30 used a series of elegant experiments with hummingbird feeders and LED lights to demonstrate for the first time that hummingbirds can distinguish non-spectral colors distributed throughout the tetrachromatic color space. However, the presence of this remarkably proficient four-color vision in hummingbirds poses an interesting evolutionary conundrum. Recent phylogenetic analyses have established that hummingbirds and swifts are phylogenetically embedded within the nocturnal caprimulgiforms31,32. The most parsimonious hypothesis is that the immediate ancestors of swifts and hummingbirds were extensively nocturnal for approximately 8 million years before they re-evolved diurnal ecology and behavior31. Given that an evolutionary history of nocturnality can lead to the degradation or loss of opsin genes33,34, it should be a high priority to establish what effect that ancestral nocturnality may have had on the molecular physiology and anatomy of the hummingbird color visual system.Our attempt to document the color diversity of an avian family has revealed that current estimates of the total avian color gamut are likely inaccurately low. Similar studies sampling from other color-diverse families, such as sunbirds (Nectariniidae), parrots (Psittacidae), tanagers (Thraupidae), birds of paradise (Paradiseidae), manakins (Pipridae), and starlings (Sturnidae), most of which have already been studied for their plumage coloration35,36,37,38,39, would help us obtain a better estimate of the true avian color gamut. More

In their original description of V. rhodanica, Fischer & Riou16 determined that the previously undescribed genus was a Jurassic relative of V. infernalis. This assignment was based on the configuration of the arm crown and armature, fin type, presence of luminous organs, lateral eyes, and the absence of an ink sac. Assuming this assignment is correct, then V. rhodanica is a member of the suborder Vampyromorphina, which includes the family Vampyroteuthidae22,29.Reappraisal of the anatomy shows that V. rhodanica and V. infernalis both have 8 arms and uniserial suckers flanked by cirri. They both possess V. infernalis-like sucker attachments34,36, which are broader at the base and taper up to a radially symmetrical sucker.Both species have distinctly modified arms though the morphology differs in each. V. infernalis, has retractable filaments in the position of arm pair II27,33,34, though there is no evidence of these appendages in V. rhodanica. Instead, the species has elongate dorsal arms (arm pair I) with a unique configuration of suckers and cirri on the distal section.The suckers and cirri of V. rhodanica are more numerous than those of V. infernalis27,37. They are also more closely positioned. Proportionally, the suckers of both species have a consistent ratio to mantle length37, though the diameter of the cirri and infundibulum are greater in V. rhodanica. The V. infernalis-like attachment1,3,34 is present in both species, though in V. rhodanica, the distal part of the neck protrudes into the acetabular cavity. Of note, the sucker stalks on the dorsal arms of V. rhodanica are more elongate than those on the other arms (Figs. 2b,c, and 3a,b). This variation in suckers and their attachments suggests a specialized function between the dorsal and sessile appendages. On the longer dorsal arms, the larger sucker diameter, and more elongate stalks (Figs. 2b and 4) indicate the potential for increased mobility over their extant relatives, and possibly facilitated additional manipulation and prey capture capability.Figure 4Hypothesised reconstruction of V. rhodanica based on the data from this study (A. Lethiers, CR2P). The scale is based on measurements from the holotype (MNHN.B.74247) and the arm crown is completed using dimensions from MNHN.B.74244.Full size imageIn addition to the arm crown specialization, V. rhodanica has a more streamlined shape than V. infernalis, which is caused by a proportionally narrower head. Their muscular body is narrower and more elongate than the gelatinous V. infernalis16,27,37 suggesting a higher energy locomotory style. This is consistent with increased predation relative to the modern form. Observations in this study support many assertions of Fischer & Riou16 about the characters in V. rhodanica, though the presence of luminous organs cannot be confirmed. Rather than luminous organs much larger than those present in the deep-sea, extant V. infernalis, it is possible that these structures represent displaced cartilage prior to fossilization (Supplementary Fig. 6).Two other genera from the La Voulte-sur-Rhône locality, Gramadella and Proteroctopus are, like V. rhodanica, considered to be Incertae sedis Vampyromorpha22. All three share morphological similarities that include an elongated mantle fused with the head, and a longer dorsal arm pair with armature on the distal ends1,16,22,38. Neither the second nor fourth arm pair have been modified. Each has one pair of fins. In Gramadella, the fins are lateral and skirt-like16,38. In V. rhodanica and Proteroctopus these fins are located posteriorly1,16.V. rhodanica shows the greatest length variation between the dorsal and sessile arms (Fig. 4), though proportionally, Gramadella, and Proteroctopus have longer dorsal arms1,31. Fischer & Riou31 and Kruta et al.1 described biserial suckers in their descriptions of Gramadella, and Proteroctopus, respectively. In Proteroctopus, these suckers have a proportionally smaller diameter than the uniserial row in V. rhodanica, and do not exhibit the same tapered pattern.None of these specimens shows evidence of an ink sac, though it is present in contemporaneous genera from the same assemblage (Mastigophora, Rhomboteuthis and Romaniteuthis)8,16. That this character occurs only in some taxa from the same assemblage suggests variation in ecology, possibly associated with the steep, bathymetric relief in the La Voulte-sur-Rhône paleoenvironment11. The mosaic of characters found within the coleoid taxa at La Voulte-sur-Rhône suggests that Mesozoic vampyromorphs co-occurred in different ecological niches during the mid-Jurassic.Today, extant V. infernalis is uniquely adapted to a low-energy, deep-sea mode of life27,28,29,39, though the timing of character acquisition and progression of this ecology is unclear24. It is hypothesised that the vampyromorph Necroteuthis Kretzoi 1942 was already exploiting this niche by the Oligocene29, and that the initial shift to offshore environments was possibly driven by onshore competition24,29. The data obtained here suggests that V. rhodanica, the purportedly oldest-known genus of the Vampyromorphina group, was an active predator following a pelagic mode of life.Indeed, several anatomical details, mainly found in the brachial crown, seem to support this hypothesis. Though we cannot directly compare functionality of the arm crown elements with other Jurassic taxa, we can infer function based on observation in modern forms. In Octopoda, the sister group to Vampyromorpha, suckers are attached to the arm by a cylindrical layer of muscle, encircling oblique musculature40,41, that connects the arm musculature and the lateral margin of the acetabulum34,40,41,42. This facilitates a variety of functions including locomotion, manipulation, and prey retention43. The sucker attaches by flattening the infundibulum against the surface and then the encircling epithelium creates a watertight seal36,40,41,42,43,44,45. Contraction of the radial acetabular muscles provides the pressure differential required to create the suction force43,44,46.The stalked sucker attachments2,34 of decabrachians (Fig. 3d, and Supplementary Fig. 4) are muscular35 and connect the musculature of the arm with the base of the sucker, forming part of the acetabulum33,34. Tension on the sucker stretches this muscular attachment, which pulls locally on the acetabular base. This facilitates a greater pressure differential inside the sucker, allowing the teeth on the sucker ring to maintain the hold47.Extant V. infernalis lack decabrachian-like stalks2,18 and the neck of the attachment joins to the base of the acetabulum (Fig. 3c, and Supplementary Fig. 4), rather than being inserted into it18. The infundibulum is not distinct, and the suckers do not provide strong suction27. Instead, suckers function by secreting mucus to coat detritus—marine snow captured by retractable filaments—which is then moved to the mouth by cirri7,27.A mosaic of these characters is present in V. rhodanica (Fig. 3a,b), therefore, suggesting their potential for increased attachment and hold on prey over extant V. infernalis. These include a larger infundibular diameter, a neck attachment integrated with the acetabular muscles, and the elongated stalks of the dorsal suckers.Additionally, the paired, filamentous cirri observed in extant cirrates48 are present in V. rhodanica (Fig. 4, and Supplementary Fig. 2). In extant forms they are understood to have a sensory function and are used in the detection and capture of prey48. In V. infernalis, they serve to transport the food proximally along the arms to the mouth27. The greater diameters of cirri, and placement along the entire arm in V. rhodanica (Fig. 4), suggests an increased sensory function in these fossil forms.The shape of the arms also contributes to the suction potential49 in coleoids. Functional analysis in Octopoda highlights a positive correlation between distal tapering of the arms and their flexibility. A tapered, flexible arm facilitates more precise adhesion than a cylindrical-shaped one and requires a greater force for sucker detachment49. Suckers detach sequentially, rather than the more simultaneous release observed in models of arms with less taper variation. The tapered diameter of the suckers, like those seen on the sessile arms of V. rhodanica, potentially facilitated this kind of sequential detachment49 allowing them more adherence force and flexibility. Though V. rhodanica has just two suckers on the distal tips of their dorsal arms, the most distal is marginally smaller in diameter than the proximal one. On the dorsal arms, this tapering is observed in conjunction with a well-developed axial nerve cord (Fig. 2b). In extant forms, the nerve cord facilitates complex motor functions42. The combination of these characters in V. rhodanica suggests their arms had increased potential to be actively used in prey capture50 over extant V. infernalis.Though arm crown characters offer insight on the ecology of V. rhodanica, in fossil coleoid phylogenies only a few characters are based on the suckers1, 3. Two studies that have attempted to create a phylogeny using morphological characters that include both fossil and extant taxa return V. rhodanica and V. infernalis as sister taxa1,3. These matrices are, by necessity, heavily influenced by the gladius51 and more than 50% of the characters are based on this feature1,3. Indeed, the authors1 note that the lack of gladius data for some fossil forms, including V. rhodanica, creates an inherent bias in the phylogenetic matrix. Fischer & Riou16 suggested that V. rhodanica and V. infernalis are related on the basis of the observable morphological characters in the family Vampyroteuthidae, though without morphological information on the gladius, a recent systematic synthesis of fossil Octobrachia22 positioned V. rhodanica as Vampyromorpha Incertae sedis.X-ray CT analysis in this study did not allow a reconstruction of the gladius. Nevertheless, it does provide new data on soft tissues, and permits comparisons between extant and fossil taxa. Specifically, we can add distinct states to 4 of the 132 characters in the existing phylogenetic matrix from Sutton et al.3 that was modified and used in Kruta et al.1. These four characters (#89–#92) represent the suckers, and sucker attachments. Detailed examination revealed that the sessile and dorsal arms have the Vampyroteuthis-like attachment. In the dorsal arms, this is more elongated, though it cannot be considered pedunculate like those seen in modern decabrachians. Indeed, the attachment type (plug and base34) is the same, only the length varies. As previously discussed, this variation may have functional implications.When updated with these new data, the matrix from this study returns the same topology seen in Kruta et al.1 that supports the positioning of V. rhodanica and V. infernalis as sister taxa. Further, it strengthens their relationship as they both share a sucker attachment that is not clearly attached to the arm muscles, a state that was previously considered autapomorphic in V. infernalis. However, it is important to note that no additional characters were added for the gladius, which is the cornerstone of coleoid systematics52. Indeed, just 29 of the 132 matrix characters can so far be coded for V. rhodanica, with only 9 of these relating to the 74 states of the gladius.Assuming the phylogenetic work so far is correct, then both species belong to the family Vampyromorphina, and are joined by the Oligocene fossil Necroteuthis hungarica29. While the lack of gladius characters precludes a full phylogenetic understanding of this group, preservation and observation of the soft tissues allow us to infer information regarding palaeobiology.The data obtained in this study demonstrates that the characters observed in V. infernalis, including the sucker attachments and lack of ink sac, were present in Jurassic Vampyromorpha. Comparative anatomy of V. rhodanica and extant V. infernalis revealed that the fossil taxon displayed more morphological variation and were more diversified than previously understood. The assemblage of characters observed in V. rhodanica are consistent with a pelagic predatory lifestyle and corroborate the likelihood of a distinctly different ecological niche. These findings support the hypothesis that a shift towards a deep-sea environment occurred prior to the Oligocene5,29. More

Life cycle cost modelAn analysis of life cycle costs is an economic analysis of the assessment of the total cost of acquisition, ownership and liquidation of a product. It is applicable during the entire life cycle of the product or a life cycle stage or combination of different stages21 and22.There are five period phases of the vehicle life cycle:Generally, the total costs for the above listed phases are acquisition costs, ownership costs and liquidation costs21 and22. For the LCC model, I recommend to divide the life cycle costs into four categories:$$LCC={C}_{P}+{C}_{M}+{C}_{O}+{C}_{D},$$
(1)
$${LCC}_{s}=frac{LCC}{t},$$
(2)

where LCC—the life cycle cost of vehicles, LCCs—the specific life cycle cost of vehicles, CP—the vehicle purchase cost, CM—the maintenance cost, CO—operating state of vehicle cost, CD—the vehicle disposal cost, t—the time of vehicle operation.The model for evaluating the economic viability of products is based on the general LCC model which is based on acquisition and ownership costs$$LCC={C}_{P}+{C}_{OW},$$
(3)

where CP—purchase cost, COW—ownership costs.Acquisition cost (CP) is represented by the purchase price at the time of acquisition of the assessed passenger vehicle.Ownership cost (COW) is significant during the life cycle of a motor vehicle and varies according to the type of the vehicle. This cost includes the costs of maintenance and operation time can be defined as follows10$${C}_{Ow}={C}_{M}+{C}_{O},$$
(4)

where CM—cost of maintenance, CO—operation cost.The cost of ownership a vehicle (COW) can be defined as follows$${C}_{OW}={C}_{O}+{C}_{MC}+{C}_{MP},$$
(5)

where CO—operation cost, CMC—corrective maintenance cost, CMP—preventive maintenance cost.The cost of ownership (COW) may include the operating and maintenance costs which consist of the corrective maintenance cost (CMC) and the cost of preventive maintenance (CMP) of a motor vehicle.Calculation of operating costsOperating cost CO is determined by the price and amount consumed of conventional or alternative types of fuel. It cover the cost of fuel CF, operating fluids, oils and lubricants COL that are supplied during vehicle operation (not during service inspection), tyres CT, accumulator batteries CAB, vehicle insurance fee and road tax or other mandatory fees CIRT, cost of the motorway tax sticker CMT, mandatory vehicle inspection and emission measurement in special vehicles CETC. The costs are calculated according to$${C}_{O}={C}_{F}+{C}_{OL}+{C}_{T}+{C}_{AB}+{C}_{IRT}+{C}_{MT}+{C}_{ETC}.$$
(6)
Fuel costs (CF) are affected by the average consumption of a given type of propulsion vehicle. Then the comparative fuel costs (CF) can be expressed by the equation$${C}_{F}=frac{{bar{c}}_{aF}}{100}{p}_{F}{t}_{l},$$
(7)

where CF—total fuel costs (EUR), (bar{c})aF—average fuel consumption (l/100 km), pF—fuel price (EUR/l), tl—service life of a passenger vehicle (km).Costs for operating fluids, oils and lubricants (COL) are any costs for operating fluids, oils and lubricants that are replenished during operation and not during service maintenance; it can be expressed by the equation$${C}_{OL}=frac{{bar{c}}_{aOL}}{100}{p}_{OL}{t}_{l},$$
(8)

where (bar{c})aOL—average consumption of oil and lubricant (l/100 km), pOL—price of oil and lubricant (EUR/l).The cost of tyres (CT) can be expressed by the equation$${C}_{T}=frac{{t}_{l}}{{bar{d}}_{aT}}{n}_{T}{p}_{T},$$
(9)

where (bar{d})aT—average life of a passenger vehicle tyre (km), nt—number of tyres on the passenger vehicle (pc), pT—price of one piece of tyre (EUR).Accumulator battery costs (CAB) —can be expressed by the equation$${C}_{AB}=frac{{t}_{l}}{{bar{d}}_{aAB}}{n}_{AB}{p}_{AB},$$
(10)

where (bar{d}_{aB})—average life of one accumulator battery (km), nAB—number of accumulator batteries in the passenger vehicle (pc), pAB—price of an accumulator battery (EUR).Costs arising from laws (CIRT) are the costs of motor vehicle insurance (compulsory liability, accident insurance, or other). Some of them can be omitted in case of the same costs due to the simplification of the model. Otherwise, they can be expressed by the equation$${C}_{IRT}=left({C}_{SI}+{C}_{AI}+{C}_{RT}+{C}_{R}right){t}_{la},$$
(11)
where CS1—price of mandatory annual insurance of a passenger vehicle (EUR), CA1—price of the annual accident insurance of a passenger vehicle (EUR), CRT—price of annual road tax (EUR), CR—price of statutory fee (EUR), tla—operating time of the passenger vehicle until decommissioning (years).The cost of obtaining a motorway sticker (CMT) may be omitted if the same type of passenger vehicle is compared. Otherwise, the cost of a motorway sticker (CMT) can be expressed by the equation$${C}_{MT}={c}_{MT}{t}_{la},$$
(12)

where cMT—price of annual motorway sticker for the passenger vehicle (EUR).The costs of the mandatory vehicle inspection and emission measurement (CETC) include the costs incurred for the measurement of emissions of the drive engine unit (CE) and for the technical inspection of the passenger vehicle (CTC). For the proposed model, the costs of the mandatory technical inspections and emission measurements can be expressed by the equation$${C}_{ETC}=left({C}_{E}+{C}_{TC}right)frac{{y}_{n}}{{t}_{la}},$$
(13)

where CE—costs related to the measurement of passenger vehicle emissions (EUR), CTC—costs of mandatory technical inspection (EUR), yn—number of years of legal validity of emission measurement and technical condition for the given type of the passenger vehicle (years).Calculation of maintenance costThe total costs for vehicle maintenance CM consist of the cost of preventive maintenance CMP and the cost of corrective maintenance CMC10,11$${C}_{M}={C}_{MC}+{C}_{MP}.$$
(14)
Vehicle maintenance costs include the cost of material and the cost of labour$${C}_{M}={(C}_{MCM}+{C}_{MCL}+{C}_{MCF})+left({C}_{MPM}+{C}_{MPL}+{C}_{MPF}right),$$
(15)

where CM—cumulative maintenance costs, CMC—corrective maintenance costs, CMP—preventive maintenance costs, CMCM—costs of material used for corrective maintenance, CMCL—costs of labour force for corrective maintenance, CMCF—costs of workshop equipment used for corrective maintenance, CMPM—costs of material used for preventive maintenance, CMPL—costs of labour force for preventive maintenance, CMPF—costs of workshop equipment used for preventive maintenance.

Preventive maintenance costs (CMP) are costs that include all costs associated with preventive maintenance performed to reduce degradation and mitigate the likelihood of failure. At present, preventive maintenance is performed at predetermined time intervals (according to the manufacturer’s preventive maintenance program) or when a specified number of kilometres are not covered before the next service maintenance, depending on the time. In practice, for passenger cars, it is usually 1 or 2 years, depending on the use of engine oil. This mainly includes the cost of:

material consumed during preventive maintenance,

work spent on preventive maintenance,

workshop equipment, training of preventive maintenance specialists.$${C}_{MP}=frac{{t}_{l}}{MTB{M}_{p}}left({C}_{MPM}+{(bar{c}}_{p}{bar{t}}_{pm})right),$$
(17)

where MTBMp—mean operating time between preventive maintenances (km), CMPM—costs of material used for preventive maintenance (EUR), (bar{c})p—average hourly cost of labour and workshop equipment used for maintenance (EUR/hour), ̅tpm—mean time of labour-intensity per one preventive maintenance (hour).

Design of a model for the analysis of selected life cycle costs of a passenger motor vehicleThe model for performing an analysis of life cycle costs for the purchase of a new motor vehicle is based on the basic Eq. (3), (18). We will not count the costs of improvement (CE) and the costs of the decommissioning phase (CD) for the mentioned model due to the calculations of costs that are unnecessary for the analysis. Then the model can be expressed as follows$$LCC={C}_{P}+{C}_{O}+{C}_{M}.$$
(18)
Then, the following Eqs. (6), (7), (8), (9), (10), (11), (12), (13), (16) and (17) are substituted into the given equation, and the selected costs can be calculated for individual vehicles. The resulting model for calculating the LCC costs has the following form$$LCC={C}_{p}+frac{{bar{c}}_{aF}}{100}{p}_{F}{t}_{l}+frac{{bar{c}}_{aOL}}{100}{p}_{OL}{t}_{l}+frac{{t}_{l}}{{bar{d}}_{aT}}{n}_{T}{p}_{T}+frac{{t}_{l}}{{bar{d}}_{aAB}}{n}_{AB}{p}_{AB}+{C}_{SI}{t}_{la}+{c}_{MT}{t}_{la}+left({C}_{E}+{C}_{TC}right)frac{{y}_{n}}{{t}_{la}}+frac{{t}_{l}}{MTBF}left({bar{c}}_{m}+{(bar{c}}_{p}{bar{t}}_{pc})right)+frac{{t}_{l}}{MTB{M}_{p}}left({C}_{OMPM}+{bar{(c}}_{p}{bar{t}}_{pm})right).$$
(19)
It is presented in a Fig. 6.Figure 6Structure of model input parameters for LCC model calculation.Full size imageIn this way, the cumulative costs for each passenger motor vehicle are calculated. Since the passenger motor vehicles may have a different service life tl which is expressed in kilometres, it is recommended to convert this equation to specific costs which are related to one kilometre of use. The selected LCCS life cycle specific costs can be expressed by the following equation$${LCC}_{S}=frac{LCC}{{t}_{l}}.$$
(20)
LCC model input values and items affecting ownership costs for alternative drivesThe process of the calculation of selected life cycle costs for the propulsion of passenger vehicles and the structure of individual cost items is shown in Fig. 6. These are the input parameters to the LCC model.The total life cycle costs are divided into two main cost groups, which are the ownership and acquisition costs for a given drive type. Fuel costs are determined by the price and the quantity of conventional or alternative fuel consumed. For the calculation of the selected LCCs, the authors of the paper assume that the availability of conventional and alternative fuels is not limited in any way. It is assumed that the availability of fuels is ideal, which is not entirely true in practice. This is dependent on the support for each alternative fuel in each state.In practice, therefore, multiple costs may arise due to the distance to the refuelling station to provide alternative fuels such as E85, CNG, LPG and hydrogen. In addition, there is a distance to the charging station for electric drives.Another item that affects the cost of operation for hybrid passenger vehicles is the percentage of alternative fuel driving, which can have a significant impact on life cycle costs. Values for this item are given as a percentage, which is then converted into the number of kilometres driven on alternative and conventional fuel.One of the important parameters for calculating the life cycle operating costs for the hybrid-electric and electric drive is the setting of a threshold value for the capacity of the electric vehicle battery (EV battery) when the replacement is performed. For the model calculation, a limit value of 70% of the electric vehicle battery capacity at 20 °C was set.Return on investmentReturn on investment (ROI) is a performance measure used to evaluate the efficiency or profitability of an investment or compare the efficiency of a number of different investments. ROI tries to directly measure the amount of return on a particular investment, relative to the investment’s cost. To calculate ROI, the benefit (or return) of an investment is divided by the cost of the investment. The result is expressed as a percentage or a ratio12,23.For our calculation of the return on investment ROI on alternative and conventional passenger car propulsion the following formula is used, which is expressed as a percentage$$ROI=frac{{LCC}_{A}-{LCC}_{C}}{{LCC}_{C}}100,$$
(21)

where LCCA—selected live cycle costs of the alternative passenger car propulsion (EUR), LCCC—selected live cycle costs of the conventional passenger car propulsion (EUR).The return on investment of an alternative vehicle ROIAV purchase expresses after how many kilometres the increased cost of purchasing an alternative fuel vehicle compared to a conventional one is recovered. If the value is negative, the payback will not occur for various reasons. The following equation is used to calculate ROIAV$${ROI}_{AV}=frac{{C}_{{P}_{AV}}-{C}_{{P}_{CV}}}{frac{{C}_{O{W}_{CV}}-{C}_{O{W}_{AV}}}{{t}_{l}}}$$
(22)

where ({C}_{{P}_{AV}})—purchase cost on alternative vehicle (EUR), ({C}_{{P}_{CV}})—purchase cost on conventional vehicle (EUR), ({C}_{O{W}_{CV}})—ownership cost on conventional vehicle (EUR), ({C}_{O{W}_{AV}})—ownership cost on alternative vehicle (EUR), tl—service life of the passenger vehicle (km).Ownership costs on conventional vehicle are expressed by the following equation$${C}_{{OW}_{CV}}={left(frac{{bar{c}}_{aF}}{100}{p}_{F}{t}_{l}+frac{{bar{c}}_{aOL}}{100}{p}_{OL}{t}_{l}+frac{{t}_{l}}{{bar{d}}_{aT}}{n}_{T}{p}_{T}+frac{{t}_{l}}{{bar{d}}_{aAB}}{n}_{AB}{p}_{AB}+{C}_{SI}{t}_{la}+{c}_{MT}{t}_{la}+left({C}_{E}+{C}_{TC}right)frac{{y}_{n}}{{t}_{la}}+frac{{t}_{l}}{MTBF}left({bar{c}}_{m}+{(bar{c}}_{p}{bar{t}}_{pc})right)+frac{{t}_{l}}{MTB{M}_{p}}left({C}_{OMPM}+({bar{c}}_{p}{bar{t}}_{pm})right)right)}_{CV}.$$
(23)
Ownership costs on alternative vehicle are expressed by the following equation$${C}_{{OW}_{AV}}={left(frac{{bar{c}}_{aF}}{100}{p}_{F}{t}_{l}+frac{{bar{c}}_{aOL}}{100}{p}_{OL}{t}_{l}+frac{{t}_{l}}{{bar{d}}_{aT}}{n}_{T}{p}_{T}+frac{{t}_{l}}{{bar{d}}_{aAB}}{n}_{AB}{p}_{AB}+{C}_{SI}{t}_{la}+{c}_{MT}{t}_{la}+left({C}_{E}+{C}_{TC}right)frac{{y}_{n}}{{t}_{la}}+frac{{t}_{l}}{MTBF}left({bar{c}}_{m}+{(bar{c}}_{p}{bar{t}}_{pc})right)+frac{{t}_{l}}{MTB{M}_{p}}left({C}_{OMPM}+({bar{c}}_{p}{bar{t}}_{pm})right)right)}_{AV}.$$
(24)
The rate of return on investment for the purchase of an alternative vehicle depending on the kilometres travelled to is expressed by the following equation$${ROI}_{AV({t}_{o})}={(C}_{{P}_{AV}}-{C}_{{P}_{CV}})-({C}_{O{W}_{CV}left({t}_{o}right)}-{C}_{O{W}_{AV}left({t}_{o}right)}) quad text{when} ;to = (0-tl)$$
(25)

where to—operation of the passenger vehicle (km). More

Experimental designThe comparative field microcosm experiments were conducted on Laojun Mountain in China (26.6959 N; 99.7759 E) in September–October 2013, and on Balggesvarri Mountain in Norway (69.3809 N; 20.3483 E) in July 2013, designed to be broadly representative of subtropical and subarctic climatic zones, respectively, as first reported in Wang et al.29. In the Laojun Mountain region, mean annual temperatures ranged from 4.2 to 12.9 °C, with July mean temperatures of 17–25 °C. In the Balggesvarri Mountain region, mean annual temperatures ranged from −2.9 to 0.7 °C, with July mean temperatures of 8–16 °C. The experiments were characterised by an aquatic ecosystem with consistent initial DOM composition but different locally colonised microbial communities and newly produced endogenous DOM. While allowing us to minimise the complexity of natural ecosystems, the experiment provided a means for investigating DOM-microbe associations at large spatial scales by controlling the initial DOM supply. Briefly, we selected locations with five different elevations on each mountainside. The elevations were 3822, 3505, 2915, 2580 and 2286 m a.s.l. on Laojun Mountain in China, and 750, 550, 350, 170 and 20 m a.s.l. on Balggesvarri Mountain in Norway. At each elevation, we established 30 aquatic microcosms (1.5 L bottle) composed of 15 g of sterilised lake sediment and 1.2 L of sterilised artificial lake water at one of ten nutrient levels of 0, 0.45, 1.80, 4.05, 7.65, 11.25, 15.75, 21.60, 28.80 and 36.00 mg N L−1 of KNO3 in the overlying water. To compensate for nitrate additions shifting stoichiometric ratios, KH2PO4 was added to the bottles so that the N/P ratio of the initial overlying water was 14.93, which was similar to the annual average ratio in Taihu Lake during 2007 (that is, 14.49). Thus, we use “nutrient enrichment” to indicate a series of targeted nutrient levels of both nitrate and phosphate, the former of which was used to represent nutrient enrichment in the statistical analyses. Each nutrient level was replicated three times. The lake sediments were obtained from the centre of Taihu Lake, China, and were aseptically canned per bottle after autoclaving at 121 °C for 30 min. Nutrient levels for the experiments were selected based on conditions of the eutrophic Taihu Lake, and the highest nitrate concentration was based on the maximum total nitrogen in 2007 (20.79 mg L−1; Fig. S19). We chose the nutrient level of this year because a massive cyanobacteria bloom in Taihu Lake happened in May 2007 and initiated an odorous drinking water crisis in the nearby city of Wuxi.The microcosms were left in the field for one month allowing airborne bacteria to freely colonise the sediments and water. To keep the microbial dispersal events as natural as possible, we did not cover the experimental microcosms in case of rainfall. To avoid or minimize potential influence of extreme nature events, we (i) left the top 20% of each microcosm empty to prevent water from overflowing during heavy rains, and (ii) checked the experimental sites twice during each experimental period, and added sterilized water to obtain a final volume of approximately 1.2 L. The bottom of our microcosm was buried into the local soils by 10% of the bottle height, partly to reduce UV exposure to sediments. More considerations of the experimental design were detailed in the Supplementary Methods. To avoid the effects of daily temperature variation, we measured the water temperature and pH within 2 h before noon at all elevations in the day before the final sample collection. At the end of the experimental period, we aseptically sampled the water and sediments of the 300 bottles (that is, 2 mountains × 5 elevations × 10 nutrient levels × 3 replicates) for the following analyses of physiochemical variables, bacterial community and DOM composition.Physiochemical variables and bacterial communityWe measured environmental variables, namely, the total nitrogen (TN), total phosphorus (TP), dissolved nutrients (that is, NOx−, NO2−, NH4+ and PO43−), total organic carbon (TOC), dissolved organic carbon (DOC) and chlorophyll a (Chl a) in the sediments, and the NO3−, NO2−, NH4+, PO43− and pH in the overlying water (Table S2, Fig. S20), according to Wang et al.29.The sediment bacteria were examined using high-throughput sequencing of 16S rRNA genes. The sequences were processed in QIIME (v1.9)45 and OTUs were defined at 97% sequence similarity. The bacterial sequences were rarefied to 20,000 per sample. Further details on physicochemical and bacterial community analyses are available in Wang et al.29.ESI FT-ICR MS analysis of DOM samplesHighly accurate mass measurements of DOM within the sediment samples were conducted using a 15 Tesla solariX XR system, a ultrahigh-resolution Fourier transform ion cyclotron resonance mass spectrometer (FT-ICR MS, Bruker Daltonics, Billerica, MA) coupled with an electrospray ionization (ESI) interface, as demonstrated previously46 with some modifications. It should be noted that FT-ICR MS does not identify molecules, but only molecular formulae in terms of elemental composition and there can be many molecular structures sharing the same elemental compositions. DOM was solid-phase extracted (SPE) with Agilent VacElut resins before FT-ICR MS measurement47 with minor modifications. Briefly, an aliquot of 0.7 g freeze-dried sediment was sonicated with 30 ml ultrapure water for 2 h, and centrifuged at 5000 × g for 20 min. The extracted water was filtered through the 0.45 μm Millipore filter and further acidified to pH 2 using 1 M HCl. Cartridges were drained, rinsed with ultrapure water and methanol (ULC-MS grade), and conditioned with pH 2 ultrapure water. Calculated volumes of extracts were slowly passed through cartridges based on DOC concentration. Cartridges were rinsed with pH 2 ultrapure water and dried with N2 gas. Samples were finally eluted with methanol into precombusted amber glass vials, dried with N2 gas and stored at −20 °C until DOM analysis. The extracts were continuously injected into the standard ESI source with a flow rate of 2 μl min−1 and an ESI capillary voltage of 3.5 kV in negative ion mode. One hundred single scans with a transient size of 4 mega word (MW) data points, an ion accumulation time of 0.3 s, and within the mass range of m/z 150–1200, were co-added to a spectrum with absorption mode for phase correction, thereby resulting in a resolving power of 750,000 (FWHM at m/z 400). All FT-ICR mass spectra were internally calibrated using organic matter homologous series separated by 14 Da (-CH2 groups). The mass measurement accuracy was typically within 1 ppm for singly charged ions across a broad m/z range (150–1200 m/z).Data Analysis software (BrukerDaltonik v4.2) was used to convert raw spectra to a list of m/z values using FT-MS peak picker with a signal-to-noise ratio (S/N) threshold set to 7 and absolute intensity threshold to the default value of 100. Putative chemical formulae were assigned using the software Formularity (v1.0)48 following the Compound Identification Algorithm49. In total, 19,538 molecular formulas were putatively assigned for all samples (n = 300) based on the following criteria: S/N  > 7, and mass measurement error  0.80, P ≤ 0.001; Fig. S9). Similar conclusions were also obtained with either OTUs or genera when relating the pairwise distances of molecular traits with SparCC correlation coefficient ρ values among DOM molecules in Fig. 4c. To reduce type I errors in the correlation calculations created by low-occurrence genera or molecules, the majority rule was applied; that is, we retained genera or molecules that were observed in more than half of the total samples (≥75 samples) in China or Norway. The filtered table, including 1340 and 1246 DOM molecules, and 75 and 49 bacterial genera in China and Norway, respectively, was then used for pairwise correlation calculation of DOM and bacteria using SparCC with default parameters35.Finally, bipartite network analysis at a molecular level was performed to quantify the specialization of DOM-bacteria networks (Box 1). The specialization considers interaction abundance and is standardised to account for heterogeneity in the interaction strength and species richness, which describes the levels of “vulnerability” of DOM molecules and “generality” of bacterial taxa27. The threshold correlation for inclusion in bipartite networks was |ρ| = 0.30 to exclude weak interactions and we retained the adjacent matrix with only the interactions between DOM and bacteria. We then constructed two types of interaction networks (i.e., negative and positive networks) based on negative and positive correlation coefficients (SparCC ρ ≤ −0.30 and ρ ≥ 0.30, respectively). According to resource-consumer relationships, negative networks likely indicate the degradation of larger molecules into smaller structures, while positive networks may suggest the production of new molecules via degradation or biosynthetic processes. The SparCC ρ values were multiplied by 10,000 and rounded to integers, and the absolute values were taken for negative networks to enable the calculations of specialization indices. A separate negative and positive sub-network was obtained for each microcosm by selecting the DOM molecules and bacterial taxa in each sample based on its bacterial and DOM compositions. For the network level analysis, we calculated H2′, a measure of specialization27, for each network:$${H}_{2}=-mathop{sum }limits_{i{{mbox{=}}}1}^{i}mathop{sum }limits_{j{{mbox{=}}}1}^{j}{{mbox{(}}}{{{mbox{p}}}}_{{ij}}{{{{{{rm{ln}}}}}}}{{{mbox{p}}}}_{{ij}}{{mbox{)}}}$$
(2)
$${H}_{2}{prime} =frac{{H}_{2{max }}{-}{H}_{2}}{{H}_{2{max }}{-}{H}_{2{min }}}$$
(3)
where ({{{mbox{p}}}}_{{ij}}{{mbox{=}}}{{{mbox{a}}}}_{{ij}}{{mbox{/}}}m), represents the proportion of interactions in a i × j matrix. ({{{mbox{a}}}}_{{ij}}) is the number of interactions between DOM molecule i and bacterial genus j, which is also referred as “link weight”. m is the total number of interactions between all DOM molecules and bacterial genera. H2′ is the standardised H2 against the minimum (H2min) and maximum (H2max) possible for the same distribution of interaction totals. For the molecular level analysis, we calculated the specialization index Kullback–Leibler distance (d′) for DOM molecules (di′) and bacterial genera (dj′), which describes the levels of “vulnerability” of DOM molecules and “generality” of bacterial genera, respectively:$${d}_{i}=mathop{sum }limits_{j=1}^{j}left(frac{{{{mbox{a}}}}_{{ij}}}{{{{mbox{A}}}}_{i}}{{{mbox{ln}}}}frac{{{{mbox{a}}}}_{{ij}}m}{{{{mbox{A}}}}_{i}{{{mbox{A}}}}_{j}}right)$$
(4)
$${d}_{i}{prime} =frac{{d}_{i}-{d}_{{min }}}{{d}_{{max }}-{d}_{{min }}}$$
(5)
where ({A}_{i}) = (mathop{sum }limits_{j{{mbox{=}}}1}^{j}{{{mbox{a}}}}_{{ij}}) and ({A}_{j}) = (mathop{sum }limits_{i{{mbox{=}}}1}^{i}{{{mbox{a}}}}_{{ij}}), are the total number of interactions of DOM molecule i and bacterial genus j, respectively. di′ is the standardised di against the minimum (dmin) and maximum (dmax) possible for the same distribution of interaction totals. The equations of dj′ are analogous to di′, replacing j by i. Weighted means of d′ for DOM were calculated for each network as the sum of the product of d′ for each individual molecule i (di′) and relative intensity Ii divided by the sum of all intensities d′  = Ʃ(di′ × Ii)/Ʃ(Ii). Weighted means of d′ for bacteria were calculated as the sum of the d′ of each individual bacterial genus j (dj′) and relative abundance of bacterial genus Ij divided by the sum of all abundance. All calculations were performed using the R package FD V1.0.12. The observed H2′ and d′ values ranged from 0 (complete generalization) to 1 (complete specialization)28 (Fig. S21). Specifically, elevated H2′ or d′ values indicate a high degree of specialization, while lower values suggest increased generalization, that is, higher vulnerability of DOM and/or higher generality of microbes. To directly compare the network indices across the elevations or nutrient enrichment levels, we used a null modelling approach. We standardised the three observed specialization indices (Sobserved; that is, H2′, d′ of DOM, and d′ of bacteria) by calculating their z-scores63 using the equation:$${z}_{S}=({S}_{{{{{{rm{observed}}}}}}}-overline{{{S}}_{{{{{{rm{null}}}}}}}})/({sigma }_{S_{{{{{rm{null}}}}}}})$$
(6)
where (overline{{{S}}_{{{{{{rm{null}}}}}}}}) and ({sigma }_{S_{{{{{rm{null}}}}}}}) were, respectively, the mean and standard deviation of the null distribution of S (Snull). One hundred randomised null networks were generated for each bipartite network to derive Snull using the swap.web algorithm, which keeps species richness and the number of interactions per species constant along with network connectance. This null model analysis indicates that interactions between DOM and bacteria were non-random as the observed network specialization indices were generally significantly lower than expected by chance (P  0.05), which tests whether the model structure differs from the observed data, high comparative fit index (CFI  > 0.95) and low standardised root mean squared residual (SRMR  More

Wellnhofer, P. The Illustrated Encyclopedia of Pterosaurs (Crescent Books, 1991).Unwin, D. M. The pterosaurs from deep time (Pi Press, 2005).Witton, M. P. Pterosaurs: Natural History (Anatomy (Princeton University Press, 2013).Book 

Google Scholar 
Williams, C. J. et al. Helically arranged cross struts in azhdarchid pterosaur cervical vertebrae and their biomechanical implications. iScience 24, 102338 (2021).PubMed 
PubMed Central 
Article 

Google Scholar 
Bestwick, J., Unwin, D. M., Butler, R. J. & Purnell, M. A. Dietary diversity and evolution of the earliest flying vertebrates revealed by dental microwear texture analysis. Nat. Commun. 11, 5293 (2020).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Ryang, W. H. Characteristics of strike-slip basin formation and sedimentary fills and the Cretaceous small basins of the Korean Peninsula. J. Geo. Soc. Korea 49, 31–45 (2013).CAS 

Google Scholar 
Kim, B. G. & Park, B. G. Geological report of the Dongbok sheet (1:50,000) (Geological Survey of Korea, Seoul, 1966).Lee, H., Sim, M. S. & Choi, T. Stratigraphic evolution of the northern part of the Cretaceous Neungju basin South Korea. Geosci. J. 23, 849–865 (2019).CAS 
Article 

Google Scholar 
Paik, I. S., Huh, M., So, Y. H., Lee, J. E. & Kim, H. J. Traces of evaporites in Upper Cretaceous lacustrine deposits of Korea: Origin and paleoenvironmental implications. J. Asian Earth Sci. 30, 93–107 (2007).Article 

Google Scholar 
Cohen, K. M., Finney, S. M., Gibbard, P. L. & Fan, J.-X. The ICS international Chronostratigraphic chart. Episodes 36, 199–204 (2013).Article 

Google Scholar 
Calvo, J. O. & Lockley, M. G. The first pterosaur tracks from Gondwana. Cretac. Res. 22, 585–590 (2001).Article 

Google Scholar 
Kukihara, R. & Lockley, M. G. Fossil footprints from the dakota group (Cretaceous) john martin reservoir, bent county, Colorado: New insights into the paleoecology of the Dinosaur freeway. Cretac. Res. 33, 165–182 (2012).Article 

Google Scholar 
Lockley, M. & Schumacher, B. A new pterosaur swim tracks locality from the Cretaceous Dakota Group of eastern Colorado: implications for pterosaur swim track behavior. Fossil Footprints of Western North America. Bull. NM Mus. Nat. Hist. Sci, 365–371 (2014).Smith, R. E., Martill, D. M., Unwin, D. M. & Steel, L. Edentulous pterosaurs from the Cambridge Greensand (Cretaceous) of eastern England with a review of Ornithostoma Seeley, 1871. Proc. Geol. Assoc. (2020).Ibrahim, N., Unwin, D. M., Martill, D. M., Baidder, L. & Zouhri, S. A new pterosaur (Pterodactyloidea: Azhdarchidae) from the Upper Cretaceous of Morocco. PLoS ONE 5, e10875 (2010).PubMed 
PubMed Central 
Article 

Google Scholar 
Martill, D. M. & Ibrahim, N. An unusual modification of the jaws in cf. Alanqa, a mid-Cretaceous azhdarchid pterosaur from the Kem Kem beds of Morocco. Cretac. Res. 53, 59–67 (2015).Article 

Google Scholar 
Jacobs, M. L., Martill, D. M., Ibrahim, N. & Longrich, N. A new species of Coloborhynchus (Pterosauria, Ornithocheiridae) from the mid-Cretaceous of North Africa. Cretac. Res. 95, 77–88 (2019).Article 

Google Scholar 
Jacobs, M. L. et al. New toothed pterosaurs (Pterosauria: Ornithocheiridae) from the middle Cretaceous Kem Kem beds of Morocco and implications for pterosaur palaeobiogeography and diversity. Cretac. Res. 110, 104413 (2020).Article 

Google Scholar 
McPhee, J. et al. A new ? Chaoyangopterid (Pterosauria: Pterodactyloidea) from the Cretaceous Kem Kem beds of southern Morocco. Cretac. Res. 110, 104410 (2020).Article 

Google Scholar 
Martill, D. M. et al. A new tapejarid (Pterosauria, Azhdarchoidea) from the mid-Cretaceous Kem Kem beds of Takmout, southern Morocco. Cretac. Res. 112, 104424 (2020).Article 

Google Scholar 
Martill, D. M., Unwin, D. M., Ibrahim, N. & Longrich, N. A new edentulous pterosaur from the Cretaceous Kem Kem beds of south eastern Morocco. Cretac. Res. 84, 1–12 (2018).Article 

Google Scholar 
Smith, R. E. et al. Small, immature pterosaurs from the Cretaceous of Africa: implications for taphonomic bias and palaeocommunity structure in flying reptiles. Cretac. Res. 130, 105061 (2022).Article 

Google Scholar 
Smith, R. E., Martill, D. M., Kao, A., Zouhri, S. & Longrich, N. A long-billed, possible probe-feeding pterosaur (Pterodactyloidea: ?Azhdarchoidea) from the mid-Cretaceous of Morocco North Africa. Cretac. Res. 118, 104643 (2021).Article 

Google Scholar 
Kellner, A. W. A. et al. First complete pterosaur from the Afro-Arabian continent: insight into pterodactyloid diversity. Sci. Rep. 9, 17875 (2019).PubMed 
PubMed Central 
Article 

Google Scholar 
Elgin, R. A. & Frey, E. A new azhdarchoid pterosaur from the Cenomanian (Late Cretaceous) of Lebanon. Swiss J. Geosci. 104, 21–33 (2011).Article 

Google Scholar 
Averianov, A. O., Kurochkin, E. N., Pervushov, E. M. & Ivanov, A. V. Two bone fragments of ornithocheiroid pterosaurs from the Cenomanian of Volgograd Region, southern Russia. Acta Palaeontol. Pol. 50 (2005).Averianov, A. & Kurochkin, E. A new pterosaurian record from the Cenomanian of the Volga region. Paleontol. J. 44, 695–697 (2010).Article 

Google Scholar 
Nessov, L. Flying reptiles from the Jurassic and cretaceous of the USSR and significance of their remains for the reconstruction of paleogeographical conditions. Vestn. Leningr. Gos. Univ. Ser. 7, 28 (1990).
Google Scholar 
Bakhurina, N. N. & Unwin, D. M. A survey of pterosaurs from the Jurassic and Cretaceous of the former Soviet Union and Mongolia. (1995).Averianov, A. O. New records of azhdarchids (Pterosauria, Azhdarchidae) from the Late Cretaceous of Russia, Kazakhstan, and Central Asia. Paleontol. J. 41, 189–197 (2007).Article 

Google Scholar 
Averianov, A. Mid-Cretaceous ornithocheirids (Pterosauria, Ornithocheiridae) from Russia and Uzbekistan. Paleontol. J. 41, 79–86 (2007).Article 

Google Scholar 
Huh, M., Paik, I. S., Chung, C. H., Hwang, K. G. & Kim, B. S. Theropod tracks from Seoyuri in Hwasun, Jeollanamdo, Korea: occurrence and paleontological significance. J. Geo. Soc. Korea 39, 461–478 (2003).CAS 

Google Scholar 
Huh, M. et al. Well-preserved theropod tracks from the Upper Cretaceous of Hwasun County, southwestern South Korea, and their paleobiological implications. Cretac. Res. 27, 123–138 (2006).Article 

Google Scholar 
Lockley, M. G., Huh, M. & Kim, B. S. Ornithopodichnus and pes-only sauropod Trackways from the Hwasun tracksite Cretaceous of Korea. Ichnos 19, 93–100 (2012).Article 

Google Scholar 
Hwang, K. G., Huh, M. & Paik, I. S. A unique trackway of small theropod from Seoyu-ri, Hwasun-gun Jeollanam province. J. Geo. Soc. Korea 42, 69–78 (2006).CAS 

Google Scholar 
Kim, B. S. & Huh, M. Analysis of the acceleration phase of a theropod dinosaur based on a Cretaceous trackway from Korea. Palaeogeogr. Palaeoclimatol. Palaeoecol. 293, 1–8 (2010).Article 

Google Scholar 
Marchetti, L. et al. Defining the morphological quality of fossil footprints. Problems and principles of preservation in tetrapod ichnology with examples from the Palaeozoic to the present. Earth-Sci. Rev. 193, 109–145 (2019).Article 

Google Scholar 
Rodríguez-de La Rosa, R. A. Pterosaur tracks from the latest Campanian Cerro del Pueblo formation of southeastern Coahuila. Mexico. Geol. Soc. Spec. Publ. 271, 275–282 (2003).Article 

Google Scholar 
Lockley, M. G. & Meyer, C. Crocodylomorph trackways from the Jurassic to early cretaceous of North America and Europe: Implications for Ichnotaxonomy. Ichnos 11, 167–178 (2004).Article 

Google Scholar 
Ambroggi, R. & De Lapparent, A. Les empreintes de pas fossiles du Maestrichtien d’Agadir. Notes du Service Géologique du Maroc 10, 43–57 (1954).
Google Scholar 
Stokes, W. L. Pterodactyl tracks from the Morrison Formation. J. Paleontol. 31, 952–954 (1957).
Google Scholar 
Delair, J. Note on Purbeck fossil footprints, with descriptions of two hitherto unknown forms from Dorset. Proceedings of the Dorset Natural History and Archaeological Society. 92–100 (1963).Hwang, K.-G., Huh, M. I. N., Lockley, M. G., Unwin, D. M. & Wright, J. L. New pterosaur tracks (Pteraichnidae) from the Late Cretaceous Uhangri Formation, southwestern Korea. Geol. Mag. 139, 421–435 (2002).Article 

Google Scholar 
Mazin, J.-M. & Pouech, J. The first non-pterodactyloid pterosaurian trackways and the terrestrial ability of non-pterodactyloid pterosaurs. Geobios 58, 39–53 (2020).Article 

Google Scholar 
Masrour, M., de Ducla, M., Billon-Bruyat, J.-P. & Mazin, J.-M. Rediscovery of the Tagragra tracksite (Maastrichtian, Agadir, Morocco): Agadirichnus elegans Ambroggi and Lapparent 1954 is Pterosaurian Ichnotaxon. Ichnos 25, 285–294 (2018).Article 

Google Scholar 
Wright, J. L., Unwin, D. M., Lockley, M. G. & Rainforth, E. C. Pterosaur tracks from the Purbeck limestone formation of Dorset England. Proc. Geol. Assoc. 108, 39–48 (1997).Article 

Google Scholar 
Lockley, M. G. et al. The fossil trackway Pteraichnusis pterosaurian, not crocodilian: Implications for the global distribution of pterosaur tracks. Ichnos 4, 7–20 (1995).Article 

Google Scholar 
Billon-Bruyat, J.-P. & Mazin, J.-M. The systematic problem of tetrapod ichnotaxa: the case study of Pteraichnus Stokes, 1957 (Pterosauria, Pterodactyloidae). Geol. Soc. Spec. Publ. 217, 315–324 (2003).Article 

Google Scholar 
Pascual Arribas, C. & Sanz Pérez, E. Huellas de Pterosaurios en el grupo Oncala (Soria, España). Pteraichnus palaciei-saenzi, nov. icnosp. Estudios Geol. 56, 73–100 (2000).
Google Scholar 
Calvo, M. M., Vidarte, C. F., Fuentes, F. M. & Fuentes, M. M. Huellas de Pterosaurios en la Sierra de Oncala (Soria, España). Nuevas icnoespecies: pteraichnus vetustior, Pteraichnus parvus. Pteraichnus manueli. Celtiberia 54, 471–490 (2004).
Google Scholar 
Fuentes Vidarte, C., Meijide Calvo, M., Meijide Fuentes, F. & Meijide Fuentes, M. Pteraichnus longipodus nov. icnosp. en la Sierra de Oncala (Soria, España). Studia Geologica Salmanticensia, 103–114 (2004).Peng, B.-X., Du, Y.-S., Li, D.-Q. & Bai, Z.-C. The first discovery of the early Cretaceous Pterosaur track and its significance in Yanguoxia, Yongjing County, Gansu Province. Earth Sci.-J. China Univ. Geosci. 29, 21–24 (2004).
Google Scholar 
Lee, Y.-N., Lee, H.-J., Lü, J. & Kobayashi, Y. New pterosaur tracks from the Hasandong formation (Lower Cretaceous) of Hadong County South Korea. Cretac. Res. 29, 345–353 (2008).Article 

Google Scholar 
Lee, Y.-N., Azuma, Y., Lee, H.-J., Shibata, M. & Lü, J. The first pterosaur trackways from Japan. Cretac. Res. 31, 263–273 (2010).Article 

Google Scholar 
Chen, R. et al. Pterosaur tracks from the early late cretaceous of Dongyang City, Zhejiang Province China. Geol. Bull. China. 32, 693–698 (2013).CAS 

Google Scholar 
Li, Y., Wang, X. & Jiang, S. A new pterosaur tracksite from the Lower Cretaceous of Wuerho, Junggar Basin, China: inferring the first putative pterosaur trackmaker. PeerJ 9, e11361 (2021).PubMed 
PubMed Central 
Article 

Google Scholar 
Ha, S. et al. Diminutive pterosaur tracks and trackways (Pteraichnus gracilis ichnosp. Nov.) from the lower Cretaceous Jinju formation, Gyeongsang basin. Korea. Cretac. Res. 131, 105080 (2021).Article 

Google Scholar 
Sánchez-Hernández, B., Przewieslik, A. G. & Benton, M. J. A reassessment of the Pteraichnus ichnospecies from the early Cretaceous of Soria Province Spain. J. Vertebr. Paleontol. 29, 487–497 (2009).Article 

Google Scholar 
Zhou, X. et al. A new darwinopteran pterosaur reveals arborealism and an opposed thumb. Curr. Biol. 31, 2429-2436.e2427 (2021).CAS 
PubMed 
Article 

Google Scholar 
Lü, J. et al. Dragons of the Skies (recent advances on the study of pterosaurs from China) (Zhejiang Science and Technology Press, 2013).
Google Scholar 
Beccari, V. et al. Osteology of an exceptionally well-preserved tapejarid skeleton from Brazil: Revealing the anatomy of a curious pterodactyloid clade. PLoS ONE 16, e0254789 (2021).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Lü, J. A new boreopterid pterodactyloid pterosaur from the Early Cretaceous Yixian Formation of Liaoning Province, northeastern China. Acta Geologica Sinica-English Edition 84, 241–246 (2010).Article 

Google Scholar 
Bennett, S. C. Terrestrial locomotion of pterosaurs: A reconstruction based on Pteraichnus trackways. J. Vertebr. Paleontol. 17, 104–113 (2010).Article 

Google Scholar 
Wang, X. & Lü, J. Discovery of a pterodactylid pterosaur from the Yixian Formation of western Liaoning China. Chin. Sci. Bull. 46, A3–A8 (2001).Article 

Google Scholar 
Frey, E. et al. A new specimen of nyctosaurid pterosaur, cf. Muzquizopteryx sp. from the Late Cretaceous of northeast Mexico. Revista mexicana de ciencias geológicas 29, 131–139 (2012).
Google Scholar 
Wu, W.-H., Zhou, C.-F. & Andres, B. The toothless pterosaur Jidapterus edentus (Pterodactyloidea: Azhdarchoidea) from the Early Cretaceous Jehol Biota and its paleoecological implications. PLoS ONE 12, e0185486 (2017).PubMed 
PubMed Central 
Article 

Google Scholar 
Lü, J. et al. The toothless pterosaurs from China. Acta Geol. Sin. 90, 2513–2525 (2016).
Google Scholar 
Zhang, X., Jiang, S., Cheng, X. & Wang, X. New Material of Sinopterus (Pterosauria, Tapejaridae) from the Early Cretaceous Jehol Biota of China. An. Acad. Bras. Cienc. 91 (2019).Bestwick, J., Unwin, D. M., Butler, R. J., Henderson, D. M. & Purnell, M. A. Pterosaur dietary hypotheses: A review of ideas and approaches. Biol. Rev. 93, 2021–2048 (2018).PubMed 
Article 

Google Scholar 
Chen, H. et al. New anatomical information on Dsungaripterus weii Young, 1964 with focus on the palatal region. PeerJ 8, e8741 (2020).PubMed 
PubMed Central 
Article 

Google Scholar 
Li, D. et al. A manus dominated pterosaur track assemblage from Gansu, China: Implications for behavior. Sci. Bull. 60, 264–272 (2015).Article 

Google Scholar 
Masrour, M., Pascual-Arribas, C., de Ducla, M., Hernández-Medrano, N. & Pérez-Lorente, F. Anza palaeoichnological site. Late Cretaceous. Morocco. Part I. The first African pterosaur trackway (manus only). J. African Earth Sci. 134, 766–775 (2017).Article 

Google Scholar 
Bramwell, C. D. & Whitfield, G. R. Biomechanics of Pteranodon. Phil. Trans. R. Soc. Lond. B. 267, 503–581 (1974).Article 

Google Scholar 
Bennett, S. C. Terrestrial locomotion of pterosaurs: a reconstruction based on Pteraichnus trackways. J. Vertebr. Paleontol. 17, 104–113 (1997).Article 

Google Scholar 
Mazin, J.-M., Billon-Bruyat, J.-P., Hantzpergue, P. & Lafaurie, G. Ichnological evidence for quadrupedal locomotion in pterodactyloid pterosaurs: Trackways from the Late Jurassic of Crayssac (southwestern France). Geol. Soc. Spec. Publ. 217, 283–296 (2003).Article 

Google Scholar 
Henderson, D. M. Pterosaur body mass estimates from three-dimensional mathematical slicing. J. Vertebr. Paleontol. 30, 768–785 (2010).Article 

Google Scholar 
Lockley, M. G. & Wright, J. L. Pterosaur swim tracks and other ichnological evidnce of behaviour and ecology. Geol. Soc. Spec. Publ. 217, 297–313 (2003).Article 

Google Scholar 
Lockley, M., Mitchell, L. & Odier, G. P. Small Theropod track assemblages from middle Jurassic Eolianites of eastern Utah: Paleoecological insights from dune Ichnofacies in a transgressive sequence. Ichnos 14, 131–142 (2007).Article 

Google Scholar 
Fiorillo, A. R., Hasiotis, S. T., Kobayashi, Y. & Tomsich, C. S. A pterosaur manus track from Denali National park, Alaska Range, Alaska United States. Palaios 24, 466–472 (2009).Article 

Google Scholar 
Bell, P. R., Fanti, F. & Sissons, R. A possible pterosaur manus track from the late Cretaceous of Alberta. Lethaia 46, 274–279 (2013).Article 

Google Scholar 
Stinnesbeck, W. et al. Theropod, avian, pterosaur, and arthropod tracks from the uppermost Cretaceous Las Encinas Formation, Coahuila, northeastern Mexico, and their significance for the end-Cretaceous mass extinction. Geol. Soc. Am. Bull. 129, 331–348 (2017).Article 

Google Scholar 
Xing, L. et al. Late Cretaceous ornithopod-dominated, theropod, and pterosaur track assemblages from the Nanxiong Basin, China: New discoveries, ichnotaxonomy, and paleoecology. Palaeogeogr. Palaeoclimatol. Palaeoecol. 466, 303–313 (2017).Article 

Google Scholar 
Lockley, M. G., Gierlinski, G. D., Adach, L., Schumacher, B. & Cart, K. Newly discovered tetrapod ichnotaxa from the Upper Blackhawk Formation Utah. Bull. N. M. M. Nat. Hist. Sci. 79, 469–480 (2018).
Google Scholar 
Lockley, M. G. & Gillette, D. Pterosaur and bird tracks from a new Late Cretaceous locality in Utah. Verteb. Paleontol. Utah 99, 355–359 (1999).
Google Scholar 
Bennett, S. C. The ontogeny of Pteranodon and other pterosaurs. Paleobiology 19, 92–106 (1993).Article 

Google Scholar 
Bennett, S. C. Year-classes of pterosaurs from the Solnhofen Limestone of Germany: taxonomic and systematic implications. J. Vertebr. Paleontol. 16, 432–444 (1996).Article 

Google Scholar 
Chiappe, L. M., Codorniú, L., Grellet-Tinner, G. & Rivarola, D. Argentinian unhatched pterosaur fossil. Nature 432, 571–572 (2004).CAS 
PubMed 
Article 

Google Scholar 
Codorniú, L., Chiappe, L. & Rivarola, D. Neonate morphology and development in pterosaurs: evidence from a Ctenochasmatid embryo from the Early Cretaceous of Argentina. Geol. Soc. Spec. Publ. 455, 83–94 (2018).Article 

Google Scholar 
Mickelson, D. L., Lockley, M. G., Bishop, J. & Kirkland, J. A New Pterosaur Tracksite from the Jurassic Summerville Formation, near Ferron Utah. Ichnos 11, 125–142 (2004).Article 

Google Scholar  More

Model 1: Evolution of multiple mating and mate monopolisation under ancestral monogamyIn all models, I assume a large population with a 1:1 sex ratio. I begin with what is possibly the simplest model set-up for deriving Bateman functions in a scenario that is completely symmetrical aside from gamete number. Assume a monogamous, externally fertilising population where parents pair up and release their gametes into a nest. That is, every individual in the initial population participates in exactly one fertilisation event (the equivalent of a mating). Now consider a mutant individual that can attract multiple mates of the opposite type to release gametes into its nest, with no competition from other individuals of its own type. This simple set-up avoids asymmetries arising from internal fertilisation, and the complication of direct gamete competition for the multiply mating mutant individual (which is examined in Models 2–3), placing focus directly on the core of the problem: the asymmetry arising in fertilisation from imbalanced gamete numbers. All gametes are released in one burst by all individuals, but the focal individual may achieve ‘multiple matings’ simply by monopolising multiple mates at its nest. The reproductive success of the focal individual is then equivalent to the number of fertilisations that take place in that nest. Our aim is to understand how the reproductive success of an individual deviating from the monogamous population strategy and instead mating with (hat{m}) individuals of the opposite type is altered. A strong positive relationship between (hat{m}) and reproductive success then indicates a steep Bateman gradient. If Bateman’s assertion is correct, the resulting gradient should be steeper for the type that produces the larger number of gametes. Note that there is a game-theoretical25 flavour to this setting, where the focus is on the fitness of a rare mutant in a population with a fixed resident strategy.The two types are labelled with x and y, which could correspond to the two sexes, depending on what gamete numbers are assigned to them. The number of gametes produced by a single individual is labelled nx and ny, and the total number of gametes in a nest (or more generally, a fertilisation arena which could be internal or external) is labelled with Nx and Ny. To compute the number of fertilisations in a nest with a total of Nx and Ny gametes, I use a fertilisation function first derived by Togashi et al.24 purely from biophysical principles, treating the two gamete types symmetrically, with no pre-existing assumptions about differences between females and males or their gametes (for a broader context and comparison to other functions, see Table 1 and function F7 in19). Any sex-specific differences arise only retrospectively after different gamete numbers are assigned to x and y of which either one could be male or female. The fertilisation function is (fleft({N}_{x},{N}_{y}right)={N}_{x}{N}_{y}frac{{e}^{a{N}_{x}}-{e}^{a{N}_{y}}}{{{N}_{x}e}^{a{N}_{x}}-{N}_{y}{e}^{a{N}_{y}}}), where a is a parameter controlling fertilisation efficiency (for the special case Nx = Ny the function is defined as (fleft({N}_{x},{N}_{y}right)=frac{a{N}_{x}^{2}}{1+a{N}_{x}})19,24, which is also the limit of f when Ny → Nx).In a monogamous resident pair, we have simply Nx = nx and Ny = ny. But if a mutant individual of type x is able to attract (hat{m}) fertilisation partners of type y, then for that individual ({N}_{y}=hat{m}{n}_{y}), and the corresponding Bateman function is$${b}_{x}left(hat{m}right)=fleft({N}_{x},{N}_{y}right)=fleft({n}_{x},hat{m}{n}_{y}right)$$
(1)
where the fertilisation function f is as described above. Because of symmetry, the corresponding function for y is found simply by swapping x and y. This function can reproduce the characteristic Bateman gradient asymmetry as gamete numbers diverge (progressing from isogamy to anisogamy in Fig. 1), showing how Bateman’s assertion follows from biophysical effects that arise from unequal numbers of fusing particles: the fertilisation function f is derived solely from such biophysical effects, not from any sex-specific assumptions. Equation (1) makes no reference to sexes, and they only become specified when values are assigned to nx and ny. For example, if nx = 10 and ny = 10,000, the female Bateman function is ({b}_{x}left(hat{m}right)) and the male Bateman function ({b}_{y}left(hat{m}right)), where for the latter all xs in Eq. (1) are replaced with ys and vice versa. The labels x and y are truly just labels. While there are inevitably assumptions built into the equations, crucially we can be certain there are no sex-specific assumptions. Yet the typical shapes reminiscent of Bateman gradients arise from the model when different values are specified for nx and ny (Fig. 1).Fig. 1: The Bateman function of Eq. (1).This figure shows how the basic Bateman gradient asymmetry arises from simple biophysics and mathematics of fertilisation. The population is monogamous aside from a mutant individual, whose number of fertilisation partners (‘matings’) varies on the horizontal axes within panels. a–d show the effect of variation in sex-specific gamete numbers under efficient fertilisation, while e–h show the effect of variation in sex-specific gamete numbers under inefficient fertilisation. Parameter values used are shown in the figure. Females (gamete number nx) are indicated by blue crosses and connecting lines, while males (gamete number ny) are indicated by black dots and connecting lines. Under isogamy, females and males are undefined, and the two colours overlap. The typical sex-specific shapes of Bateman gradients arise from a single equation (which itself is not sex-specific) when a difference in gamete numbers is assigned to nx and ny, confirming Bateman’s intuition that the primary cause of the difference in selection is that females produce fewer gametes than males. The entire range of gamete number ratios presented in the figure is observed in nature, from equal gamete size in many unicellular organisms39 to vertebrates, where sperm count per ejaculate can commonly exceed 109 (see ref. 40 and Supplementary Information therein).Full size imageGamete limitation changes the results quantitatively so that under conditions of poor fertilisation efficiency a larger imbalance in gamete numbers is needed for Bateman gradients to diverge to a similar extent. However, even under inefficient fertilisation, the Bateman gradients do not reverse.Model 2: An external fertiliser model with population-level polygamy and gamete competitionModel 1 presented the simplest possible scenario, where all individuals except a rare mutant mate only once, and gamete competition (sperm competition26, but without assigning either gamete type to be sperm) was thus excluded for the focal mutant individual. Now I generalise from this to a situation that remains entirely symmetrical, but where the resident number of matings can take on any value, and then derive the Bateman function for a rare mutant that deviates from this population-level value. This set-up allows for gamete competition for the focal mutant individual, a crucial addition because of the empirical and theoretical importance of sperm competition26, as well as earlier theory suggesting that polyandry decreases the sex difference in Bateman gradients2.The biological set-up is such that there is a large population and a large number of patches (fertilisation arenas) where multiple individuals of both sexes can release their gametes for fertilisation. After all individuals have released their gametes, those in each patch mix freely and fertilisations take place randomly. Set up in this way, the model is again identical from the perspective of both sexes, and gamete number can be isolated as the sole possible causal factor in any subsequent differences that may arise, extending from the initially monogamous and gamete competition-free scenario of Model 1. All individuals of both sexes are assumed to initially have the same strategy: to divide their nx or ny gametes equally between m patches, and distribute themselves in such a way that gametes from m individuals of each type release gametes into each patch (the number of individuals of each sex per patch need not necessarily be strictly equal to m, but this is the simplest assumption to account for the fact that gamete competition tends to increase with multiple ‘matings’). Now, if a rare x mutant divides its gametes evenly into (hat{m}) randomly selected patches, its gamete number per patch and consequently competitiveness in each patch is altered. Therefore, gametes of a mutant of type x will gain, on average, a fraction ({c}_{x}=left({n}_{x}/hat{m}right)/{N}_{x}) of the fertilisations in that patch, where ({N}_{x}={n}_{x}/hat{m}+(m-1){n}_{x}/m). To compute the number of realised fertilisations in a patch, I use the same fertilisation function as in Model 1, where the mutant number of gametes in a patch is Nx as above and the number of gametes of the opposite type is ({N}_{y}=mfrac{{n}_{y}}{m}={n}_{y}). All the components are now in place to write down the Bateman function corresponding to this scenario, for a mutant of type x:$${b}_{x}left(hat{m},mright)=hat{m}{c}_{x}fleft({N}_{x},{N}_{y}right)$$
(2)
where cx, Nx and Ny are as defined above, and the fertilisation function f is as in Model 1. For completeness, define bx(0, m) = 0, which is necessarily true, but useful to define separately because division by 0 renders Eq. (2) formally undefined when (hat{m}=0).As in Model 1, Eq. (2) makes no reference to sexes, and they only become specified when values are assigned to nx and ny (Fig. 2).Fig. 2: The Bateman function of Eq. (2) for an externally fertilising population with potential for population-wide polygamy and gamete competition.Results are shown for two values of resident matings (m = 1 and m = 2). a–h show the effect of variation in sex-specific gamete numbers and in fertilisation efficiency with m = 1, while i–p show the same with m = 2. Parameter values used are shown in the figure. The value m = 2 is used here because it is comparable to the mean number of matings in Bateman’s1 work (see Fig. 3 for corresponding results with internal fertilisation, but note that the aim of the models is not to quantitatively reproduce Bateman’s results). Females (gamete number nx) are indicated by blue crosses and connecting lines, while males (gamete number ny) are indicated by black dots and connecting lines. Under isogamy, females and males are undefined, and the two colours overlap. Further variation in m is examined in Fig. 4.Full size imageModel 3: An internal fertiliser modelModels 1–2 were set up with the central aim of full symmetry and exclusion of any sex-specific assumptions. Internal fertilisation breaks this symmetry by introducing a sex-specific assumption other than gamete number. Bateman gradients are, however, most commonly applied to situations with internal fertilisation where females are gamete recipients and males are gamete donors27. I therefore construct a model accounting for internal fertilisation. Where Eqs. (1) and (2) allowed no sex differences aside from gamete number, here I additionally consider the fact that females receive gametes while males donate them.As in model 2, there is a very large population, and I assume that in the resident population, all females and males mate exactly m times. It is then considered how a rare mutant individual’s (of either sex) fitness depends on its number of matings (hat{m}).I use the same fertilisation function as in Models 1-2. Consider first the female perspective (labelled with x). A female produces nx gametes and retains them internally. Each female mates with m males, who also mate with m females, dividing their gametes evenly over these matings. Therefore a mutant female receives (hat{m}frac{{n}_{y}}{m}) male gametes, and her reproductive success is$${b}_{x}left(hat{m},mright)=fleft({n}_{x},hat{m}frac{{n}_{y}}{m}right)$$
(3)
A mutant male, on the other hand, mates with (hat{m}) females, each of which mates with m−1 additional males. Therefore, the mutant male’s mating partners will receive a total of ({{N}_{y}=n}_{y}/hat{m}+(m-1){n}_{y}/{m}) male gametes. Thus, the mutant male gains a fraction ({c}_{y}=left({n}_{y}/hat{m}right)/{N}_{y}) of the fertilisations with each female, while the total reproductive success per female is f(nx,Ny). The mutant male’s reproductive success is therefore$${b}_{y}left(hat{m},mright)=hat{m}{c}_{y}fleft({n}_{x},{N}_{y}right)$$
(4)
To avoid division by 0, we can again define by (0, m) = 0, analogous to Model 2. In contrast to Models 1–2, there are now separate equations for each sex because of the additional sex-specific assumption of internal fertilisation, but no further sex-specific assumptions are used in their derivation. Visually the Bateman functions (Fig. 3) are nevertheless very similar to Model 2, and again reproduce the sex-specific shapes first proposed by Bateman1 when fertilisation is efficient. However, an interesting exception arises when relatively weak asymmetry in gamete numbers is combined with inefficient fertilisation and gamete limitation. When these conditions are combined with internal fertilisation, Bateman gradients can theoretically be reversed.Fig. 3: The Bateman functions of Eqs. (3) and (4) for internal fertilisation.Where Figs. 1 and 2 show that the sex-specific shapes of Bateman functions are ultimately caused by differences in gamete number, Fig. 3 shows that internal fertilisation does not invalidate this outcome when fertilisation is efficient. As in Fig. 2, results are shown for two values of resident matings (1 and 2), and the value m = 2 is used because it is comparable to the mean number of matings in Bateman’s1 work. a–h show the effect of variation in sex-specific gamete numbers and in fertilisation efficiency with m = 1, while i–p show the same with m = 2. Parameter values used are shown in the figure. Inefficient fertilisation combined with relatively low asymmetry in gamete numbers and the added asymmetry of internal fertilisation can in principle reverse the Bateman gradients (second and fourth row). Females (gamete number nx) are indicated by blue crosses and connecting lines, while males (gamete number ny) are indicated by black dots and connecting lines.Full size image More