Ecology

Nesting recordsDipogon nests were created singly per cane, because there were no examples in which wasps of two species emerged from the same cane in the study site. Thus, we designate “utilized canes” as “nests”.In the four years, in pine forests in Takarazuka, Hyogo, Japan, we collected a total of 419 nests with 1033 cells from which species of Dipogon emerged (Fig. 1; Table 1; Supplementary Table S1). The numbers of nests and cells and the average and SD of the number of cells per nest for each species are shown in Table 1. Other wasps, bees, and parasitic wasps and flies also emerged from our trap nests (Supplementary Table S2), but we did not consider their nesting in the following analyses. Among 1033 cells, D. sperconsus emerged from 623 cells, D. inconspersus from 26 cells, and D. bifasciatus (Geoffroy) from 4 cells, while rearing failure occurred in 380 cells (Table 1), the owners of which we designate as “unknown Dipogon spp.” Based on the total cells of Dipogon, the proportion of cells constructed by D. sperconsus was 60.3% (623/1033*100), that of D. inconspersus was 2.5% (26/1033*100), and that of D. bifasciatus was 0.39% (4/1033*100). Based on the cells of the identified species, the proportion of cells constructed by D. sperconsus was 95.4% (623/(623 + 26 + 4)*100), that of D. inconspersus was 4.0% (26/(623 + 26 + 4)*100), and that of D. bifasciatus was 0.6% (4/(623 + 26 + 4)*100). From these proportions, we can estimate the number of cells constructed by the three species of Dipogon in the total 1033 Dipogon cells as ca. 985.5 cells (1033*0.954) by D. sperconsus, ca. 41.3 cells (1033*0.04) by D inconspersus, and ca. 6.2 cells (1033*0.006) by D. bifasciatus.Figure 1The study site in Kirihata, Takarazuka City, Hyogo Pref., Japan, and trap nests. (a) An old pine forest in which trap nests were installed. (b) A set of trap nests (cane bundle), 15 mixed-size bamboo canes bound vertically with vinyl-covered wires like a screen, attached to a tree trunk approximately 1.5 m above the ground. (c) A nest of D. sperconsus; this cane was installed in Shibutani, Takarazuka, Hyoto Pref. about 1 km southeast of the present study site on 29 July 2007 and was withdrawn on 6 August 2007. (d) A nest (6–5-5–1) of D. sperconsus; this cane was installed in Kirihata, Takarazuka, Hyoto Pref. about 500 m west-southwest of the present study site on 25 August 2010 and was withdrawn on 27 August 2010 (prey spider, Agelena limbata Thorell). (e) A nest of D. sperconsus; this cane was installed in Najio, Nishinomiya, Hyoto Pref. about 10 km southwest of the present study site on 15 July 2007 and was withdrawn on 25 July 2007. The minimum grid in the background graph paper of (c)–(e) is 1 mm. All photos taken by Y. Nishimoto.Full size imageTable 1 The numbers of the collected nests and brood cells, and the mean number of cells per nest in three species of Dipogon (Deuteragenia).Full size tableBecause multiple cells were often constructed in a single nest, the number of nests was much smaller than the number of constructed cells. Among the 419 nests, 221 nests belonged to D. sperconsus, 7 nests belonged to D. inconspersus, and a single nest belonged to D. bifasciatus, but the remaining 190 nests could not be identified because of rearing failure (Table 1). The proportions of the nests in the three Dipogon species were calculated as follows: 96.5% (221/(221 + 7 + 1)*100) in D. sperconsus, 3.1% (7/(221 + 7 + 1)*100) in D inconspersus, and 0.4% (1/(221 + 7 + 1)*100) in D. bifasciatus. Thus, the estimated number of nests in each species was ca. 404.3 (419*0.965) in D. sperconsus, ca. 13.0 (419*0.031) in D inconspersus, and ca. 1.7 (419*0.004) in D. bifasciatus.Next, we considered whether the cane bundles were used randomly. Based on the yearly frequency distributions of nests (Supplementary Tables S3–S6), we developed a null hypothesis assuming the nests are randomly distributed over bundles, where a negative binomial distribution is expected (Supplementary Tables S7–S8). Our yearly data indicate that the null hypothesis was rejected and that nests were more or less aggregated in a few bundles (Supplementary Figure S1; test statistics, Supplementary Table S8). This aggregation tendency (e.g., no nests in some bundles) may imply that some selected sites for bundles are not appropriate for D. sperconsus, for some unknown behavioral reasons. Further studies are needed to verify the habitat use of this species.Yearly frequency distributions of the number of cells show that the range of cells constructed by D. sperconsus and unknown D. spp. combined were 1–10 cells, and the median was 2 cells (Supplementary Table S3–S6, Supplementary Figure S2). Most of the nests included 1–3 cells, and five or more cells were very rare. Most of the nests with many cells (e.g., 7–10 cells) were likely to be constructed by a single wasp because these wasps avoid interactions with other spider wasps. The average number of D. sperconsus cells per nest was 2.82 for four years, varying from 2.21 (2014) to 3.16 (2016) (Table 1), and the yearly differences were significant (Kruskal–Wallis test, (chi ^{2} = 7.70), df = 3, p = 0.05). In contrast, the average number of cells per nest of D. sayi sayi was slightly greater than that of D. sperconsus: 3.2 (1–6, SD = 1.47, n = 41) in the first generation and 4.7 (1–13, SD = 2.52, n = 107) in the second generation in Wisconsin, USA8; and 6.2 (1961), 4.0 (1962) and 3.0 (1963) in the summer generation and 7.5 (1961) and 3.2 (1962) in the overwintering generation in Northwestern Ontario9.Life history of Dipogon sperconsusDevelopmental periodThe developmental period of reared wasps was estimated in the summer and overwintering generations separately (Table 2, Supplementary Figure S3, Supplementary Tables S9–S12). In the summer generation, both females and males developed from egg to adult over approximately three weeks (23.1 days for females and 21.6 days for males; Table 2). There was no significant difference between sexes (t-test, after adjustment by Bonferroni method: p  > 0.05). In the overwintering generation, approximately eight months were required from egg to adult (246 days for females and 247 days for males). There was also no significant difference between sexes (t-test, after adjustment by Bonferroni method: p  > 0.05). In females, all developmental periods were significantly longer in the overwintering generation than in the summer generation (t-test, after adjustment by Bonferroni method: p  0.05 for egg and larval periods; p  0.1). Among the 40 coelotid female spiders, the sex ratio of wasp eggs was even: 20 female wasp eggs and 20 males. However, the female spiders on which female wasp eggs were laid were significantly greater in cephalothorax width than those on which male eggs were laid (t = 3.98, p  More

1.Schoener, T. W. Field experiments on interspecific competition. Am. Nat. 122, 240–285 (1983).Article 

Google Scholar 
2.Keddy, P. A. Competition 2nd edn. (Kluwer, 2001).Book 

Google Scholar 
3.Kotler, B. P. & Brown, J. S. Environmental heterogeneity and the coexistence of desert rodents. Annu. Rev. Ecol. Syst. 19, 281–307 (1988).Article 

Google Scholar 
4.Kronfeld-Schor, N. & Dayan, T. Partitioning of time as an ecological resource. Annu. Rev. Ecol. Evol. Syst. 34, 153–181 (2003).Article 

Google Scholar 
5.Connell, J. H. On the prevalence and relative importance of interspecific competition: evidence from field experiments. Am. Nat. 122, 661–696 (1983).Article 

Google Scholar 
6.Adler, P. B. et al. Competition and coexistence in plant communities: intraspecific competition is stronger than interspecific competition. Ecol. Lett. 21, 1319–1329 (2018).PubMed 
Article 
PubMed Central 

Google Scholar 
7.Morris, D. W. Toward an ecological synthesis: a case for habitat selection. Oecologia 136, 1–13 (2003).ADS 
PubMed 
Article 
PubMed Central 

Google Scholar 
8.Barkae, E. D., Abramsky, Z. & Ovadia, O. Can models of density-dependent habitat selection be applied for trap-building predators?. Popul. Ecol. 56, 175–184 (2014).Article 

Google Scholar 
9.Halliday, W. D. & Blouin-Demers, G. Red flour beetles balance thermoregulation and food acquisition via density-dependent habitat selection. J. Zool. 294, 198–205 (2014).Article 

Google Scholar 
10.Tregenza, T. Building on the ideal free distribution. Adv. Ecol. Res. 26, 253–307 (1995).Article 

Google Scholar 
11.Kingsolver, J. G. & Pfennig, D. W. Individual-level selection as a cause of Cope’s rule of phyletic size increase. Evolution 58, 1608–1612 (2004).PubMed 
Article 
PubMed Central 

Google Scholar 
12.Alatalo, R. V. & Moreno, J. Body size, interspecific interactions, and use of foraging sites in tits (Paridae). Ecology 68, 1773–1777 (1987).PubMed 
Article 
PubMed Central 

Google Scholar 
13.Honěk, A. Intraspecific variation in body size and fecundity in insects: a general relationship. Oikos 66, 483–492 (1993).Article 

Google Scholar 
14.Sokolovska, N., Rowe, L. & Johansson, F. Fitness and body size in mature odonates. Ecol. Entomol. 25, 239–248 (2000).Article 

Google Scholar 
15.Werner, E. E. & Anholt, B. R. Ecological consequences of the trade-off between growth and mortality rates mediated by foraging activity. Am. Nat. 142, 242–272 (1993).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
16.Blanckenhorn, W. U. The evolution of body size: What keeps organisms small?. Q. Rev. Biol. 75, 385–407 (2000).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
17.Gotthard, K. Increased risk of predation as a cost of high growth rate: an experimental test in a butterfly. J. Anim. Ecol. 69, 896–902 (2000).PubMed 
Article 
PubMed Central 

Google Scholar 
18.Van Buskirk, J. Competition, cannibalism, and size class dominance in a dragonfly. Oikos 65, 455–464 (1992).Article 

Google Scholar 
19.Fincke, O. M. Larval behaviour of a giant damselfly: Territoriality or size-dependent dominance?. Anim. Behav. 51, 77–87 (1996).Article 

Google Scholar 
20.Hopper, K. R., Crowley, P. H. & Kielman, D. Density dependence, hatching synchrony, and within-cohort cannibalism in young dragonfly larvae. Ecology 77, 191–200 (1996).Article 

Google Scholar 
21.Eitam, A., Blaustein, L. & Mangel, M. Density and intercohort priority effects on larval Salamandra salamandra in temporary pools. Oecologia 146, 36–42 (2005).ADS 
PubMed 
Article 
PubMed Central 

Google Scholar 
22.Barkae, E. D., Scharf, I. & Ovadia, O. A stranger is tastier than a neighbor: cannibalism in Mediterranean and desert antlion populations. Behav. Ecol. 28, 69–76 (2017).Article 

Google Scholar 
23.Alford, R. A. & Wilbur, H. M. Priority effects in experimental pond communities: competition between Bufo and Rana. Ecology 66, 1097–1105 (1985).Article 

Google Scholar 
24.Dayton, G. H. & Fitzgerald, L. A. Priority effects and desert anuran communities. Can. J. Zool. 83, 1112–1116 (2005).Article 

Google Scholar 
25.Louette, G. & De Meester, L. Predation and priority effects in experimental zooplankton communities. Oikos 116, 419–426 (2007).Article 

Google Scholar 
26.Geange, S. W. & Stier, A. C. Order of arrival affects competition in two reef fishes. Ecology 90, 2868–2878 (2009).PubMed 
Article 
PubMed Central 

Google Scholar 
27.Huey, R. B. & Pianka, E. R. Ecological consequences of foraging mode. Ecology 62, 991–999 (1981).Article 

Google Scholar 
28.Shine, R. & Li-Xin, S. Arboreal ambush site selection by pit-vipers Gloydius shedaoensis. Anim. Behav. 63, 565–576 (2002).Article 

Google Scholar 
29.Clark, R. W. Feeding experience modifies the assessment of ambush sites by the timber rattlesnake, a sit-and-wait predator. Ethology 110, 471–483 (2004).Article 

Google Scholar 
30.Tsairi, H. & Bouskila, A. Ambush site selection of a desert snake (Echis coloratus) at an oasis. Herpetologica 60, 13–23 (2004).Article 

Google Scholar 
31.Scharf, I., Lubin, Y. & Ovadia, O. Foraging decisions and behavioural flexibility in trap-building predators: a review. Biol. Rev. 86, 626–639 (2011).PubMed 
Article 

Google Scholar 
32.Blamires, S. J. Biomechanical costs and benefits of sit-and-wait foraging traps. Isr. J. Ecol. Evol. 66, 5–14 (2020).Article 

Google Scholar 
33.Simberloff, D. et al. Holes in the doughnut theory: the dispersion of ant-lions. Brenesia 14, 13–46 (1978).
Google Scholar 
34.Farji-Brener, A. G., Carvajal, D., Gei, M. G., Olano, J. & Sanchez, J. D. Direct and indirect effect of soil structure on the density of an antlion larva in a tropical dry forest. Ecol. Entomol. 33, 183–188 (2008).Article 

Google Scholar 
35.Lucas, J. R. Metabolic rates and pit-construction costs of two antlion species. J. Anim. Ecol. 54, 295–309 (1985).Article 

Google Scholar 
36.Tanaka, K. Energetic cost of web construction and its effect on web relocation in the web-building spider Agelena limbata. Oecologia 81, 459–464 (1989).ADS 
PubMed 
Article 
PubMed Central 

Google Scholar 
37.Lubin, Y., Ellner, S. & Kotzman, M. Web relocation and habitat selection in desert widow spider. Ecology 74, 1915–1928 (1993).Article 

Google Scholar 
38.Loria, R., Scharf, I., Subach, A. & Ovadia, O. The interplay between foraging mode, habitat structure, and predator presence in antlions. Behav. Ecol. Sociobiol. 62, 1185–1192 (2008).Article 

Google Scholar 
39.Griffiths, D. Interference competition in ant-lion (Macroleon quinquemaculatus) larvae. Ecol. Entomol. 17, 219–226 (1992).Article 

Google Scholar 
40.Heiling, A. M. & Herberstein, M. E. The importance of being larger: intraspecific competition for prime web sites in orb-web spiders (Araneae, Araneidae). Behaviour 136, 669–677 (1999).Article 

Google Scholar 
41.Rayor, L. S. & Uetz, G. W. Trade-offs in foraging success and predation risk with spatial position in colonial spiders. Behav. Ecol. Sociobiol. 27, 77–85 (1990).Article 

Google Scholar 
42.Wilson, D. S. Prey capture and competition in the ant lion. Biotropica 6, 187–193 (1974).Article 

Google Scholar 
43.Rao, D. Experimental evidence for the amelioration of shadow competition in an orb-web spider through the ‘ricochet’ effect. Ethology 115, 691–697 (2009).Article 

Google Scholar 
44.Scharf, I. Factors that can affect the spatial positioning of large and small individuals in clusters of sit-and-wait predators. Am. Nat. 195, 649–663 (2020).PubMed 
Article 
PubMed Central 

Google Scholar 
45.Matsura, T. & Takano, H. Pit-relocation of antlion larvae in relation to their density. Res. Popul. Ecol. 31, 225–234 (1989).Article 

Google Scholar 
46.Griffiths, D. Intraspecific competition in larvae of the ant-lion Morter sp. and interspecific interactions with Macroleon quinquemaculatus. Ecol. Entomol. 16, 193–201 (1991).Article 

Google Scholar 
47.Wise, D. H. Cannibalism, food limitation, intraspecific competition, and the regulation of spider populations. Annu. Rev. Entomol. 51, 441–465 (2006).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
48.Klokočovnik, V., Veler, E. & Devetak, D. Antlions in interaction: confrontation of two competitors in limited space. Isr. J. Ecol. Evol. 66, 73–81 (2020).Article 

Google Scholar 
49.Buddle, C. M., Walker, S. E. & Rypstra, A. L. Cannibalism and density-dependent mortality in the wolf spider Pardosa milvina (Araneae: Lycosidae). Can. J. Zool. 81, 1293–1297 (2003).Article 

Google Scholar 
50.Ovadia, O., Scharf, I., Barkae, E. D., Levi, T. & Alcalay, Y. Asymmetrical intra-guild predation and niche differentiation in two pit-building antlions. Isr. J. Ecol. Evol. 66, 82–90 (2020).Article 

Google Scholar 
51.Devetak, D. Wormlion Vermileo vermileo (L.) (Diptera: Vermileonidae) in Slovenia and Croatia. Ann. Ser. Hist. Nat. 18, 283–286 (2008).
Google Scholar 
52.Dor, R., Rosenstein, S. & Scharf, I. Foraging behaviour of a neglected pit-building predator: the wormlion. Anim. Behav. 93, 69–76 (2014).Article 

Google Scholar 
53.Miler, K., Yahya, B. E. & Czarnoleski, M. Substrate moisture, particle size and temperature preferences of trap-building larvae of sympatric antlions and wormlions from the rainforest of Borneo. Ecol. Entomol. 44, 488–493 (2019).Article 

Google Scholar 
54.Miler, K., Yahya, B. E. & Czarnoleski, M. Different predation efficiencies of trap-building larvae of sympatric antlions and wormlions from the rainforest of Borneo. Ecol. Entomol. 43, 255–262 (2018).Article 

Google Scholar 
55.Franks, N. R., Worley, A., Falkenberg, M., Sendova-Franks, A. B. & Christensen, K. Digging the optimum pit: antlions, spirals and spontaneous stratification. Proc. R. Soc. B 286, 20190365 (2019).PubMed 
Article 
PubMed Central 

Google Scholar 
56.Scharf, I., Daniel, A., MacMillan, H. A. & Katz, N. The effect of fasting and body reserves on cold tolerance in 2 pit-building insect predators. Curr. Zool. 63, 287–294 (2017).PubMed 
PubMed Central 

Google Scholar 
57.Devetak, D. Substrate particle size-preference of wormlion Vermileo vermileo (Diptera: Vermileonidae) larvae and their interaction with antlions. Eur. J. Entomol. 105, 631–635 (2008).Article 

Google Scholar 
58.Adar, S., Dor, R. & Scharf, I. Habitat choice and complex decision making in a trap-building predator. Behav. Ecol. 27, 1491–1498 (2016).Article 

Google Scholar 
59.Scharf, I. et al. The contribution of shelter from rain to the success of pit-building predators in urban habitats. Anim. Behav. 142, 139–145 (2018).Article 

Google Scholar 
60.Katz, N., Pruitt, J. N. & Scharf, I. The complex effect of illumination, temperature, and thermal acclimation on habitat choice and foraging behavior of a pit-building wormlion. Behav. Ecol. Sociobiol. 71, 137 (2017).Article 

Google Scholar 
61.Bar-Ziv, M. A., Bega, D., Subach, A. & Scharf, I. Wormlions prefer both fine and deep sand but only deep sand leads to better performance. Curr. Zool. 65, 393–400 (2019).PubMed 
Article 
PubMed Central 

Google Scholar 
62.Abramoff, M. D., Magalhaes, P. J. & Ram, S. J. Image processing with ImageJ. Biophoton. Int. 11, 36–42 (2004).
Google Scholar 
63.Ovadia, O. & Abramsky, Z. Density-dependent habitat selection: evaluation of the isodar method. Oikos 73, 86–94 (1995).Article 

Google Scholar 
64.Jensen, W. E. & Cully, J. F. Density-dependent habitat selection by brown-headed cowbirds (Molothrus ater) in tallgrass prairie. Oecologia 142, 136–149 (2005).ADS 
PubMed 
Article 
PubMed Central 

Google Scholar 
65.Whitham, T. G. The theory of habitat selection: examined and extended using Pemphigus aphids. Am. Nat. 115, 449–466 (1980).Article 

Google Scholar 
66.van Beest, F. M. et al. Increasing density leads to generalization in both coarse-grained habitat selection and fine-grained resource selection in a large mammal. J. Anim. Ecol. 83, 147–156 (2014).PubMed 
Article 
PubMed Central 

Google Scholar 
67.Mathis, A. Territoriality in a terrestrial salamander: the influence of resource quality and body size. Behaviour 112, 162–175 (1990).Article 

Google Scholar 
68.Croy, M. I. & Hughes, R. N. Effects of food supply, hunger, danger and competition on choice of foraging location by the fifteen-spined stickleback, Spinachia spinachia L. Anim. Behav. 42, 131–139 (1991).Article 

Google Scholar 
69.Davey, A. J. H., Hawkins, S. J., Turner, G. F. & Doncaster, C. P. Size-dependent microhabitat use and intraspecific competition in Cottus gobio. J. Fish Biol. 67, 428–443 (2005).Article 

Google Scholar 
70.Abrahams, M. V. Patch choice under perceptual constraints: a cause for departures from an ideal free distribution. Behav. Ecol. Sociobiol. 19, 409–415 (1986).Article 

Google Scholar 
71.Sutherland, W. J., Townsend, C. R. & Patmore, J. M. A test of the ideal free distribution with unequal competitors. Behav. Ecol. Sociobiol. 23, 51–53 (1988).Article 

Google Scholar 
72.McClure, M. S. Spatial distribution of pit-making ant lion larvae (Neuroptera: Myrmeleontidae): density effects. Biotropica 8, 179–183 (1976).Article 

Google Scholar 
73.Rayor, L. S. & Uetz, G. W. Age-related sequential web building in the colonial spider Metepeira incrassata (Araneidae): an adaptive spacing strategy. Anim. Behav. 59, 1251–1259 (2000).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
74.Yip, E. C., Levy, T. & Lubin, Y. Bad neighbors: hunger and dominance drive spacing and position in an orb-weaving spider colony. Behav. Ecol. Sociobiol. 71, 128 (2017).Article 

Google Scholar 
75.Murcia, C. Edge effects in fragmented forests: implications for conservation. Trends Ecol. Evol. 10, 58–62 (1995).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
76.Minias, P., Janiszewski, T. & Lesner, B. Center-periphery gradients of chick survival in the colonies of Whiskered Terns Chlidonias hybrida may be explained by the variation in the maternal effects of egg size. Acta Ornithol. 48, 179–186 (2013).Article 

Google Scholar 
77.Geange, S. W. & Stier, A. C. Priority effects and habitat complexity affect the strength of competition. Oecologia 163, 111–118 (2010).ADS 
PubMed 
Article 
PubMed Central 

Google Scholar 
78.Hallander, H. Prey, cannibalism and microhabitat selection in the wolf spiders Pardosa chelata OF Müller and P. pullata Clerck. Oikos 21, 337–340 (1970).Article 

Google Scholar 
79.Skevington, J. H. & Dang, P. T. Exploring the diversity of flies (Diptera). Biodiversity 3, 3–27 (2002).Article 

Google Scholar 
80.Scharf, I., Silberklang, A., Avidov, B. & Subach, A. Do pit-building predators prefer or avoid barriers? Wormlions’ preference for walls depends on light conditions. Sci. Rep. 10, 10928 (2020).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar  More

Russia has been reporting almost no changes in forested area, growing stock volume (GSV) and biomass to the United Nations Framework Convention on Climate Change (UNFCCC)1 and the Food and Agriculture Organization of the United Nations (FAO) Forest Resources Assessment (FRA)2 since the collapse of the USSR and the decline in the Soviet Forest Inventory and Planning (FIP) system. According to the State Forest Register (SFR)3, which is the main repository of forest information, and national reporting to the FAO FRA2, the GSV and the above ground biomass (AGB) increased by 1.1% and 0.6% (Table S1), respectively, during 1990–2015, yet studies using remote sensing (RS) indicate increased vegetation productivity4, tree cover (annual rate + 0.417% over 1982–2016)5, increased AGB (+ 329 Tg C yr−1 over 2000–20076), total biomass (annual rate + 0.44% or + 153 Tg C yr−1 over 1990–20077), and forest ecosystem carbon pools (ca + 470 Tg C yr−1 over 2001–20198). This inconsistency in estimates can be explained by an information gap that appeared when Russia decided to move from the FIP to another system for the collection of forest information at the national scale – the National Forest Inventory (NFI).The FIP involves revisiting every forest stand (on the ground for managed forests or using RS techniques for remote non-commercial forests) on a 10–15-year interval, with the measurement of forest parameters combined with the formulation of forest management directives. After the collapse of the USSR, the inventory within the FIP system slowed down substantially. For example, more than 50% of the forest area was surveyed by the FIP more than 25 years ago9. For these reasons, the reliability of information on forests in Russia has deteriorated since 1988, which is the year when FIP-based reporting10 involved the largest inventory efforts in recent decades. According to this report10, the total GSV of Russian forests was 81.7 × 109 m3 (without shrubland, bias corrected11). This value is used here as a reference to quantify biomass stock changes in Russia with respect to the current decade.In contrast, NFI is a state-of-the-art inventory system based on a statistical sampling method. It was initiated in 2007 and the first cycle was completed in 2020. The NFI data processing is ongoing, but the first official press release12 suggests that Russian forest accumulated 102 × 109 m3 over its lifespan until 2014. Once finalized, the NFI will be verified before adoption as the official source of information to the SFR and national reporting. The NFI has received some criticism13 because of the relatively sparse sampling employed and the stratification method used, which is partially based on outdated FIP data.In Russia, the long intervals between consecutive surveys and the difficulty in accessing very remote regions in a timely manner by an inventory system make satellite RS an essential tool for capturing forest dynamics and providing a comprehensive, wall-to-wall perspective on biomass distribution. However, observations from current RS sensors are not suited for producing accurate biomass estimates unless the estimation method is calibrated with a dense network of measurements from ground surveys14. Here we calibrated models relating two global RS biomass data products (GlobBiomass GSV15 and CCI Biomass GSV16) and additional RS data layers (forest cover mask9, the Copernicus Global Land Cover CGLS‐LC100 product17) with ca 10,000 ground plots (see Material and Methods) to reduce nuances in the individual input maps due to imperfections in the RS data and approximations in the retrieval procedure18,19. The combination of these two sources of information, i.e., ground measurements and RS, utilizes the advantages of both sources in terms of: (i) highly accurate ground measurements and (ii) the spatially comprehensive coverage of RS products and methods. The amount of ground plots currently available may be insufficient for providing an accurate estimate of GSV for the country when used alone, but they are the key to obtaining unbiased estimates when used to calibrate RS datasets20. The map merging procedure was preferred over a plot-aided direct estimation of GSV or AGB from the RS data because of the usually poor association between biomass measured at inventory plots and remote sensing observables21. In addition, models relating biomass and remote sensing observables that are trained with spatially inhomogeneous datasets (Figure S1) tend to be biased in regions not represented by the dataset of the reference biomass measurements.We estimate the total GSV of Russia for the year 2014 for the official forested area (713.1 × 106 ha) to be 111 ± 1.3 × 109 m3, which is 39% higher than the 79.9 × 109 m3 (excluding shrubland) figure reported in the SFR3 for the same year. An additional 7.1 × 109 m3 or 9% were found due to the larger forested area (+ 45.7 106 ha) recognized by RS9, following the expansion of forests to the north22, to higher elevations, in abandoned arable land23, as well as the inclusion of parks, gardens and other trees outside of forest, which were not counted as forest in the SFR. Based on cross-validation, our estimate at the regional level (81 regions of Russia – Table S2, Figure S2) is unbiased. The standard error varied from 0.6 to 17.6% depending on the region. The median error was 1.6%, while the area weighted error was 1.2%. The predicted GSV (Fig. 1) with associated uncertainties is available here (https://doi.org/10.5281/zenodo.3981198) as a GeoTiff at a spatial resolution of 3.2 arc sec. (ca 0.5 ha).Figure 1Predicted mean forest growing stock volume (m3 ha-1) for the year ca 2014 (Generated by Esri ArcGIS Desktop v.10.7, URL: https://desktop.arcgis.com/en/arcmap/).Full size imageHoughton et al.24 estimated forest biomass based on RS and FIP data in Russia for the year 2000. Average forest biomass density varied between 80.6 and 88.2 Mg ha-1 depending on which forest mask was used. Our estimate for the year 2014 of 107 Mg ha-1 (using the conversion factor of GSV to AGB from24 0.6859) is 21–33% higher than the one by Houghton et al., but this is consistent with expected biomass increases over time, i.e., 14 years after the Houghton et al. estimate.Assuming an unchanged total forest area (721.7 × 106 ha) in 1988 and 2014, we conclude that Russian forests have accumulated 1,163 × 106 m3 yr-1 or 407 Tg C yr-1 in live biomass of trees on average over 26 years. This gives an average GSV change rate of + 1.61 m3 ha-1 yr-1 or + 0.56 t C ha-1 yr-1. The sequestration rate obtained, however, should be treated with caution because different methods have been applied in 1988 and 2014 (see “Caveats and Limitations” section). To provide some context for the magnitude of these numbers, one can compare the Russian forest gain to the net GSV losses in tropical forests over the period 1990–2015 according to FAO FRA25 (-1,033 × 106 m3 yr-1 in the regions with a negative trend: South and Central America, South and Southeast Asia, and Africa). A similar divergence in the carbon sink between Tropical and Boreal forest was recognized by Tagesson et al.26.In terms of carbon stock change, our estimates are substantially higher than those reported by Pan et al.7 for 1990–2007 (+ 153 Tg C yr-1) based on FIP data. The biomass carbon estimates by Liu et al.6 are instead in line with our results. There is an increase in the annual rate of AGB in Russia of + 329 Tg C yr−1 (annual variation from 214 to 400 Tg C yr−1) over 2000–20076. Interestingly, another boreal country – Canada – has demonstrated neutral or negative trends (from 0 to -14 Tg C yr−1) for the same time span using the same estimation method6.We can observe different spatial patterns in the change in the GSV density between 1988 (FIP10, bias corrected11) and 2014 (our estimate), which can be explained by climate change, CO2 fertilisation and changes in disturbance regimes (Fig. 2). The average linear trend in the annual temperature increase during 1976–2014 in Russia is + 0.45 °C per 10 years27. The temperature increase is statistically significant in every region except for western Siberia (Fig. 2–3). Significantly increased temperature extremes and an increase in the number of days without precipitation is observed in the south of European Russia, Baikal, Kamchatka, and Chukotka27 (Fig. 2–1). Some regions in the south of the European part of Russia are colored in dark blue, but they, as a rule, have a small share of forested area, which is often linked to water bodies and, therefore, suffers less from increased drought (Fig. 2–1). Central and eastern Siberia suffer from an increase in disturbances, which offsets the climate stimulation effect (Fig. 2–4). The forested area in the Nenets region (Fig. 2–2) is 4 times larger in 2014 based on the RS forest mask compared to the SFR in 1988 (where forest was accounted for up until a certain latitude at that time), where the increase in area resulted in a decrease in the average GSV.Figure 2Change in growing stock volume (m3 ha-1) from 1988 to 2014 (average over administrative regions) (Generated by Esri ArcGIS Desktop v.10.7, URL: https://desktop.arcgis.com/en/arcmap/). These changes can be categorized into: 1—significant increase in air temperature and drought; 2—substantially increased forest area, which lowers the average GSV density; 3—least (not significant) temperature increase; 4—increase of disturbances: wildfire and harvest (southern part), which offsets the climate stimulation effect.Full size imageFocusing specifically on national reporting of managed forest to the UNFCCC, 72% of forested area in Russia is considered to be managed1. We multiplied the GSV density by the managed forest area for each administrative region (Table S3). The difference in GSV estimation (between ours and the one from the SFR report) is 23.6 × 109 m3 (Table S3) or 33% higher. From the GSV of managed forests in 2014 and based on the same area in 1988, we can estimate the sequestration rate of live biomass of managed forests as 354 Tg C yr-1 , which is considerably higher than the figure of 230 Tg C yr-1 in the current report1.This proof of concept demonstrates the relevance of complementing recent NFI data with remote sensing map products. Our study demonstrates that the already considerable value of forest inventory data can be further enhanced in a forest resources mapping scenario. In addition, we seek to promote greater access to these data by opening up their access to the larger scientific community. Through the integration of RS estimates of GSV and forest inventory data from Russia, we confirm that carbon stocks increased substantially during the last few decades in contrast to the figures provided in official national reporting. Russian forests play an even more important global role in carbon sequestration than previously thought, where the increase in growing stock is of the same magnitude as the net losses in tropical forests over the same time period. More

1.Falconer, D. S. & Mackay, T. F. C. Introduction to Quantitative Genetics (Longman, 1996).2.Lande, R. & Arnold, S. J. The measurement of selection on correlated characters. Evolution 37, 1210–1226 (1983).Article 

Google Scholar 
3.Kingsolver, J. G., Diamond, S. E., Siepielski, A. M. & Carlson, S. M. Synthetic analyses of phenotypic selection in natural populations: lessons, limitations and future directions. Evol. Ecol. 26, 1101–1118 (2012).Article 

Google Scholar 
4.Barrett, R. D. H. & Schluter, D. Adaptation from standing genetic variation. Trends Ecol. Evol. 23, 38–44 (2008).PubMed 
Article 
PubMed Central 

Google Scholar 
5.Kulbaba, M. W., Sheth, S. N., Pain, R. E., Eckhart, V. M. & Shaw, R. G. Additive genetic variance for lifetime fitness and the capacity for adaptation in an annual plant. Evolution 73, 1746–1758 (2019).PubMed 
Article 
PubMed Central 

Google Scholar 
6.Lande, R. & Shannon, S. The role of genetic variation in adaptation and population persistence in a changing environment. Evolution 50, 434–437 (1996).PubMed 
Article 
PubMed Central 

Google Scholar 
7.Etterson, J. R. & Shaw, R. G. Constraint to adaptive evolution in response to global warming. Science 294, 151–154 (2001).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
8.Anderson, J. T., Inouye, D. W., McKinney, A. M., Colautti, R. I. & Mitchell-Olds, T. Phenotypic plasticity and adaptive evolution contribute to advancing flowering phenology in response to climate change. Proc. R. Soc. B 279, 3843–3852 (2012).PubMed 
Article 
PubMed Central 

Google Scholar 
9.Steffen, W., Crutzen, P. J. & McNeil, J. R. The Anthropocene: are humans now overwhelming the great forces of nature? Ambio 36, 614–621 (2007).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
10.Zhang, X.-S. & Hill, W. G. Genetic variability under mutation selection balance. Trends Ecol. Evol. 20, 468–470 (2005).PubMed 
Article 
PubMed Central 

Google Scholar 
11.McGuigan, K., Aguirre, J. D. & Blows, M. W. Simultaneous estimation of additive and mutational genetic variance in an outbred population of Drosophila serrata. Genetics 201, 1239–1251 (2015).PubMed 
PubMed Central 
Article 

Google Scholar 
12.Huang, W. et al. Spontaneous mutations and the origin and maintenance of quantitative genetic variation. eLife 5, e14625 (2016).PubMed 
PubMed Central 
Article 

Google Scholar 
13.Mitchell-Olds, T., Willis, J. H. & Goldstein, D. B. Which evolutionary processes influence natural genetic variation for phenotypic traits? Nat. Rev. Genet. 8, 845–856 (2007).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
14.Charlesworth, B. Causes of natural variation in fitness: evidence from studies of Drosophila populations. Proc. Natl Acad. Sci. USA 112, 1662–1669 (2015).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
15.Subramaniam, B. & Rausher, M. D. Balancing selection on a floral polymorphism. Evolution 54, 691–695 (2000).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
16.Charlesworth, D. Balancing selection and its effects on sequences in nearby genome regions. PLoS Genet. 2, e64 (2006).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
17.Hedrick, P. W. & Thomson, G. Evidence for balancing selection at HLA. Genetics 104, 449–456 (1983).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
18.Troth, A., Puzey, J. R., Kim, R. S., Willis, J. H. & Kelly, J. K. Selective trade-offs maintain alleles underpinning complex trait variation in plants. Science 361, 475–478 (2018).CAS 
PubMed 
PubMed Central 

Google Scholar 
19.Delph, L. F. & Kelly, J. K. On the importance of balancing selection in plants. N. Phytol. 201, 45–56 (2014).Article 

Google Scholar 
20.Anderson, J. T., Wagner, M. R., Rushworth, C. A., Prasad, K. V. S. K. & Mitchell-Olds, T. The evolution of quantitative traits in complex environments. Heredity 112, 4–12 (2014).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
21.Anderson, J. T. & Wadgymar, S. M. Climate change disrupts local adaptation and favours upslope migration. Ecol. Lett. 23, 181–192 (2020).PubMed 
Article 
PubMed Central 

Google Scholar 
22.Agrawal, A. A. & Fishbein, M. Plant defense syndromes. Ecology 87, S132–S149 (2006).PubMed 
Article 
PubMed Central 

Google Scholar 
23.Carmona, D., Lajeunesse, M. J. & Johnson, M. T. Plant traits that predict resistance to herbivores. Funct. Ecol. 25, 358–367 (2011).Article 

Google Scholar 
24.DeLucia, E. H., Nabity, P. D., Zavala, J. A. & Berenbaum, M. R. Climate change: resetting plant–insect interactions. Plant Physiol. 160, 1677–1685 (2012).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
25.Mithöfer, A. & Boland, W. Plant defense against herbivores: chemical aspects. Annu. Rev. Plant Biol. 63, 431–450 (2012).PubMed 
Article 
CAS 
PubMed Central 

Google Scholar 
26.Prasad, K. V. S. K. et al. A gain-of-function polymorphism controlling complex traits and fitness in nature. Science 337, 1081–1084 (2012).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
27.Bergelson, J., Dwyer, G. & Emerson, J. J. Models and data on plant–enemy coevolution. Annu. Rev. Genet. 35, 469–499 (2001).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
28.Hodgins, K. A. & Barrett, S. C. H. Female reproductive success and the evolution of mating-type frequencies in tristylous populations. N. Phytol. 171, 569–580 (2006).Article 

Google Scholar 
29.Trotter, M. V. & Spencer, H. G. Complex dynamics occur in a single-locus, multiallelic model of general frequency-dependent selection. Theor. Popul. Biol. 76, 292–298 (2009).PubMed 
Article 
PubMed Central 

Google Scholar 
30.Tuinstra, M. R., Ejeta, G. & Goldsbrough, P. B. Heterogeneous inbred family (HIF) analysis: a method for developing near-isogenic loci that differ at quantitative traits. Theor. Appl. Genet. 95, 1005–1011 (1997).CAS 
Article 

Google Scholar 
31.Salehin, M. et al. Auxin-sensitive Aux/IAA proteins mediate drought tolerance in Arabidopsis by regulating glucosinolate levels. Nat. Commun. 10, 4021 (2019).PubMed 
PubMed Central 
Article 

Google Scholar 
32.Hossain, M. S. et al. Glucosinolate degradation products, isothiocyanates, nitriles, and thiocyanates, induce stomatal closure accompanied by peroxidase-mediated reactive oxygen species production in Arabidopsis thaliana. Biosci. Biotechnol. Biochem. 77, 977–983 (2013).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
33.Mitchell-Olds, T. & Schmitt, J. Genetic mechanisms and evolutionary significance of natural variation in Arabidopsis. Nature 441, 947–952 (2006).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
34.Wang, B. et al. Ancient polymorphisms contribute to genome-wide variation by long-term balancing selection and divergent sorting in Boechera stricta. Genome Biol. 20, 126 (2019).PubMed 
PubMed Central 
Article 

Google Scholar 
35.Bloom, T. C., Baskin, J. M. & Baskin, C. C. Ecological life history of the facultative woodland biennial Arabis laevigata variety laevigata (Brassicaceae): seed dispersal. J. Torrey Bot. Soc. 129, 21–28 (2002).Article 

Google Scholar 
36.Song, B.-H. et al. Multilocus patterns of nucleotide diversity, population structure, and linkage disequilibrium in Boechera stricta, a wild relative of Arabidopsis. Genetics 181, 1021–1033 (2009).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
37.Mackay, T., Stone, E. & Ayroles, J. The genetics of quantitative traits: challenges and prospects. Nat. Rev. Genet. 10, 565–577 (2009).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
38.Hedrick, P. W. Genetic polymorphism in heterogeneous environments: a decade later. Annu. Rev. Ecol. Syst. 17, 535–566 (1986).Article 

Google Scholar 
39.Hedrick, P. W. Antagonistic pleiotropy and genetic polymorphism: a perspective. Heredity 82, 126–133 (1999).Article 

Google Scholar 
40.Turelli, M. & Barton, N. H. Polygenic variation maintained by balancing selection: pleiotropy, sex-dependent allelic effects and G × E interactions. Genetics 166, 1053–1079 (2004).PubMed 
PubMed Central 
Article 

Google Scholar 
41.Gillespie, J. H. & Langley, C. H. A general model to account for enzyme variation in natural populations. Genetics 76, 837–848 (1974).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
42.Anderson, J. T., Willis, J. H. & Mitchell-Olds, T. Evolutionary genetics of plant adaptation. Trends Genet. 27, 258–266 (2011).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
43.Anderson, J. T., Lee, C.-R., Rushworth, C. A., Colautti, R. I. & Mitchell-Olds, T. Genetic trade-offs and conditional neutrality contribute to local adaptation. Mol. Ecol. 22, 699–708 (2013).PubMed 
Article 
PubMed Central 

Google Scholar 
44.Oakley, C. G., Ågren, J., Atchison, R. A. & Schemske, D. W. QTL mapping of freezing tolerance: links to fitness and adaptive trade-offs. Mol. Ecol. 23, 4304–4315 (2014).PubMed 
Article 
PubMed Central 

Google Scholar 
45.Price, N. et al. Combining population genomics and fitness QTLs to identify the genetics of local adaptation in Arabidopsis thaliana. Proc. Natl Acad. Sci. USA 115, 5028–5033 (2018).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
46.Kettunen, J. et al. Genome-wide association study identifies multiple loci influencing human serum metabolite levels. Nat. Genet. 44, 269–276 (2012).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
47.Abuelsoud, W., Hirschmann, F. & Papenbrock, J. in Drought Stress in Plants Vol. 1 (eds Hossain, M. A. et al.) 227–248 (Springer, 2016).48.Nguyen, D., Rieu, I., Mariani, C. & van Dam, N. M. How plants handle multiple stresses: hormonal interactions underlying responses to abiotic stress and insect herbivory. Plant Mol. Biol. 91, 727–740 (2016).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
49.Shani, E. M. et al. Plant stress tolerance requires auxin-sensitive Aux/IAA transcriptional repressors. Curr. Biol. 27, 437–444 (2017).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
50.Hopkins, R. J., van Dam, N. M. & van Loon, J. J. A. Role of glucosinolates in insect–plant relationships and multitrophic interactions. Annu. Rev. Entomol. 54, 57–83 (2009).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
51.Burow, M., Müller, R., Gershenzon, J. & Wittstock, U. Altered glucosinolate hydrolysis in genetically engineered Arabidopsis thaliana and its influence on the larval development of Spodoptera littoralis. J. Chem. Ecol. 32, 2333–2349 (2006).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
52.Wagner, M. R. & Mitchell-Olds, T. Plasticity of plant defense and its evolutionary implications in wild populations of Boechera stricta. Evolution 72, 1034–1049 (2018).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
53.Purcell, S. et al. PLINK: a tool set for whole-genome association and population-based linkage analyses. Am. J. Hum. Genet. 81, 559–575 (2007).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
54.Pagès, H., Aboyoun, P., Gentleman, R. & DebRoy, S. Biostrings: Efficient manipulation of biological strings. R package version 2.56.0 (2020).55.Wang et al. Correction to: Ancient polymorphisms contribute to genome-wide variation by long-term balancing selection and divergent sorting in Boechera stricta. Genome Biol. 20, 16 (2019).Article 

Google Scholar 
56.Carley, L. et al. Data to accompany: Ecological factors influence balancing selection on leaf chemical profiles of a wildflower. Dryad Data https://doi.org/10.5061/dryad.7h44j0zsr (2021).57.Atkinson, N. J., Lilley, C. J. & Urwin, P. E. Identification of genes involved in the response of Arabidopsis to simultaneous biotic and abiotic stresses. Plant Physiol. 162, 2028–2041 (2013).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
58.Sharma, A. et al. Comprehensive analysis of plant rapid alkalization factor (RALF) genes. Plant Physiol. Biochem. 106, 82–90 (2016).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
59.Dutilleul, C., Jourdain, A., Bourguignon, J. & Hugouvieux, V. The Arabidopsis putative selenium-binding protein family: expression study and characterization of SBP1 as a potential new player in cadmium detoxification processes. Plant Physiol. 147, 239–251 (2008).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
60.Jiang, S.-C. et al. Crucial roles of the pentatricopeptide repeat protein SOAR1 in Arabidopsis response to drought, salt and cold stresses. Plant Mol. Biol. 88, 369–385 (2015).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
61.Wen, J., Vanek-Krebitz, M., Hoffmann-Sommergruber, K., Scheiner, O. & Breitender, H. The potential of Betv1 homologues, a nuclear multigene family, as phylogenetic markers in flowering plants. Mol. Phylogenet. Evol. 8, 317–333 (1997).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
62.Koo, A. J., Fulda, M., Browse, J. & Ohlrogge, J. B. Identification of a plastid acyl‐acyl carrier protein synthetase in Arabidopsis and its role in the activation and elongation of exogenous fatty acids. Plant J. 44, 620–632 (2005).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
63.Henrissat, B. et al. Conserved catalytic machinery and the prediction of a common fold for several families of glycosyl hydrolases. Proc. Natl Acad. Sci. USA 92, 7090–7094 (1995).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar  More

A uniform response of the pollen grains towards storage conditions was registered in all five shrubs investigated with a conspicuous decline in germination percentage and pollen tube length after storage. Pollen tube growth reacted more sensitively to storage than germination. The most profound reductions in pollen viability traits were observed in samples stored at + 4 °C. The germination percentage of freshly collected pollen of individual shrubs ranged between 67.3 and 88.6%, whereas that in stored pollen was between 18.0 and 39.6%. In relative terms, storage represented a 49.3–73.2% decline in germination (Fig. 1). The same tendency was also observed in pollen tube growth, when freshly collected pollen possessed 248.0–367.3 µm long pollen tubes, and pollen stored at + 4 °C was characterised by 93.9–218.5 µm long pollen tubes. The corresponding decline reached 32.5–68.7%.Figure 1Graphical illustrations of variation in pollen germination percentage (a) and pollen tube length (b) of individual shrubs revealed in fresh pollen and in pollen under storage. Different letters refer to the statistical significance of the differences between tested individuals and storage variants, resulting from Duncan’s pairwise tests.Full size imageContrary to storage at + 4 °C, pollen stored at − 20 °C had an increased germination by 0.3% in shrub no. 1 and 0.6% in shrub no. 5 as compared with fresh pollen. A more conspicuous increase in pollen germinability was registered in individual no. 4, exhibiting 70.0% germination in fresh pollen and 93.6% in pollen stored at − 20 °C. In the remaining two shrubs (no. 2, 3), only a negligible decline in pollen germination was recorded. The deviation from freshly collected pollen varied within 0.5–16.8%. In general, the germination characteristics of pollen stored at − 20 °C were comparable with those of the fresh pollen and varied between 67.6 and 93.6%. As a second viability trait, pollen tube growth deviated more profoundly from that of fresh pollen than germination. On average, the pollen tube length of pollen stored at − 20 °C ranged from 163.0 to 286.6 µm, which represents a 11.4–45.7% decline compared to fresh pollen (Figs. 1, S1). ANOVA and Duncan`s grouping confirmed the highly significant differences between tested shrubs in both pollen germination percentage (P  More

1.Kallmeyer J, Pockalny R, Adhikari RR, Smith DC, D’Hondt S. Global distribution of microbial abundance and biomass in subseafloor sediment. Proc Natl Acad Sci USA. 2012;109:16213–6.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
2.Lloyd KG, May MK, Kevorkian RT, Steen AD. Meta-analysis of quantification methods shows that Archaea and Bacteria have similar abundances in the subseafloor. Appl Environ Microbiol. 2013;79:7790–9.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
3.Hoshino T, Inagaki F. Abundance and distribution of Archaea in the subseafloor sedimentary biosphere. ISME J. 2019;13:227–31.PubMed 
Article 
PubMed Central 

Google Scholar 
4.Lipp JS, Morono Y, Inagaki F, Hinrichs K-U. Significant contribution of Archaea to extant biomass in marine subsurface sediments. Nature. 2008;454:991–4.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
5.Vuillemin A, Wankel SD, Coskun ÖK, Magritsch T, Vargas S, Estes ER, et al. Archaea dominate oxic subseafloor communities over multimillion-year time scales. Sci Adv. 2019;5:eaaw4108.PubMed 
PubMed Central 
Article 

Google Scholar 
6.Zhao R, Hannisdal B, Mogollon JM, Jørgensen SL. Nitrifier abundance and diversity peak at deep redox transition zones. Sci Rep. 2019;9:8633.PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
7.Hiraoka S, Hirai M, Matsui Y, Makabe A, Minegishi H, Tsuda M, et al. Microbial community and geochemical analyses of trans-trench sediments for understanding the roles of hadal environments. ISME J. 2020;14:740–56.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
8.Hoshino T, Doi H, Uramoto GI, Wörmer L, Adhikari RR, Xiao N, et al. Global diversity of microbial communities in marine sediment. Proc Natl Acad Sci. 2020;117:27587–97.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
9.Durbin AM, Teske A. Archaea in organic-lean and organic-rich marine subsurface sediments: an environmental gradient reflected in distinct phylogenetic lineages. Front Microbiol. 2012;3:168.PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
10.Biddle JF, Lipp JS, Lever MA, Lloyd KG, Sørensen KB, Anderson R, et al. Heterotrophic archaea dominate sedimentary subsurface ecosystems off Peru. Proc Natl Acad Sci USA. 2006;103:3846–51.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
11.Lloyd KG, Schreiber L, Petersen DG, Kjeldsen KU, Lever MA, Steen AD, et al. Predominant archaea in marine sediments degrade detrital proteins. Nature. 2013;496:215–8.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
12.Yu T, Wu W, Liang W, Lever MA, Hinrichs K-U, Wang F. Growth of sedimentary Bathyarchaeota on lignin as an energy source. Proc Natl Acad Sci. 2018;115:6022–7.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
13.Zaremba-Niedzwiedzka K, Caceres EF, Saw JH, Bäckström D, Juzokaite L, Vancaester E, et al. Asgard archaea illuminate the origin of eukaryotic cellular complexity. Nature. 2017;541:353–8.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
14.Spang A, Saw JH, Jørgensen SL, Zaremba-Niedzwiedzka K, Martijn J, Lind AE, et al. Complex archaea that bridge the gap between prokaryotes and eukaryotes. Nature. 2015;521:173–9.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
15.Spang A, Caceres EF, Ettema TJG. Genomic exploration of the diversity, ecology, and evolution of the archaeal domain of life. Science. 2017;357:eaaf3883.PubMed 
Article 
CAS 
PubMed Central 

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

Google Scholar 
17.Manoharan L, Kozlowski JA, Murdoch RW, Löffler FE, Sousa FL, Schleper C. Metagenomes from coastal marine sediments give insights into the ecological role and cellular features of Loki-and Thorarchaeota. mBio. 2019;10:e02039–02019.CAS 
PubMed 
PubMed Central 

Google Scholar 
18.Imachi H, Nobu MK, Nakahara N, Morono Y, Ogawara M, Takaki Y, et al. Isolation of an archaeon at the prokaryote–eukaryote interface. Nature. 2020;577:519–25.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
19.Seitz KW, Dombrowski N, Eme L, Spang A, Lombard J, Sieber JR, et al. Asgard archaea capable of anaerobic hydrocarbon cycling. Nat Commun. 2019;10:1822.PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
20.Farag IF, Zhao R, Biddle JF. “Sifarchaeota” a novel Asgard phylum from Costa Rican sediment capable of polysaccharide degradation and anaerobic methylotrophy. Appl Environ Microbiol. 2021;87:e02584–02520.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
21.Spang A, Stairs CW, Dombrowski N, Eme L, Lombard J, Caceres EF, et al. Proposal of the reverse flow model for the origin of the eukaryotic cell based on comparative analyses of Asgard archaeal metabolism. Nat Microbiol. 2019;4:1138–48.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
22.Laso-Pérez R, Wegener G, Knittel K, Widdel F, Harding KJ, Krukenberg V, et al. Thermophilic archaea activate butane via alkyl-coenzyme M formation. Nature. 2016;539:396–401.PubMed 
Article 
CAS 
PubMed Central 

Google Scholar 
23.Chen S-C, Musat N, Lechtenfeld OJ, Paschke H, Schmidt M, Said N, et al. Anaerobic oxidation of ethane by archaea from a marine hydrocarbon seep. Nature. 2019;568:108–11.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
24.Wang Y, Wegener G, Hou J, Wang F, Xiao X. Expanding anaerobic alkane metabolism in the domain of Archaea. Nat Microbiol. 2019;4:595–602.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
25.Laso-Pérez R, Hahn C, van Vliet DM, Tegetmeyer HE, Schubotz F, Smit NT, et al. Anaerobic degradation of non-methane alkanes by “Candidatus Methanoliparia” in hydrocarbon seeps of the Gulf of Mexico. mBio. 2019;10:e01814–01819.PubMed 
PubMed Central 
Article 

Google Scholar 
26.Krukenberg V, Harding K, Richter M, Glöckner FO, Gruber-Vodicka HR, Adam B, et al. Candidatus Desulfofervidus auxilii, a hydrogenotrophic sulfate‐reducing bacterium involved in the thermophilic anaerobic oxidation of methane. Environ Microbiol. 2016;18:3073–91.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
27.Martino A, Rhodes ME, León-Zayas R, Valente IE, Biddle JF, House CH. Microbial diversity in sub-seafloor sediments from the Costa Rica Margin. Geosciences. 2019;9:218.CAS 
Article 

Google Scholar 
28.Farag IF, Biddle JF, Zhao R, Martino AJ, House CH, León-Zayas RI. Metabolic potentials of archaeal lineages resolved from metagenomes of deep Costa Rica sediments. ISME J. 2020;14:1345–58.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
29.Barry PH, de Moor JM, Giovannelli D, Schrenk M, Hummer DR, Lopez T, et al. Forearc carbon sink reduces long-term volatile recycling into the mantle. Nature. 2019;568:487–92.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
30.Expedition 334 Scientists. Site U1379. In Vannucchi, P, Ujiie, K, Stroncik, N, Malinverno, A, and the Expedition 334 Scientists, Proc IODP, 334: Tokyo (Integrated Ocean Drilling Program Management International, Inc) (2012).31.Formolo M, Nuzzo M, Torres M, Solomon E. Expedition I Gas geochemical results from IODP Expedition 334: Influence of subsurface structure and fluid flow on gas composition. In: Proceedings of AGU Fall Meeting Abstracts) 2011.32.Boyd JA, Woodcroft BJ, Tyson GW. GraftM: a tool for scalable, phylogenetically informed classification of genes within metagenomes. Nucleic Acids Res. 2018;46:e59.PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
33.Singleton CM, McCalley CK, Woodcroft BJ, Boyd JA, Evans PN, Hodgkins SB, et al. Methanotrophy across a natural permafrost thaw environment. ISME J. 2018;12:2544–58.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
34.Borrel G, Adam PS, McKay LJ, Chen LX, Sierra-García IN, Sieber C, et al. Wide diversity of methane and short-chain alkane metabolisms in uncultured archaea. Nat Microbiol. 2019;4:603–13.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
35.Hua Z-S, Wang YL, Evans PN, Qu YN, Goh KM, Rao YZ, et al. Insights into the ecological roles and evolution of methyl-coenzyme M reductase-containing hot spring Archaea. Nat Commun. 2019;10:1–11.Article 
CAS 

Google Scholar 
36.Cai M, et al. Diverse Asgard archaea including the novel phylum Gerdarchaeota participate in organic matter degradation. Science China Life Sciences, (2020).37.Konstantinidis KT, Rosselló-Móra R, Amann R. Uncultivated microbes in need of their own taxonomy. ISME J. 2017;11:2399–406.PubMed 
PubMed Central 
Article 

Google Scholar 
38.Hahn CJ, Laso-Pérez R, Vulcano F, Vaziourakis KM, Stokke R, Steen IH, et al. “Candidatus Ethanoperedens,” a thermophilic genus of Archaea mediating the anaerobic oxidation of ethane. mBio. 2020;11:e00600–00620.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
39.Rastogi S, Liberles DA. Subfunctionalization of duplicated genes as a transition state to neofunctionalization. BMC Evolut Biol. 2005;5:28.Article 
CAS 

Google Scholar 
40.Hug LA, Baker BJ, Anantharaman K, Brown CT, Probst AJ, Castelle CJ, et al. A new view of the tree of life. Nat Microbiol. 2016;1:16048.CAS 
PubMed 
Article 
PubMed Central 

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

Google Scholar 
42.Skennerton CT, Chourey K, Iyer R, Hettich RL, Tyson GW, Orphan VJ. Methane-fueled syntrophy through extracellular electron transfer: uncovering the genomic traits conserved within diverse bacterial partners of anaerobic methanotrophic archaea. mBio. 2017;8:e00530–00517.PubMed 
PubMed Central 

Google Scholar 
43.Beulig F, Røy H, McGlynn SE, Jørgensen BB. Cryptic CH4 cycling in the sulfate–methane transition of marine sediments apparently mediated by ANME-1 archaea. ISME J. 2019;13:250–62.CAS 
PubMed 
Article 
PubMed Central 

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

Google Scholar 
45.Dombrowski N, Teske AP, Baker BJ. Expansive microbial metabolic versatility and biodiversity in dynamic Guaymas Basin hydrothermal sediments. Nat Commun. 2018;9:4999.PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
46.Dong X, Greening C, Rattray JE, Chakraborty A, Chuvochina M, Mayumi D, et al. Metabolic potential of uncultured bacteria and archaea associated with petroleum seepage in deep-sea sediments. Nat Commun. 2019;10:1816.PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
47.Brown CT, Olm MR, Thomas BC, Banfield JF. Measurement of bacterial replication rates in microbial communities. Nat Biotechnol. 2016;34:1256–63.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
48.Greening C, Biswas A, Carere CR, Jackson CJ, Taylor MC, Stott MB, et al. Genomic and metagenomic surveys of hydrogenase distribution indicate H2 is a widely utilised energy source for microbial growth and survival. ISME J. 2016;10:761–77.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
49.Shimoyama T, Kato S, Ishii SI, Watanabe K. Flagellum mediates symbiosis. Science. 2009;323:1574–1574.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
50.Valentine DL, Reeburgh WS. New perspectives on anaerobic methane oxidation: minireview. Environ Microbiol. 2000;2:477–84.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
51.Vannucchi P, Ujiie K, Stroncik N, the IESP. IODP Expedition 334: An investigation of the sedimentary record, fluid flow and state of stress on top of the seismogenic zone of an erosive subduction margin. Sci Dril. 2013;15:23–30.Article 

Google Scholar 
52.Torres ME, Muratli JM, Solomon EA Data report: minor element concentrations in pore fluids from the CRISP-A transect drilled during Expedition 334. In: Proceeding sof IODP | Volume) 2014.53.Riedinger N, Torres ME, Screaton E, Solomon EA, Kutterolf S, Schindlbeck‐Belo J, et al. Interplay of subduction tectonics, sedimentation, and carbon cycling. Geochem, Geophys, Geosyst. 2019;20:4939–55.CAS 
Article 

Google Scholar 
54.Andrews S. FastQC: a quality control tool for high throughput sequence data. https://www.bioinformatics.babraham.ac.uk/projects/fastqc/ 2010.55.Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30:2114–20.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
56.Gruber-Vodicka HR, Seah BKB, Pruesse E. phyloFlash: rapid small-subunit rRNA profiling and targeted assembly from metagenomes. mSystems. 2020;5:e00920–00920.CAS 
PubMed 
PubMed Central 
Article 

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

Google Scholar 
58.Wu YW, Simmons BA, Singer SW. MaxBin 2.0: an automated binning algorithm to recover genomes from multiple metagenomic datasets. Bioinformatics. 2016;32:605–7.CAS 
PubMed 
Article 
PubMed Central 

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

Google Scholar 
60.Seah BK, Gruber-Vodicka HR. gbtools: interactive visualization of metagenome bins in R. Front. Microbiol. 2015;6:1451.PubMed 
PubMed Central 
Article 

Google Scholar 
61.Bushnell B. BBMap: a fast, accurate, splice-aware aligner. Ernest Orlando Lawrence Berkeley National Laboratory, Berkeley, CA (US) (2014).62.Bankevich A, Nurk S, Antipov D, Gurevich AA, Dvorkin M, Kulikov AS, et al. SPAdes: A new genome assembly algorithm and its applications to single-cell sequencing. J Comput Biol. 2012;19:455–77.CAS 
PubMed 
PubMed Central 
Article 

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

Google Scholar 
64.Huerta-Cepas J, Szklarczyk D, Forslund K, Cook H, Heller D, Walter MC, et al. eggNOG 4.5: a hierarchical orthology framework with improved functional annotations for eukaryotic, prokaryotic and viral sequences. Nucleic Acids Res. 2016;44:D286–D293.CAS 
Article 

Google Scholar 
65.Kanehisa M, Sato Y, Morishima K. BlastKOALA and GhostKOALA: KEGG tools for functional characterization of genome and metagenome sequences. J Mol Biol. 2016;428:726–31.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
66.Kanehisa M, Goto S, Sato Y, Furumichi M, Tanabe M. KEGG for integration and interpretation of large-scale molecular data sets. Nucleic Acids Research. 2011;40:D109–D114.PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
67.Garcia PS, Jauffrit F, Grangeasse C. Brochier-Armanet C. GeneSpy, a user-friendly and flexible genomic context visualizer. Bioinformatics. 2018;35:329–31.Article 
CAS 

Google Scholar 
68.Badalamenti JP, Summers ZM, Chan CH, Gralnick JA, Bond DR. Isolation and genomic characterization of ‘Desulfuromonas soudanensis WTL’, a metal-and electrode-respiring bacterium from anoxic deep subsurface brine. Front Microbiol. 2016;7:913.PubMed 
PubMed Central 
Article 

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

Google Scholar 
70.Hernsdorf AW, Amano Y, Miyakawa K, Ise K, Suzuki Y, Anantharaman K, et al. Potential for microbial H2 and metal transformations associated with novel bacteria and archaea in deep terrestrial subsurface sediments. ISME J. 2017;11:1915–29.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
71.Jain C, Rodriguez-R LM, Phillippy AM, Konstantinidis KT, Aluru S. High throughput ANI analysis of 90K prokaryotic genomes reveals clear species boundaries. Nat Commun. 2018;9:5114.PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
72.Buchfink B, Xie C, Huson DH. Fast and sensitive protein alignment using DIAMOND. Nat Methods. 2015;12:59–60.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
73.Sorek R, Zhu YW, Creevey CJ, Francino MP, Bork P, Rubin EM. Genome-wide experimental determination of barriers to horizontal gene transfer. Science. 2007;318:1449–52.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
74.Campbell JH, O’Donoghue P, Campbell AG, Schwientek P, Sczyrba A, Woyke T, et al. UGA is an additional glycine codon in uncultured SR1 bacteria from the human microbiota. Proc Natl Acad Sci USA. 2013;110:5540–5.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
75.Eren AM, Esen ÖC, Quince C, Vineis JH, Morrison HG, Sogin ML, et al. Anvi’o: an advanced analysis and visualization platformfor ‘omics data. PeerJ. 2015;3:e1319.PubMed 
PubMed Central 
Article 

Google Scholar 
76.Edgar RC. MUSCLE: multiple sequence alignment with high accuracy and high throughput. Nucleic Acids Res. 2004;32:1792–7.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
77.Capella-Gutierrez S, Silla-Martinez JM, Gabaldon T. trimAl: a tool for automated alignment trimming in large-scale phylogenetic analyses. Bioinformatics. 2009;25:1972–3.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
78.Nguyen LT, Schmidt HA, von Haeseler A, Minh BQ. IQ-TREE: a fast and effective stochastic algorithm for estimating maximum-likelihood phylogenies. Mol Biol Evolut. 2015;32:268–74.CAS 
Article 

Google Scholar 
79.Kalyaanamoorthy S, Minh BQ, Wong TKF, von Haeseler A, Jermiin LS. ModelFinder: fast model selection for accurate phylogenetic estimates. Nat Methods. 2017;14:587–9.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
80.Hoang DT, Chernomor O, von Haeseler A, Minh BQ, Vinh LS. UFBoot2: improving the ultrafast bootstrap approximation. Mol Biol Evolut. 2018;35:518–22.CAS 
Article 

Google Scholar 
81.Altschul SF, Madden TL, Schäffer AA, Zhang J, Zhang Z, Miller W, et al. Gapped BLAST and PSI-BLAST: a new generation of protein database search programs. Nucleic Acids Res. 1997;25:3389–402.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
82.Katoh K, Standley DM. MAFFT multiple sequence alignment software version 7: improvements in performance and usability. Mol Biol Evolut. 2013;30:772–80.CAS 
Article 

Google Scholar 
83.Okonechnikov K, Golosova O, Fursov M, Team U. Unipro UGENE: a unified bioinformatics toolkit. Bioinformatics. 2012;28:1166–7.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
84.Matheus Carnevali PB, Schulz F, Castelle CJ, Kantor RS, Shih PM, Sharon I, et al. Hydrogen-based metabolism as an ancestral trait in lineages sibling to the Cyanobacteria. Nat Commun. 2019;10:463.CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
85.Kessler AJ, Chen YJ, Waite DW, Hutchinson T, Koh S, Popa ME, et al. Bacterial fermentation and respiration processes are uncoupled in anoxic permeable sediments. Nat Microbiol. 2019;4:1014–23.CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
86.R Development Core Team. R: a language and environment for statistical computing.). R foundation for statistical computing, Vienna, Austria (2011). More

Data for non-human primatesWe obtained 30 datasets for six genera of non-human primates: sifaka (Propithecus spp), gracile capuchin monkey (Cebus spp), guenon (Cercopithecus spp), baboon (Papio spp), gorilla (Gorilla spp), and chimpanzee (Pan troglodytes) (Supplementary Data 1). Of these, 17 datasets correspond to long-term projects in the wild, while 13 were contributed by the non-profit Species360 from ZIMS18, which is the most extensive database of life history information for animals under human care.Basic demographic functionsLet X be a random variable for ages at death, with observations x ≥ 0, and let μ (x|θ) be a continuous, non-negative parametric hazards rate or mortality function defined as$$mu left(x,|,{boldsymbol{theta }}right)=mathop{{rm{lim}}}limits_{Delta xto 0}frac{{{Pr }}(x < Xle x+Delta x|X > x)}{Delta x},$$
(2)
given that the limit exists, where ({boldsymbol{theta }}in {{mathbb{R}}}^{p}) is a p-dimensional vector of mortality parameters. The cumulative hazards rate is$$Uleft(x|{boldsymbol{theta }}right)=int_{0}^{x}mu (t|{boldsymbol{theta }}){dt},$$
(3)
which results in the survival function$$S(x|{boldsymbol{theta }})={{exp }}[-U(x|{boldsymbol{theta }})].$$
(4)
The Cumulative distribution function (CDF) of ages at death is F (x | θ) = 1 – S (x | θ), and the probability density function (PDF) of ages at death is f (x | θ) = μ (x | θ) S (x | θ), for x ≥ 0. The remaining life expectancy after age x is calculated as$$eleft(x{rm{|}}{boldsymbol{theta }}right) = frac{{int }_{x}^{{{infty }}}{tf}(t|{boldsymbol{theta }}){dt}}{Fleft({{infty }}right)-Fleft(xright)}\ =frac{{int }_{x}^{{{infty }}}S(t|{boldsymbol{theta }}){dt}}{Sleft(xright)},$$
(5)
which yields a life expectancy at birth given by$$eleft(0{rm{|}}{boldsymbol{theta }}right)={int }_{0}^{{{infty }}}S(x|{boldsymbol{theta }}){dx}.$$
(6)
The lifespan inequality at birth, as proposed by Demetrius16,36 and later by Keyfitz17, is given by$$H(0|{boldsymbol{theta }}) =-frac{{int }_{0}^{{{infty }}}S(x|{boldsymbol{theta }}){{log }}[S(x|{boldsymbol{theta }})]{dx}}{e(0|{boldsymbol{theta }})}\ = frac{{int }_{0}^{{{infty }}}S(x|{boldsymbol{theta }})U(x|{boldsymbol{theta }}){dx}}{e(0|{boldsymbol{theta }})}.$$
(7)
Following Colchero et al.13, we define the lifespan equality as$$varepsilon (x|{boldsymbol{theta }})=-{{log }}[H(x|{boldsymbol{theta }})].$$
(8)
For simplicity, henceforth we note the life expectancy, lifespan inequality and lifespan equality at birth as e(0 | θ) = e, H (0 | θ) = H, and ε (0 | θ) = ε, respectively.Survival analysisTo estimate age-specific survival for all the wild populations of non-human primates, we modified the Bayesian model developed by Colchero et al.13 and Barthold et al.37. This model is particularly appropriate for primate studies that follow individuals continuously within a study area and when individuals of one or both sexes can permanently leave the study area (out-migration), while other individuals can join the study population from other areas (in-migration). Thus, it allowed us to make inferences on age-specific survival (or mortality) and on the age at out-migration.Here we use the five parameter Siler mortality function25, as in Eq. (1) where θ = [a0, a0, c, b0, b1] is a vector of parameters to be estimated, and where a0, b0 ({mathbb{in }}{mathbb{R}}) and a1, c, b1 ≥ 0. For all species we studied, individuals of one or both sexes often leave their natal groups to join other neighbouring groups in a process commonly identified as natal dispersal. For some species, individuals who have undergone natal dispersal can then disperse additional times, described as secondary dispersal. Although dispersal within monitored groups (i.e. those belonging to the study area) does not affect the estimation of mortality, the fate of individuals that permanently leave the study area to join unmonitored groups can be mistaken for possible death. We identify this process as “out-migration”, which we classify as natal or immigrant out-migration, the first for natal and the second for secondary dispersals to unmonitored groups. This distinction is particularly relevant because not all out-migrations are identified as such, and therefore the fate of some individuals is unknown after their last detection. For these individuals we define a latent out-migration state at the time they were last detected, given by the random variable indicator O, with observations oij ∈ {0,1}, where oij = 1 if individual i out-migrated and oij = 0 otherwise, and where j = 1 denotes natal out-migration and j = 2 for immigrant out-migration. For known out-migrations, we automatically assign oij = 1. The model therefore estimates the Bernoulli probability of out-migration, πj, such that Oij ~ Bern(πj). Those individuals assigned as exhibiting out-migration, as well as known emigrants and immigrants, contribute to the estimation of the distribution of ages at out-migration. Here, we define a gamma-distributed random variable V for ages at out-migration, with realisations v ≥ 0, where Vj | Oj = 1 ~ Gam(γj1, γj2) and where γj1, γj2  > 0 are parameters to be estimated with j defined as above. The probability density function for the gamma distribution is gV(v | γj1, γj2) for v ≥ 0, with v = xl – αj, where xl is the age at last detection and αj is the minimum age at natal or immigrant out-migration.In addition, since not all individuals have known birth dates, the model samples the unknown births bi as xil = til – bi, where til is the time of last detection for individual i. The likelihood is then defined as$$p({x}_{{il}},{x}_{{if}},|,{boldsymbol{theta}},{boldsymbol{gamma}}_{1},{boldsymbol{gamma}}_{2},{pi }_{j},{o}_{ij})=left{begin{array}{cc}frac{fleft({x}_{il}right)}{Sleft({x}_{if}right)}({1-pi }_{j})hfill& {text{if}}; o_{{ij}}=0\ frac{Sleft({x}_{{il}}right)}{Sleft({x}_{{if}}right)}{pi }_{j}{g}_{V}({x}_{{il}}-{alpha }_{j})& {text{if}}; o_{{ij}}=1end{array}right.,$$
(9)
where xif is the age at first detection, given by xif = tif – bi, with tif as the corresponding time of first detection. The parameter vectors γ1 and γ2 are for natal and immigrant out-migration, respectively. In other words, individuals with oij = 0 are assumed to have died shortly after the last detection, while those with oij = 1 are censored and contribute to the estimation of the distribution of ages at out-migration. The full Bayesian posterior is then given by$$pleft({boldsymbol{theta }}{boldsymbol{,}}{{boldsymbol{gamma }}}_{1},{{boldsymbol{gamma }}}_{2},{boldsymbol{pi }},{{bf{b}}}_{u},{{bf{o}}}_{u},|,{{bf{b}}}_{k},{{bf{o}}}_{k},{{bf{t}}}_{f},{{bf{t}}}_{l}right) propto ; pleft({{bf{x}}}_{l},{{bf{x}}}_{f},|,{boldsymbol{theta }},{{boldsymbol{gamma }}}_{1},{{boldsymbol{gamma }}}_{2},{boldsymbol{pi }},{bf{d}}right)\ , times pleft({boldsymbol{theta }}right)pleft({{boldsymbol{gamma }}}_{1}right)pleft({{boldsymbol{gamma }}}_{2}right)pleft({boldsymbol{pi }}right),$$
(10)
where the first term on the right-hand-side of Eq. (10) is the likelihood in Eq. (9), and the following terms are the priors for the unknown parameters. The vector π = [π1, π2] is the vector of probabilities of out-migration while the subscripts u and k refer to unknown and known, respectively.Following Colchero et al.13, we used published data, expert information and an agent-based model to estimate the mortality and out-migration prior parameters for each population. We assumed a normal (or truncated normal distribution depending on the parameter’s support) for all the parameters. We used vague priors for the mortality and natal out-migration parameters (sd = 10), and informative priors for the immigrant out-migration parameters (sd = 0.5). We ran six MCMC parallel chains for 25 000 iterations each with a burn-in of 5000 iterations for each population, and assessed convergence using potential scale reduction factor38.For the zoo data we used a simplified version of the model described above, which omitted all parts that related to out-migration. In order to produce Supplementary Figs. 1 and 2, we used the same method as for the zoo data on the human life tables. To achieve this, we created an individual level dataset from the lx column of each population, and then fitted the Siler model to this simulated data. It is important to note that the Siler model provides a close fit to the nonhuman primate data and to high-mortality human populations, although it does not provide the best fit to low-mortality human populations, in part due to the late life mortality plateau common among human populations39 (Supplementary Fig. 6). It is therefore possible that the values of the mortality parameter b1 we report in Supplementary Data 2 for the human populations are under-estimated. Nonetheless, and for the purposes of our analyses, the Siler fits to the human populations we considered here are reasonable (Supplementary Fig. 6) and we can therefore confidently state that the limitations of the Siler model do not affect the generality of our results.Estimation of life expectancy and lifespan equalityBased on the results of the Bayesian inference models, we calculated life expectancy at birth as$$e= int_{0}^{{infty}}Sleft(t| {hat{boldsymbol{theta }}}right){dt},$$
(11)
where S (x) is the cumulative survival function as defined in Eq. (4) and where (hat{{boldsymbol{theta }}}) is the vector of mortality parameters calculated as the mean of the conditional posterior densities from the survival analysis described above. We calculated the lifespan inequality17,36, H, as$$H=-frac{1}{e}int_{0}^{{{infty }}}Sleft(x{rm{|}}hat{{boldsymbol{theta }}}right){{log }}left[Sleft(x|hat{{boldsymbol{theta }}}right)right]{dx},$$
(12)
from which we calculated lifespan equality, ε, as in Eq. (8). We calculated both measures for each of the study populations, and performed weighted least squares regressions for each genus, with weights given by the reciprocal of the standard error of the estimated life expectancies.Sensitivities of life expectancy and lifespan equality to mortality parametersAs we mentioned above, for simplicity of notation, we will express all demographic functions by their variable notation (e.g. e = e (0 | θ), S = S (x | θ), etc.), while we will alternatively note first partial derivatives, for instance the derivative of e with respect to a given mortality parameter θ ∈ θ, as eθ or ∂e / ∂θ.Proposition: If ({S:}{{mathbb{R}}}_{ge 0}to left[{mathrm{0,1}}right]) is a continuous non-increasing parametric survival function with parameter vector ({boldsymbol{theta }}{boldsymbol{in }}{{mathbb{R}}}^{{boldsymbol{p}}}), with continuous differentiable cumulative hazards function ({U:}{{mathbb{R}}}_{ge 0}to {{mathbb{R}}}_{ge 0}), and with life expetancy at birth, lifespan inequality and lifespan equality as in Eqs. (4)-(6), respectively, then the sensitivity of life expectancy, e, to a given parameter θ ∈ θ is$${e}_{theta }=frac{partial e}{partial theta }=int_{0}^{{{infty }}}{S}_{theta }{dx},$$
(13)
while the sensitivity of lifespan equality to θ is$${varepsilon }_{theta }=frac{partial varepsilon }{partial theta }=frac{{e}_{theta }left(1+{H}^{-1}right)-{H}^{-1}{int }_{0}^{{{infty }}}{S}_{theta }{Udx}}{e},$$
(14)
where$${S}_{theta }=frac{partial }{partial theta }S(x|{boldsymbol{theta }})$$
(15)
is the sensitivity of the survival function at age x to changes in parameter θ.Proof. The sensitivity of lifespan equality to changes in θ is derived from$${e}_{theta }=frac{partial }{partial theta }int_{0}^{{{infty }}}{Sdx},$$
(16)
which, by Leibnitz’s rule, Eq. (16) becomes$${e}_{theta }=int_{0}^{{{infty }}}frac{partial S}{partial theta }{dx}=int_{0}^{{{infty }}}{S}_{theta }{dx}.$$
(17)
The sensitivity of lifespan equality to changes in θ can be calculated as$${varepsilon }_{theta } =frac{partial }{partial theta }left(-{{log }}, Hright)\ =-frac{partial }{partial theta }{{log }}, H\ =-frac{1}{H}frac{partial H}{partial theta }\ =-frac{1}{H}frac{partial }{partial theta }left(frac{{int }_{0}^{{{infty }}}{SUdx}}{e}right).$$
(18)
By the quotient and Leibnitz’s rules, Eq. (18) can be modified as$${varepsilon }_{theta } =-frac{1}{H{e}^{2}}left[frac{partial }{partial theta }left(int _{0}^{{{infty }}}{SUdx}right)e-left(int _{0}^{{{infty }}}{SUdx}right)frac{partial e}{partial theta }right]\ =-frac{1}{{He}}int _{0}^{{{infty }}}frac{partial }{partial theta }left({SU}right){dx}+frac{1}{{He}}frac{int _{0}^{{{infty }}}{SUdx}}{e}frac{partial e}{partial theta }.$$
(19)
The first term in Eq. (19) can be further decomposed by the product rule, while the second term can be modified following the equality for H in Eq. (7), which yields$${varepsilon }_{theta } =-frac{1}{{He}}int _{0}^{{{infty }}}left(frac{partial S}{partial theta }U+Sfrac{partial U}{partial theta }right){dx}+frac{1}{e}{e}_{theta }\ =-frac{1}{{He}}left(int _{0}^{{{infty }}}{S}_{theta }{Udx}+int _{0}^{{{infty }}}Sfrac{partial U}{partial theta }dxright)+frac{1}{e}{e}_{theta }.$$
(20)
By the chain rule, we have that (frac{partial U}{partial theta }=-frac{partial }{partial theta }{{log }},S=-frac{1}{S}frac{partial S}{partial theta }), which modifies Eq. (20) as$${varepsilon }_{theta } = , -frac{1}{{He}}left(int _{0}^{infty }{S}_{theta }{Udx}-int _{0}^{infty }frac{partial S}{partial theta }{dx}right)+frac{1}{e}{e}_{theta }\ = , -frac{1}{{He}}left(int _{0}^{infty }{S}_{theta }{Udx}-{e}_{theta }right)+frac{1}{e}{e}_{theta }\ = , -frac{int _{0}^{infty }{S}_{theta }{Udx}}{{He}}+frac{{e}_{theta }}{e}left(1+frac{1}{H}right)\ =, frac{{e}_{theta }left(1+{H}^{-1}right)-{H}^{-1}int _{0}^{infty }{S}_{theta }{Udx}}{e},$$
(21)
hence completing the proof. ∎Changes in parameters along the genus linesFrom the results in Eqs. (13) and (14), we calculated the vectors of change (gradient vectors) at any point (leftlangle {e}_{j},{varepsilon }_{j}rightrangle) of the life expectancy-lifespan equality landscape, as a function of each of the Siler mortality parameters (See Fig. 2A, B).To quantify the amount of change of each parameter along the genus lines, we derived the sensitivities of a given mortality parameter θ to changes in life expectancy and lifespan equality, namely (frac{partial theta }{partial e}=frac{1}{{e}_{theta }}) for ({e}_{theta }, ne, 0,) and (frac{partial theta }{partial varepsilon }=frac{1}{{varepsilon }_{theta }}) for ({varepsilon }_{theta }, ne, 0). With these sensitivities we calculated the gradient vector$$nabla theta =leftlangle frac{partial theta }{partial e},frac{partial theta }{partial varepsilon }rightrangle$$
(22)
for any parameter at any point along the genus lines. Here we find a linear relationship between life expectancy and lifespan equality, given by$$mleft({e}_{{ik}}right)={hat{varepsilon }}_{{ik}}={beta }_{0k}+{beta }_{1k}{e}_{{ik}},$$
(23)
for i = 1, …, nk, where nk is the number of populations for genus k, and ({hat{varepsilon }}_{{ik}})is the fitted value of lifespan equality for population i in genus k, and β0k and β1k are linear regresssion parameters for genus k. To estimate the amount of change in parameter θ along the line for genus k, we can solve the path integral$${Theta }_{k}=int _{{C}_{k}}nabla theta d{bf{r}},$$
(24)
where path Ck is determined by the linear model for genus k and (d{bf{r}}=leftlangle {de},dhat{varepsilon }rightrangle =leftlangle {de},{d; m}left(eright)rightrangle) is the rate of change in the velocity vector ({bf{r}}=leftlangle e,hat{varepsilon }rightrangle =leftlangle e,mleft(eright)rightrangle).In order to compare results between the different mortality parameters in vector θ, we use the transformation g(θ) = log θ, which yields the following partial derivatives$$frac{partial }{partial e}gleft(theta right)=frac{1}{theta }frac{partial theta }{partial e}$$
(25)
and$$frac{partial }{partial varepsilon }gleft(theta right)=frac{1}{theta }frac{partial theta }{partial varepsilon }.$$
(26)
Thus the gradient vector becomes$$nabla theta =leftlangle frac{partial }{partial e}gleft(theta right),frac{partial }{partial varepsilon }gleft(theta right)rightrangle$$
(27)
while the path integral in Eq. (24) is modified accordingly. In short, the path integral ({Theta }_{j}) provides a measure of the relative change in parameter θ along the genus line (Fig. 3). To allow comparisons between all genera, we scaled the values of each path integral by the length of each line.Applications to the Siler mortality modelThe Cumulative hazards for the Siler mortality model in Eq. (7) is given by$$Uleft(xright)=frac{{e}^{{a}_{0}}}{{a}_{1}}left(1-{e}^{{-a}_{1}x}right)+{cx}+frac{{e}^{{b}_{0}}}{{b}_{1}}left({e}^{{b}_{1}x}-1right),$$
(28)
The sensitivities in Eqs. (13) and (14) require calculating Sθ for all θ ∈ θ. Treating S (x) as the function composition W (V), where W = exp(x) and V = – U, then Sθ is$${S}_{theta }=frac{{dW}}{{dV}}{V}_{theta }=-S{U}_{theta },$$
(29)
where Uθ is the first derivative of U(x | θ) with respect to θ. For each of the Siler mortality parameters, we then have$${S}_{{a}_{0}}=S(x|{boldsymbol{theta }})frac{{e}^{{a}_{0}}}{{a}_{1}}left({e}^{-{a}_{1}x}-1right)$$
(30)
$${S}_{{a}_{1}}=S(x|{boldsymbol{theta }})frac{{e}^{{a}_{0}}}{{a}_{1}}left[frac{1}{{a}_{1}}-{e}^{-{a}_{1}x}left(x+frac{1}{{a}_{1}}right)right]$$
(31)
$${S}_{c}=-S(x|{boldsymbol{theta }})x$$
(32)
$${S}_{{b}_{0}}=S(x|{boldsymbol{theta }})frac{{e}^{{b}_{0}}}{{b}_{1}}left(1-{e}^{{b}_{1}x}right)$$
(33)
$${S}_{{b}_{1}}=S(x|{boldsymbol{theta }})left[{e}^{{b}_{1}x}left(frac{1}{{b}_{1}}-xright)-frac{1}{{b}_{1}}right].$$
(34)
All analyses were performed in the free open source programme R40. The R functions we created for this project can be found in41.Reporting summaryFurther information on research design is available in the Nature Research Reporting Summary linked to this article. More

Spatial and temporal variations of the SRHII at daytime and nighttimeSignificant seasonal differences are observed in the SRHII in the YRDUA (Figure A1 and A2, Appendix A). In the daytime, RHI was concentrated in the Nanjing, “Su-Xi-Chang”, Ningbo, Shanghai, and Hangzhou metropolitan areas. Due to the high built-up areas and PD, the distribution of surface RHI is denser and stronger than that in the north and southwest of the YRDUA. The built-up area can absorb heat and store heat energy, which makes the surface temperatures rise rapidly. In spring and autumn, the spatial distribution of the RHI in spring or autumn was similar to that in summer except the spatial extent was tapered. However, the RHI gradually shrinks and transfers to the southern area of the YRDUA in winter, such as Linhai and Ningbo City, which is due to the relatively high solar radiation of the geographic location of the southern cities. The distance of the RHI is gradually shortened between cities and even into one piece from 2003 to 2017 due to long-term urban expansion and rapid growth of construction land (Figure A1, Appendix A). In the nighttime, the spatial pattern of the RHI is very different from that of the daytime. RHI mainly concentrates on Taihu Lake, Dianshan Lake, Ge Lake in the center part, Hongzhe Lake in the northwest, and Qiandao Lake in the southwest. Because water has a high specific heat capacity, it has the function of preserving heat at nighttime. Some cities like Shanghai, Hangzhou, and Nanjing have the strongest heat island in winter and the weakest heat island in summer. Urban areas usually have dense buildings, PD, and energy emissions, so there are more energy emissions at night. High surface albedo in urban areas at night leads to lower heat storage4,40 and ultimately resulting in smaller UHI at nighttime (Figure A2, Appendix A).From spring to summer and then summer to winter, RHI increases first and then decreases, and it reaches a peak in summer. For example, the proportion of the RHI was 12.65%, 31.03%, 21.12%, and 5.49% in spring, summer, autumn, and winter in 2017, respectively (Fig. 2d). An upward trend in the area of the RHI is observed from 2003 to 2017 in summer. In detail, the proportion of the heat island zone is 21.74%, 22.17%, and 31.03% in the summer of 2003, 2010, and 2017, respectively (Fig. 2d). It is because the urban areas of YRDUA have increased from 3571.01 km2 to 8760.26 km2 in 2003 and 2017, respectively (Figure B1, Appendix B). Moreover, the area of the medium heat island and strong heat island increased by 41.08% and 66.40% from 2003 to 2017 (Fig. 2b,c). A gradual decreasing trend is observed for the four grades of the SRHII (2–4 °C, 4–6 °C,  > 6 °C,  > 2 °C) in winter from 2003 to 2017 (Fig. 2a–d). The area of the RHI in winter was 18,481 km2, 8640 km2, and 6280 km2 in 2003, 2010, and 2017, respectively (Fig. 2d). Vegetation coverage is low in winter and bare soil is formed after harvest. It leads to the RHI decrease in winter. The above results indicated that the SRHII became increasingly hot in summer and increasingly cold in winter and that the trend became more obvious as the SRHII increased in the ranges of 2–4 °C, 4–6 °C,  > 6 °C. However, the seasonal variation of the RHI in the nighttime is opposite to that in the daytime. From spring to summer and then to winter, the area of the RHI decreases first and then increases, and it falls in the lowest value in summer (Fig. 2e–g). For example, the area of RHI is 19,209 km2, 5659 km2, 34,621 km2, and 38,596 km2 in spring, summer, autumn, and winter in 2017, respectively (Fig. 2h). The annual average of RHI regular increases, with values of 17,510 km2, 20,042 km2, and 20,097 km2 in 2003, 2010, and 2017, respectively (Fig. 2h).Figure 2Seasonal and inter-annual variations of the SRHII during the daytime (a–d) and nighttime (e–h) of the YRDUA.Full size imageRelationship between the SRHII and influencing factorsResults showed surface biophysical factors have a higher correlation with SRHII than socio-economic factors and climate factors in the day and night. NDBI and EVI have a stronger effect on SRHII than other biophysical factors in the day. NDBI showed a significant positive correlation with SRHII, while EVI showed a negative correlation with SRHII. In detail, NDBI (r = 0.567, p  autumn  > winter. The dominant influencing factor was the MNDWI in spring, autumn and winter, while EVI had the largest contribution in summer at night. More