Ecology

Di Giacomo, G., Hadrich, J., Hutchison, W. D., Peterson, H. & Rogers, M. Economic impact of spotted wing drosophila (Diptera: Drosophilidae) yield loss on minnesota raspberry farms: A grower survey. J. Integr. Pest Manag. https://doi.org/10.1093/jipm/pmz006 (2019).Article 

Google Scholar 
Farnsworth, D. et al. Economic analysis of revenue losses and control costs associated with the spotted wing drosophila, Drosophila suzukii (Matsumura), in the California raspberry industry. Pest. Manag. Sci. 73, 1083–1090. https://doi.org/10.1002/ps.4497 (2017).CAS 
Article 
PubMed 

Google Scholar 
Beers, E. H., Van Steenwyk, R. A., Shearer, P. W., Coates, W. W. & Grant, J. A. Developing Drosophila suzukii management programs for sweet cherry in the western United States. Pest Manag. Sci. 67, 1386–1395. https://doi.org/10.1002/ps.2279 (2011).CAS 
Article 
PubMed 

Google Scholar 
Tait, G. et al. Drosophila suzukii (Diptera: Drosophilidae): A decade of research towards a sustainable integrated pest management program. J. Econ. Entomol. 114, 1950–1974. https://doi.org/10.1093/jee/toab158 (2021).Article 
PubMed 

Google Scholar 
Daane, K. M. et al. First exploration of parasitoids of Drosophila suzukii in South Korea as potential classical biological agents. J. Pest Sci. 89, 823–835. https://doi.org/10.1007/s10340-016-0740-0 (2016).ADS 
Article 

Google Scholar 
Abram, P. K. et al. New records of Leptopilina, Ganaspis, and Asobara species associated with Drosophila suzukii in North America, including detections of L. japonica and G. brasiliensis. J. Hymenoptera Res. 78, 1–17. https://doi.org/10.3897/jhr.78.55026 (2020).Article 

Google Scholar 
Chabert, S., Allemand, R., Poyet, M., Eslin, P. & Gibert, P. Ability of European parasitoids (Hymenoptera) to control a new invasive Asiatic pest, Drosophila suzukii. Biol. Control 63, 40–47. https://doi.org/10.1016/j.biocontrol.2012.05.005 (2012).Article 

Google Scholar 
Gonzalez-Cabrera, J., Moreno-Carrillo, G., Sanchez-Gonzalez, J. A., Mendoza-Ceballos, M. Y. & Arredondo-Bernal, H. C. Single and Combined Release of Trichopria drosophilae (Hymenoptera: Diapriidae) to Control Drosophila suzukii (Diptera: Drosophilidae). Neotrop. Entomol. 48, 949–956. https://doi.org/10.1007/s13744-019-00707-3 (2019).CAS 
Article 
PubMed 

Google Scholar 
Rossi Stacconi, M. V., Grassi, A., Ioriatti, C. & Anfora, G. Augmentative releases of Trichopria drosophilae for the suppression of early season Drosophila suzukii populations. Biocontrol 64, 9–19. https://doi.org/10.1007/s10526-018-09914-0 (2018).CAS 
Article 

Google Scholar 
Poyet, M. et al. The wide potential trophic niche of the asiatic fruit fly Drosophila suzukii: The key of its invasion success in temperate Europe?. PLoS ONE 10, e0142785. https://doi.org/10.1371/journal.pone.0142785 (2015).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Mazzetto, F. et al. Drosophila parasitoids in northern Italy and their potential to attack the exotic pest Drosophila suzukii. J. Pest Sci. 89, 837–850. https://doi.org/10.1007/s10340-016-0746-7 (2016).Article 

Google Scholar 
Wang, X. G., Kacar, G., Biondi, A. & Daane, K. M. Foraging efficiency and outcomes of interactions of two pupal parasitoids attacking the invasive spotted wing drosophila. Biol. Control 96, 64–71. https://doi.org/10.1016/j.biocontrol.2016.02.004 (2016).Article 

Google Scholar 
Rossi Stacconi, M. V. et al. Host location and dispersal ability of the cosmopolitan parasitoid Trichopria drosophilae released to control the invasive spotted wing Drosophila. Biol. Control 117, 188–196. https://doi.org/10.1016/j.biocontrol.2017.11.013 (2018).Article 

Google Scholar 
Esteban-Santiago, J. M., Rodríguez-Leyva, E., Lomeli-Flores, J. R. & González-Cabrera, J. Demographic parameters of Trichopria drosophilae in three host species. Entomol. Exp. Appl. 169, 330–337. https://doi.org/10.1111/eea.13026 (2021).CAS 
Article 

Google Scholar 
Häussling, B. J. M., Lienenluke, J. & Stokl, J. The preference of Trichopria drosophilae for pupae of Drosophila suzukii is independent of host size. Sci. Rep. 11, 995. https://doi.org/10.1038/s41598-020-80355-5 (2021).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Wang, X. G., Kacar, G., Biondi, A. & Daane, K. M. Life-history and host preference of Trichopria drosophilae, a pupal parasitoid of spotted wing drosophila. Biocontrol 61, 387–397. https://doi.org/10.1007/s10526-016-9720-9 (2016).CAS 
Article 

Google Scholar 
Woltz, J. M. & Lee, J. C. Pupation behavior and larval and pupal biocontrol of Drosophila suzukii in the field. Biol. Control 110, 62–69. https://doi.org/10.1016/j.biocontrol.2017.04.007 (2017).Article 

Google Scholar 
Ballman, E. S., Collins, J. A. & Drummond, F. A. Pupation behavior and predation on Drosophila suzukii (Diptera: Drosophilidae) Pupae in Maine wild blueberry fields. J. Econ. Entomol. 110, 2308–2317. https://doi.org/10.1093/jee/tox233 (2017).Article 
PubMed 

Google Scholar 
Guillén, L., Aluja, M. N., Equihua, M. & Sivinski, J. Performance of two fruit fly (Diptera: Tephritidae) pupal parasitoids (Coptera haywardi [Hymenoptera: Diapriidae] and Pachycrepoideus vindemiae [Hymenoptera: Pteromalidae]) under different environmental soil conditions. Biol. Control 23, 219–227. https://doi.org/10.1006/bcon.2001.1011 (2002).Article 

Google Scholar 
Yi, C. et al. Life history and host preference of Trichopria drosophilae from Southern China, one of the effective pupal parasitoids on the Drosophila species. Insects 11, 103. https://doi.org/10.3390/insects11020103 (2020).Article 
PubMed Central 

Google Scholar 
BoychevaWoltering, S., Romeis, J. & Collatz, J. Influence of the rearing host on biological parameters of Trichopria drosophilae, a potential biological control agent of Drosophila suzukii. Insects. https://doi.org/10.3390/insects10060183 (2019).Article 

Google Scholar 
Otto, M. & Mackauer, M. The developmental strategy of an idiobiont ectoparasitoid, Dendrocerus carpenteri: Influence of variations in host quality on offspring growth and fitness. Oecologia 117, 353–364. https://doi.org/10.1007/s004420050668 (1998).ADS 
Article 
PubMed 

Google Scholar 
Bates, D., Machler, M., Bolker, B. M. & Walker, S. C. Fitting linear mixed-effects models using lme4. J. Stat. Softw. 67, 1–48. https://doi.org/10.18637/jss.v067.i01 (2015).Article 

Google Scholar 
R: A Language and Environment for Statistical Computing (Vienna, Austria, 2008).Johnson, S. N. & Gregory, P. J. Chemically-mediated host-plant location and selection by root-feeding insects. Physiol. Entomol. 31, 1–13. https://doi.org/10.1111/j.1365-3032.2005.00487.x (2006).CAS 
Article 

Google Scholar 
Bezerra Da Silva, C. S., Park, K. R., Blood, R. A. & Walton, V. M. Intraspecific competition affects the pupation behavior of spotted-wing drosophila (Drosophila suzukii). Sci. Rep. 9, 7775. https://doi.org/10.1038/s41598-019-44248-6 (2019).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Renkema, J. M. & Devkota, S. Pupation depth of spotted wing drosophila (Drosophila suzukii) and effects of field sanitation in Florida strawberries. Viii Int. Strawberry Symp. 1156, 849–855. https://doi.org/10.17660/ActaHortic.2017.1156.125 (2017).Article 

Google Scholar 
Tsitsipis, J. A. & Papanicolaou, E. P. Pupation depth in artificially reared olive fruits-flies Dacus-oleae (Diptera, Tephritidae), as affected by several physical characteristics of the substrates. Annales De Zoologie Ecologie Animale 11, 31–40 (1979).
Google Scholar 
Dimou, I., Koutsikopoulos, C., Economopoulos, A. P. & Lykakis, J. Depth of pupation of the wild olive fruit fly, Bactrocera (Dacus) oleae (Gmel.) (Dipt., Tephritidae), as affected by soil abiotic factors. J. Appl. Entomol. 127, 12–17. https://doi.org/10.1046/j.1439-0418.2003.00686.x (2003).Article 

Google Scholar 
de Belle, J. S., Hilliker, A. J. & Sokolowski, M. B. Genetic localization of foraging (for): A major gene for larval behavior in Drosophila melanogaster. Genetics 123, 157–163. https://doi.org/10.1093/genetics/123.1.157 (1989).Article 
PubMed 

Google Scholar 
Sokolowski, M. B. et al. Ecological genetics and behaviour of Drosophila melanogaster larvae in nature. Anim. Behav. 34, 403–408. https://doi.org/10.1016/S0003-3472(86)80109-9 (1986).Article 

Google Scholar 
Rodriguez, L., Sokolowski, M. B. & Shore, J. S. Habitat selection by Drosophila melanogaster larvae. J. Evol. Biol. 5, 61–70. https://doi.org/10.1046/j.1420-9101.1992.5010061.x (1992).Article 

Google Scholar 
McIntosh, H., Atucha, A., Townsend, P. A., Hills, W. B. & Guédot, C. Plastic mulches reduce adult and larval populations of Drosophila suzukii in fall-bearing raspberry. J. Pest. Sci. 95, 525–536. https://doi.org/10.1007/s10340-021-01456-2 (2021).Article 

Google Scholar 
Ballman, E. & Drummond, F. Larval movement of spotted wing drosophila, Drosophila suzukii (Matsumura) (Diptera: Drosophilidae). J. Kansas Entomol. Soc. 92, 412–421. https://doi.org/10.2317/0022-8567-92.1.412 (2019).Article 

Google Scholar  More

In this study, the five main characteristics of green spaces that were measured were area, perimeter, perimeter-area ratio, leaf area index, and canopy density. The structure of parameter between them is shown in Table 3.Table 3 Parameter structure of the cooling and humidification effect based on the spatial characteristics of green spaces.Full size tableCorrelation between various spatial characteristics and cooling and humidifying intensity in green spacesSmall-size green spacesFigures 4 and 6 shows the results of linear regressions between spatial characteristics and the cooling effect in small-size green spaces. There were relatively weak correlations between area, perimeter, perimeter-area ratio, leaf area index and cooling intensity, and a strong correlation between canopy density and cooling intensity. Small-size green space has the weakest positive correlation between perimeter-area ratio and cooling intensity (R2 = 0.11), and its canopy density and cooling intensity have the strongest positive correlation (R2 = 0.64). Meanwhile, small-size green space has weakest negative correlation between perimeter and humidifying intensity (R2 = 0.17), and its leaf area index and humidifying intensity have significant positive correlation (R2 = 0.42). Figures 4a and 5a show that for every 1 ha increase in area of small-size green spaces, the cooling intensity increased by 1.026 °C, and the humidifying intensity decreased by 1.56%. Figures 4b and 5b show that for every 100 m increase in perimeter, the cooling intensity decreases by 1.06 °C, and the humidifying intensity decreased by 1.19%. Figures 4c and 5c show that for every 0.01 increase in the perimeter-area ratio, the cooling intensity increases by 1.12 °C, and the humidifying intensity increased by 1.46%. Figures 4d and 5d show that for every 0.1 increase in the leaf area index, the cooling intensity increases by 1.11 °C, and the humidifying intensity increased by 1.12%. Figures 4e and 5e show that each 0.01 increase in the canopy density, the cooling intensity increases by 1.60 °C, and each 0.1 increase in canopy density, the humidifying intensity increased by 1.15% (Fig. 6).
Figure 4Linear regressions between spatial characteristics and cooling intensity of small-size green spaces.Full size imageFigure 5Linear regressions of spatial characteristics and humidifying intensity of small-size green spaces.Full size imageFigure 6The correlation between the spatial characteristics of small-size green spaces and the intensity of cooling and humidifying (GA means green area; GP means green perimeter; GPAR means green perimeter-area ratio; LAI means leaf area index; CD means canopy density).Full size imageMedium-size green spacesFigures 7 and 9 shows the linear regressions between spatial characteristics and cooling intensity in medium-size green spaces. There was an extremely significant positive correlation between area and cooling intensity, an insignificant positive correlation between the leaf area index and cooling intensity, and a relatively weak negative correlation between the other three characteristics and cooling intensity. Medium-size green space has the weakest negative correlation between canopy density and cooling intensity (R2 = 0.12), and its green area and cooling intensity have the strongest positive correlation (R2 = 0.83). Meanwhile, medium-size green space has weakest negative correlation between perimeter-area ratio and humidifying intensity (R2 = 0.41), and its area and humidifying intensity have most significant positive correlation (R2 = 0.81). Figures 7a and 8a show that for every 1 ha increase in area of medium-size green spaces, the cooling intensity increased by 1.19 °C, and the humidifying intensity increased by 1.24%. Figures 7b and 8b show that for every 100 m increase in perimeter, the cooling intensity decreases by 1.02 °C, and the humidifying intensity increased by 1.17%. Figures 7c and 8c show that for every 0.01 increase in the perimeter-area ratio, the cooling intensity decreases by 1.29 °C, and the humidifying intensity decreased by 2.40%. Figures 7d and 8d show that for every 0.1 increase in the leaf area index, the cooling intensity increases by 1.37 °C, and the humidifying intensity decreased by 1.92%. Figures 7e and 8e show that each 0.01 increase in the canopy density, increases the cooling intensity decreases by 1.23 °C, and the humidifying intensity decreased by 6.48% (Fig. 9).Figure 7Linear regressions between spatial characteristics and cooling intensity of medium-size green spaces.Full size imageFigure 8Linear regressions of spatial characteristics and humidifying intensity of medium-size green spaces.Full size imageFigure 9The correlation between the spatial characteristics of medium-size green spaces and the intensity of cooling and humidifying (GA means green area; GP means green perimeter; GPAR means green perimeter-area ratio; LAI means leaf area index; CD means canopy density).Full size imageLarge-size green spacesFigures 10 and 12 shows the linear regressions between spatial characteristics and cooling intensity in large-size green spaces. There was an insignificant correlation between area and cooling intensity, a weak correlation between canopy density and cooling intensity, and a significant correlation between perimeter, perimeter-area ratio and the leaf area index and cooling intensity. Medium-size green space has the weakest negative correlation between green area and cooling intensity (R2 = 0.35), and its leaf area index and cooling intensity have the strongest positive correlation (R2 = 0.92). Meanwhile, medium-size green space has weakest negative correlation between perimeter-area ratio and humidifying intensity (R2 = 0.11), and its leaf area index and humidifying intensity have most significant positive correlation (R2 = 0.39). Figures 10a and 11a show that for every 1 ha increase in area of large-size green spaces, the cooling intensity decreased by 1.02 °C, and the humidifying intensity decreased by 1.22%. Figures 10b and 11b show that for every 100 m increase in perimeter, the cooling intensity decreases by 1.05 °C, and the humidifying intensity decreased by 1.34%. Figures 10c and 11c show that for every 0.005 increase in the perimeter-area ratio, the cooling intensity decreases by 1.43 °C, and each 0.01 increase in perimeter-area ratio, the humidifying intensity decreased by 1.27%. Figures 10d and 11d show that for every 0.1 increase in the leaf area index, the cooling intensity increases by 2.41 °C, and the humidifying intensity increased by 1.37%. Figures 10e and 11e show that each 0.1 increase in the canopy density, the cooling intensity increased by 3.69 °C, and the humidifying intensity decreased by 2.84% (Fig. 12).Figure 10Linear regressions of spatial characteristics and cooling intensity of large-size green spaces.Full size imageFigure 11Linear regressions of spatial characteristics and humidifying intensity of large-size green spaces.Full size imageFigure 12The correlation between the spatial characteristics of large-size green spaces and the intensity of cooling and humidifying (GA means green area; GP means green perimeter; GPAR means green perimeter-area ratio; LAI means leaf area index; CD means canopy density).Full size imageQuantitative analysis of the microclimatic effects of different types of green spacesQuantitative analysis of the effects of different types of green space on cooling intensityFigure 13 shows the linear regressions between the different types of green spaces and cooling intensity. There were negative correlations between green spaces a short, medium, and long distance from a water body and cooling intensity in small-size green spaces, medium-size green spaces and large-size green spaces. The negative correlation between the distance to a water body and cooling intensity in medium-size green spaces was most significant (R2 = 0.985). The greater the distance to a water body, the lower the cooling intensity. For medium-size green spaces, for every 1/4 increase in the distance ratio, the cooling intensity decreased by 0.81 °C. For small-size green spaces, for every 1/4 increase in the distance ratio, the cooling intensity decreased by 1.04 °C. For large-size green spaces, for every 1/4 increase in the distance ratio, the cooling intensity decreased by 1.36 °C. For small-, medium-, and large-size green spaces, there was a positive correlation between canopy density and cooling intensity. There was a most significant positive correlation between canopy density and cooling intensity in large-size green spaces (R2 = 0.941). The greater the canopy density, the greater the cooling intensity. For large green spaces, for every 0.5 increase in canopy density, the cooling intensity increased by 0.16 °C. For small-size green spaces, for every 0.5 increase in canopy density, the cooling effect increased by 0.15 °C. For medium-size green spaces, for every 0.5 increase in canopy density, the cooling intensity increased by 0.16 °C.Figure 13Linear regressions between the distance from different types of green spaces to water areas, canopy density and cooling intensity.Full size imageQuantitative analysis of the effects of different types of green space on humidifying intensityFigure 14 shows the linear regression between the distance of a green space from a water body, canopy density and humidifying intensity. There was a negative correlation between the distance to a water body and humidifying intensity in small, medium, and large green spaces. The negative correlation between the distance to a water body and humidifying intensity in small green spaces was most significant (R2 = 0.996). The longer the distance, the lower the humidifying intensity. For small green spaces, for every 1/4 in-crease in the distance ratio, the humidifying intensity decreased by 4.23%. For medium-size green spaces, for every 1/4 increase in the distance ratio, the humidifying intensity decreased by 3.02%. For large-size green spaces, for every 1/4 increase in the distance ratio, the humidifying intensity de-creased by 6.14%. For small, medium, and large green spaces, there was a positive correlation between canopy density and humidifying intensity. The positive correlation between canopy density and humidifying intensity in medium-size green spaces was extremely significant (R2 = 0.925). The greater the canopy density, the greater the humidifying intensity. For medium-size green spaces, for every 0.5 increase in canopy density, the humidifying intensity increased by 3.29%. For small-size green spaces, for every 0.5 increase in canopy density, the humidifying intensity increased by 3.17%. For large-size green spaces, for every 0.5 increase in canopy density, the humidifying intensity increased by 4.06% (Fig. 15).
Figure 14Linear regressions between the distance from different types of green space to water area, canopy density and humidifying intensity.Full size imageFigure 15Correlation of different green space types with water distance, canopy density and cooling and humidifying intensity.Full size imageEffect of shape and area of water bodies on microclimatic effects based on numerical simulationBanded waterWe constructed a numerical simulation model to explore the effects of a simulated increase in water body area on cooling and humidification. Figure 16 shows the simulated distribution characteristics of temperature and relative humidity after a 5% and 10% increase in water area at 14:00 when temperatures were high. The results suggest that between 7:00 and 10:00, with a 5% and 10% increase in water area, the air temperature was basically the same and the cooling effect was insignificant. However, between 12:00 and 19:00 and particularly in the hours between 13:00 and 16:00 when temperatures were highest, a 5% increase in water area produced a significant cooling effect, with a daily average value of 0.05 °C and a maximum value of 0.09 °C. A 10% increase in water area produced an extremely significant cooling effect, with a daily average value of 0.07 °C and a maximum value of 0.14 °C. From 11:00 to 19:00, a 5% increase in water area produced a significant humidifying effect, with a daily average value of 0.08% and a maximum value of 0.17%. A 10% increase produced an extremely significant humidifying effect, with a daily average value of 0.13% and a maximum value of 0.26% (See supplementary file).Figure 16Distribution characteristics of cooling and humidifying effects of simulated increase of banded water area at 14:00. (a) original cooling effect of banded water in the sample area; (b) cooling effect of 5% increase in water area; (c) cooling effect of 10% increase in water area; (d) original humidifying effect of banded water in the sample area; (e) humidifying effect of 5% increase in water area; (f) humidifying effect of 10% increase of water area.Full size imageMassive waterFigure 17 shows the simulated distribution characteristics of the cooling and humidifying effects after a 5% and 10% increase in the water area at 14:00 when temperatures were high. Between 8:00 and 19:00, a 5% and 10% increase in water area produced a significant cooling effect. At 19:00, the numerical simulation result was abnormal when the water area increased by 5% and 10%; at 13:00, the numerical simulation result was also ab-normal when the water area increased by 10%. After excluding the abnormal simulated data, a 5% increase in water area produced a cooling effect, with a daily average value of 0.06 °C and a maximum value of 0.10 °C. A 10% increase in water area produced an extremely significant cooling effect, with a daily average value of 0.10 °C and a maximum value of 0.18 °C. Between 11:00 and 19:00, a 5% increase in water area produced a significant humidifying effect, with a daily average value of 0.05% and a maximum value of 0.13%. A 10% increase in water area produced an extremely significant humidifying effect, with a daily average value of 0.13% and a maximum value of 0.27% (See supplementary file).Figure 17Distribution characteristics of cooling and humidifying effects of simulated increase of massive water area at 14:00. (a) original cooling effect of massive water in the sample area; (b) cooling effect of 5% increase in water area; (c) cooling effect of 10% increase in water area; (d) original humidifying effect of massive water in the sample area; (e) humidifying effect of 5% increase in water area; (f) humidifying effect of 10% increase of water area.Full size imageAnnular waterFigure 18 shows the simulated distribution characteristics of the cooling and humidifying effects after a 5% and 10% increase in the area of the annular water body at 14:00 when temperatures were high. Between 7:00 and 19:00, a 5% and 10% increase in water area produced a significant cooling effect. Between 11:00 and 16:00 when temperatures were high, a 5% increase in water area produced a cooling effect, with a daily average value of 0.06 °C and a maximum value of 0.14 °C°C and a 10% increase in water area produced an extremely significant cooling effect, with a daily average value of 0.13 °C and a maximum value of 0.28 °C. Between 7:00 and 19:00, a 5% and 10% increase in water area produced significant humidifying effects. Between 11:00 and 16:00 when temperatures were high, a 5% increase in water area produced an extremely significant humidifying effect, with a daily average value of 0.17% and a maximum value of 0.39% and a 10% increase in water area produced an extremely significant humidifying effect with a daily average value of 0.38% and a maximum value of 0.81% (See supplementary file).Figure 18Distribution characteristics of cooling and humidifying effects of simulated increase of annular water area at 14:00. (a) original cooling effect of annular water in the sample area; (b) cooling effect of 5% increase in water area; (c) cooling effect of 10% increase in water area; (d) original humidifying effect of annular water in the sample area; (e) humidifying effect of 5% increase in water area; (f) humidifying effect of 10% increase of water area.Full size image More

Canadell, J. G. et al. Multi-decadal increase of forest burned area in Australia is linked to climate change. Nat. Commun. 12, 6921 (2021).CAS 
Article 

Google Scholar 
Jolly, W. M. et al. Climate-induced variations in global wildfire danger from 1979 to 2013. Nat. Commun. 6, 7537 (2015).CAS 
Article 

Google Scholar 
Jain, P., Wang, X. & Flannigan, M. D. Trend analysis of fire season length and extreme fire weather in North America between 1979 and 2015. Int. J. Wildland Fire 26, 1009–1020 (2017).Article 

Google Scholar 
Wotton, B. M. & Flannigan, M. D. Length of the fire season in a changing climate. For. Chronicle 69, 187–192 (1993).
Google Scholar 
Collins, L. et al. The 2019/2020 mega-fires exposed Australian ecosystems to an unprecedented extent of high-severity fire. Environ. Res. Lett. 16, 044029 (2021).Article 

Google Scholar 
Higuera, P. E. & Abatzoglou, J. T. Record‐setting climate enabled the extraordinary 2020 fire season in the western United States. Glob. Change Biol. 27, 1–2 (2021).Article 

Google Scholar 
Nolan, R. H. et al. Limits to post-fire vegetation recovery under climate change. Plant Cell Environ. 44, 3471–3489 (2021).CAS 
Article 

Google Scholar 
Abram, N. J. et al. Connections of climate change and variability to large and extreme forest fires in southeast Australia. Commun. Earth Environ. 2, 8 (2021).Article 

Google Scholar 
Dickman, C. R. Ecological consequences of Australia’s “Black Summer” bushfires: managing for recovery. Int. Environ. Assess. Manag. 17, 1162–1167 (2021).Article 

Google Scholar 
Swain, D. L. A shorter, sharper rainy season amplifies California wildfire risk. Geophys. Res. Lett. 48, e2021GL092843 (2021).
Google Scholar 
Keeley, J. E. Fire intensity, fire severity and burn severity: a brief review and suggested usage. Int. J. Wildland Fire 18, 116–126 (2009).Article 

Google Scholar 
He, T., Pausas, J. G., Belcher, C. M., Schwilk, D. W. & Lamont, B. B. Fire-adapted traits of Pinus arose in the fiery Cretaceous. New Phytol. 194, 751–759 (2012).Article 

Google Scholar 
Bradstock, R. A. A biogeographic model of fire regimes in Australia: current and future implications. Glob. Ecol. Biogeogr. 19, 145–158 (2010).Article 

Google Scholar 
Bowman, D. M., Murphy, B. P., Neyland, D. L., Williamson, G. J. & Prior, L. D. Abrupt fire regime change may cause landscape-wide loss of mature obligate seeder forests. Glob. Change Biol. 20, 1008–1015 (2014).Article 

Google Scholar 
Keeley, J. E., Pausas, J. G., Rundel, P. W., Bond, W. J. & Bradstock, R. A. Fire as an evolutionary pressure shaping plant traits. Trends Plant Sci. 16, 406–411 (2011).CAS 
Article 

Google Scholar 
Barrett, K. et al. Postfire recruitment failure in Scots pine forests of southern Siberia. Remote Sens. Environ. 237, 111539 (2020).Article 

Google Scholar 
Miller, R. G., Fontaine, J. B., Merritt, D. J., Miller, B. P. & Enright, N. J. Experimental seed sowing reveals seedling recruitment vulnerability to unseasonal fire. Ecol. Appl. 31, e02411 (2021).
Google Scholar 
Prior, L. D., Williamson, G. J. & Bowman, D. M. Impact of high-severity fire in a Tasmanian dry eucalypt forest. Austral. J. Bot. 64, 193–205 (2016).Article 

Google Scholar 
Brewer, J. S. Long-term population changes of a fire-adapted plant subjected to different fire seasons. Nat. Areas J. 26, 267–273 (2006).Article 

Google Scholar 
Keith, D. A., Holman, L., Rodoreda, S., Lemmon, J. & Bedward, M. Plant functional types can predict decade‐scale changes in fire‐prone vegetation. J. Ecol. 95, 1324–1337 (2007).Article 

Google Scholar 
Savage, M., Mast, J. N. & Feddema, J. J. Double whammy: high-severity fire and drought in ponderosa pine forests of the Southwest. Can. J. For. Res. 43, 570–583 (2013).Article 

Google Scholar 
Miller, R. G. et al. Mechanisms of fire seasonality effects on plant populations. Trends Ecol. Evol. 34, 1104–1117 (2019).Article 

Google Scholar 
Tangney, R., Merritt, D. J., Fontaine, J. B. & Miller, B. P. Seed moisture content as a primary trait regulating the lethal temperature thresholds of seeds. J. Ecol. 107, 1093–1105 (2019).Article 

Google Scholar 
Tangney, R. et al. Seed dormancy interacts with fire seasonality mechanisms. Trends Ecol. Evol. 35, 1057–1059 (2020).Article 

Google Scholar 
Bowman, D. M. et al. The human dimension of fire regimes on Earth. J. Biogeogr. 38, 2223–2236 (2011).Article 

Google Scholar 
Knapp, E. E., Estes, B. L. & Skinner, C. N. Ecological effects of prescribed fire season: a literature review and synthesis for managers. Gen. Tech. Rep. https://doi.org/10.2737/PSW-GTR-224 (2009).Miller, R. G. et al. Fire seasonality mechanisms are fundamental for understanding broader fire regime effects. Trends Ecol. Evol. 35, 869–871 (2020).Article 

Google Scholar 
Keeley, J. E. & Syphard, A. D. Twenty-first century California, USA, wildfires: fuel-dominated vs. wind-dominated fires. Fire Ecol. 15, 24 (2019).Article 

Google Scholar 
Lamont, B. B., Enright, N. J. & He, T. Fitness and evolution of resprouters in relation to fire. Plant Ecol. 212, 1945–1957 (2011).Article 

Google Scholar 
Pausas, J. G. & Bradstock, R. A. Fire persistence traits of plants along a productivity and disturbance gradient in mediterranean shrublands of south‐east Australia. Glob. Ecol. Biogeogr. 16, 330–340 (2007).Article 

Google Scholar 
Pausas, J. G. & Keeley, J. E. Evolutionary ecology of resprouting and seeding in fire-prone ecosystems. New Phytol. 204, 55–65 (2014).Article 

Google Scholar 
Fairman, T. A., Bennett, L. T. & Nitschke, C. R. Short-interval wildfires increase likelihood of resprouting failure in fire-tolerant trees. J. Environ. Manag. 231, 59–65 (2019).Article 

Google Scholar 
Pyke, G. H. Fire-stimulated flowering: a review and look to the future. Critic. Rev. Plant Sci. 36, 179–189 (2017).Article 

Google Scholar 
Zirondi, H. L., Ooi, M. K. J. & Fidelis, A. Fire-triggered flowering is the dominant post-fire strategy in a tropical savanna. J. Veg. Sci. 32, e12995 (2021).Article 

Google Scholar 
Howe, H. F. Response of Zizia aurea to seasonal mowing and fire in a restored Prairie. Am. Midl. Nat. 141, 373–380 (1999).Article 

Google Scholar 
Thompson, K. Seeds and seed banks. New Phytol. 106, 23–34 (1987).Article 

Google Scholar 
Baskin, C. C. & Baskin, J. M. Seeds: Ecology, Biogeography, and Evolution of Dormancy and Germination 2nd edn (Academic Press, 2001).Alvarado, V. & Bradford, K. J. A hydrothermal time model explains the cardinal temperatures for seed germination. Plant Cell Environ. 25, 1061–1069 (2002).Article 

Google Scholar 
Mackenzie, B. D. E., Auld, T. D., Keith, D. A., Hui, F. K. C. & Ooi, M. K. J. The effect of seasonal ambient temperatures on fire-stimulated germination of species with physiological dormancy: a case study using boronia (Rutaceae). PLoS One 11, e0156142 (2016).Article 
CAS 

Google Scholar 
Ooi, M. K. J. Delayed emergence and post-fire recruitment success: effects of seasonal germination, fire season and dormancy type. Austral. J. Bot. 58, 248–256 (2010).Article 

Google Scholar 
Bond, W. Fire survival of Cape Proteaceae-influence of fire season and seed predators. Vegetatio 56, 65–74 (1984).Article 

Google Scholar 
Keith, D. A., Dunker, B. & Driscoll, D. A. Dispersal: the eighth fire seasonality effect on plants. Trends Ecol. Evol. 35, 305–307 (2020).Article 

Google Scholar 
Paroissien, R. & Ooi, M. K. J. Effects of fire season on the reproductive success of the post-fire flowerer Doryanthes excelsa. Environ. Exp. Bot. 192, 104634 (2021).Article 

Google Scholar 
Furlaud, J. M., Prior, L. D., Williamson, G. J. & Bowman, D. M. J. S. Bioclimatic drivers of fire severity across the Australian geographical range of giant Eucalyptus forests. J. Ecol. 109, 2514–2536 (2021).Article 

Google Scholar 
Thomsen, A. M. & Ooi, M. K. J. Shifting season of fire and its interaction with fire severity: Impacts on reproductive effort in resprouting plants. Ecol. Evol. 12, e8717 (2022).Article 

Google Scholar 
Fill, J. M. & Crandall, R. M. Stronger evidence needed for global fire season effects. Trends Ecol. Evol. 35, 867–868 (2020).Article 

Google Scholar 
Jump, A. S. & Peñuelas, J. Running to stand still: adaptation and the response of plants to rapid climate change. Ecol. Lett. 8, 1010–1020 (2005).Article 

Google Scholar 
Inouye, D. W. Climate change and phenology. Wiley Interdiscip. Rev. Clim. Change n/a, e764 (2022).
Google Scholar 
Enright, N. J., Marsula, R., Lamont, B. B. & Wissel, C. The ecological significance of canopy seed storage in fire-prone environments: a model for non-sprouting shrubs. J. Ecol. 86, 946–959 (1998).Article 

Google Scholar 
Setterfield, S. A. The impact of experimental fire regimes on seed production in two tropical eucalypt species in northern Australia. Austral. J. Ecol. 22, 279–287 (1997).Article 

Google Scholar 
Collette, J. C. & Ooi, M. K. J. Evidence for physiological seed dormancy cycling in the woody shrub Asterolasia buxifolia and its ecological significance in fire-prone systems. Plant Biol. 22, 745–749 (2020).CAS 
Article 

Google Scholar 
Setterfield, S. A. Seedling establishment in an Australian tropical savanna: effects of seed supply, soil disturbance and fire. J. Appl. Ecol. 39, 949–959 (2002).Article 

Google Scholar 
Russell-Smith, J. & Edwards, A. C. Seasonality and fire severity in savanna landscapes of monsoonal northern Australia. Int. J. Wildland Fire 15, 541–550 (2006).Article 

Google Scholar 
Whitehead, P. J., Purdon, P., Russell-Smith, J., Cooke, P. M. & Sutton, S. The management of climate change through prescribed Savanna burning: Emerging contributions of indigenous people in Northern Australia. Public Adm. Dev. 28, 374–385 (2008).Article 

Google Scholar 
Prior, L. D., Williams, R. J. & Bowman, D. M. Experimental evidence that fire causes a tree recruitment bottleneck in an Australian tropical savanna. J. Tropical Ecol. 26, 595–603 (2010).Abatzoglou, J. T., Williams, A. P. & Barbero, R. Global emergence of anthropogenic climate change in fire weather indices. Geophys. Res. Lett. 46, 326–336 (2019).Article 

Google Scholar 
Ferreira, L. N., Vega-Oliveros, D. A., Zhao, L., Cardoso, M. F. & Macau, E. E. N. Global fire season severity analysis and forecasting. Comput. Geosci. 134, 104339 (2020).Article 

Google Scholar 
Flannigan, M. et al. Global wildland fire season severity in the 21st century. For. Ecol. Manag. 294, 54–61 (2013).Article 

Google Scholar 
Ansley, R. J. & Castellano, M. J. Prickly pear cactus responses to summer and winter fires. Rangel. Ecol. Manag. 60, 244–252 (2007).Article 

Google Scholar 
Ansley, R. J., Kramp, B. A. & Jones, D. L. Honey mesquite (Prosopis glandulosa) seedling responses to seasonal timing of fire and fireline intensity. Rangel. Ecol. Manag. 68, 194–203 (2015).Article 

Google Scholar 
Armstrong, G. & Legge, S. The post-fire response of an obligate seeding Triodia species (Poaceae) in the fire-prone Kimberley, north-west Australia. Int. J. Wildland Fire 20, 974–981 (2012).Article 

Google Scholar 
Bellows, R. S., Thomson, A. C., Helmstedt, K. J., York, R. A. & Potts, M. D. Damage and mortality patterns in young mixed conifer plantations following prescribed fires in the Sierra Nevada, California. For. Ecol. Manag. 376, 193–204 (2016).Article 

Google Scholar 
Beyers, J. L. & Wakeman, C. D. Season of burn effects in southern California chaparral. In Second interface between ecology and land development in California 45–55 (Occidental College, CA, 2000).Bowen, B. J. & Pate, J. S. Effect of season of burn on shoot recovery and post‐fire flowering performance in the resprouter Stirlingia latifolia R. Br.(Proteaceae). Austral Ecol. 29, 145–155 (2004).Article 

Google Scholar 
Casals, P., Valor, T., Rios, A. & Shipley, B. Leaf and bark functional traits predict resprouting strategies of understory woody species after prescribed fires. For. Ecol. Manag. 429, 158–174 (2018).Article 

Google Scholar 
Céspedes, B., Torres, I., Luna, B., Pérez, B. & Moreno, J. M. Soil seed bank, fire season, and temporal patterns of germination in a seeder-dominated Mediterranean shrubland. Plant Ecol. 213, 383–393 (2012).Article 

Google Scholar 
Clabo, D. C. & Clatterbuck, W. K. Shortleaf pine (Pinus echinata, Pinaceae) seedling sprouting responses: Clipping and burning effects at various seedling ages and seasons. J. Torrey Bot. Soc. 146, 96–110 (2019).Article 

Google Scholar 
Drewa, P. B. Effects of fire season and intensity on Prosopis glandulosa Torr. var. glandulosa. Int. J. Wildland Fire 12, 147–157 (2003).Article 

Google Scholar 
Drewa, P. B., Platt, W. J. & Moser, E. B. Fire effects on resprouting of shrubs in headwaters of southeastern longleaf pine savannas. Ecology 83, 755–767 (2002).Article 

Google Scholar 
Drewa, P. B., Thaxton, J. M. & Platt, W. J. Responses of root‐crown bearing shrubs to differences in fire regimes in Pinus palustris (longleaf pine) savannas: exploring old‐growth questions in second‐growth systems. Appl. Veg. Sci. 9, 27–36 (2006).
Google Scholar 
Ellsworth, L. M. & Kauffman, J. B. Seedbank responses to spring and fall prescribed fire in mountain big sagebrush ecosystems of differing ecological condition at Lava Beds National Monument, California. J. Arid Environ. 96, 1–8 (2013).Article 

Google Scholar 
Fairfax, R. et al. Effects of multiple fires on tree invasion in montane grasslands. Landsc. Ecol. 24, 1363–1373 (2009).Article 

Google Scholar 
Fill, J. M., Welch, S. M., Waldron, J. L. & Mousseau, T. A. The reproductive response of an endemic bunchgrass indicates historical timing of a keystone process. Ecosphere 3, 1–12 (2012).Article 

Google Scholar 
Grant, C. Post-burn vegetation development of rehabilitated bauxite mines in Western Australia. For. Ecol. Manag. 186, 147–157 (2003).Article 

Google Scholar 
Hajny, K. M., Hartnett, D. C. & Wilson, G. W. Rhus glabra response to season and intensity of fire in tallgrass prairie. Int. J. Wildland Fire 20, 709–720 (2011).Article 

Google Scholar 
Holmes, P. A comparison of the impacts of winter versus summer burning of slash fuel in alien-invaded fynbos areas in the Western Cape. Southern African For. J. 192, 41–50 (2001).Article 

Google Scholar 
Jasinge, N., Huynh, T. & Lawrie, A. Consequences of season of prescribed burning on two spring-flowering terrestrial orchids and their endophytic fungi. Austr. J. Bot. 66, 298–312 (2018).Article 

Google Scholar 
Jasinge, N., Huynh, T. & Lawrie, A. Changes in orchid populations and endophytic fungi with rainfall and prescribed burning in Pterostylis revoluta in Victoria, Australia. Ann. Bot. 121, 321–334 (2018).CAS 
Article 

Google Scholar 
Kauffman, J. & Martin, R. Sprouting shrub response to different seasons and fuel consumption levels of prescribed fire in Sierra Nevada mixed conifer ecosystems. For. Sci. 36, 748–764 (1990).
Google Scholar 
Keyser, T. L., Greenberg, C. H. & McNab, W. H. Season of burn effects on vegetation structure and composition in oak-dominated Appalachian hardwood forests. For. Ecol. Manag. 433, 441–452 (2019).Article 

Google Scholar 
Knox, K. & Clarke, P. J. Fire season and intensity affect shrub recruitment in temperate sclerophyllous woodlands. Oecologia 149, 730–739 (2006).CAS 
Article 

Google Scholar 
Lamont, B. B. & Downes, K. S. Fire-stimulated flowering among resprouters and geophytes in Australia and South Africa. Plant Ecol. 212, 2111–2125 (2011).Article 

Google Scholar 
Lesica, P. & Martin, B. Effects of prescribed fire and season of burn on recruitment of the invasive exotic plant, Potentilla recta, in a semiarid grassland. Restoration Ecol. 11, 516–523 (2003).Article 

Google Scholar 
Moreno, J. M. et al. Rainfall patterns after fire differentially affect the recruitment of three Mediterranean shrubs. Biogeosciences 8, 3721–3732 (2011).Article 

Google Scholar 
Mulligan, M. K. & Kirkman, L. K. Burning influences on wiregrass (Aristida beyrichiana) restoration plantings: natural seedling recruitment and survival. Restor. Ecol. 10, 334–339 (2002).Article 

Google Scholar 
Nield, A. P., Enright, N. J. & Ladd, P. G. Fire-stimulated reproduction in the resprouting, non-serotinous conifer Podocarpus drouynianus (Podocarpaceae): the impact of a changing fire regime. Popul. Ecol. 58, 179–187 (2016).Article 

Google Scholar 
Norden, A. H. & Kirkman, L. K. Persistence and prolonged winter dormancy of the federally endangered Schwalbea Americana L.(Scrophulariaceae) following experimental management techniques. Nat. Areas J. 24, 129–134 (2004).
Google Scholar 
Olson, M. S. & Platt, W. J. Effects of habitat and growing season fires on resprouting of shrubs in longleaf pine savannas. Vegetatio 119, 101–118 (1995).Article 

Google Scholar 
Ooi, M. K. The importance of fire season when managing threatened plant species: a long-term case-study of a rare Leucopogon species (Ericaceae). J. Environ. Manag. 236, 17–24 (2019).Article 

Google Scholar 
Pavlovic, N. B., Leicht-Young, S. A. & Grundel, R. Short-term effects of burn season on flowering phenology of savanna plants. Plant Ecology 212, 611–625 (2011).Article 

Google Scholar 
Payton, I. J. & Pearce, H. G. Fire-Induced Changes to the Vegetation of Tall-Tussock (Chionochloa rigida) Grassland Ecosystems. (Department of Conservation Wellington, New Zealand, 2009).Peguero, G. & Espelta, J. M. Disturbance intensity and seasonality affect the resprouting ability of the neotropical dry-forest tree Acacia pennatula: do resources stored below-ground matter? J. Tropical Ecol. 28, 539–546 (2011).Risberg, L. & Granström, A. Exploiting a window in time. Fate of recruiting populations of two rare fire-dependent Geranium species after forest fire. Plant Ecol. 215, 613–624 (2014).Article 

Google Scholar 
Rodríguez-Trejo, D. A., Castro-Solis, U. B., Zepeda-Bautista, M. & Carr, R. J. First year survival of Pinus hartwegii following prescribed burns at different intensities and different seasons in central Mexico. Int. J. Wildland Fire 16, 54–62 (2007).Article 

Google Scholar 
Russell, M., Vermeire, L., Ganguli, A. & Hendrickson, J. Fire return interval and season of fire alter bud banks. Rangel. Ecol. Manag.72, 542–550 (2019).Article 

Google Scholar 
Russell-Smith, J., Whitehead, P. J., Cook, G. D. & Hoare, J. L. Response of Eucalyptus‐dominated savanna to frequent fires: lessons from Munmarlary, 1973–1996. Ecol. Monogr. 73, 349–375 (2003).Article 

Google Scholar 
Schmidt, I. B., Sampaio, A. B. & Borghetti, F. Effects of the season on sexual reproduction and population structure of Heteropterys pteropetala (Adr. Juss.), Malpiguiaceae, in areas of Cerrado sensu stricto submitted to biennial fires. Acta Bot. Brasilica 19, 927–934 (2005).Article 

Google Scholar 
Shepherd, B. J., Miller, D. L. & Thetford, M. Fire season effects on flowering characteristics and germination of longleaf pine (Pinus palustris) savanna grasses. Restor. Ecol. 20, 268–276 (2012).Article 

Google Scholar 
Spier, L. P. & Snyder, J. R. Effects of wet-and dry-season fires on Jacquemontia curtisii, a south Florida pine forest endemic. Nat. Areas J. 18, 350–357 (1998).
Google Scholar 
Tsafrir, A. et al. Fire season modifies the perennial plant community composition through a differential effect on obligate seeders in eastern Mediterranean woodlands. Appl. Veg. Sci. 22, 115–126 (2019).Article 

Google Scholar 
Vander Yacht, A. L. et al. Vegetation response to canopy disturbance and season of burn during oak woodland and savanna restoration in Tennessee. For. Ecol. Manag. 390, 187–202 (2017).Article 

Google Scholar 
Vidaller, C., Dutoit, T., Ramone, H. & Bischoff, A. Fire increases the reproduction of the dominant grass Brachypodium retusum and Mediterranean steppe diversity in a combined burning and grazing experiment. Appl. Veg. Sci. 22, 127–137 (2019).Article 

Google Scholar 
Williams, P. R., Congdon, R. A., Grice, A. C. & Clarke, P. J. Soil temperature and depth of legume germination during early and late dry season fires in a tropical eucalypt savanna of north‐eastern Australia. Austral Ecol. 29, 258–263 (2004).Article 

Google Scholar 
Williams, P. R., Congdon, R. A., Grice, A. C. & Clarke, P. J. Germinable soil seed banks in a tropical savanna: seasonal dynamics and effects of fire. Austral Ecol. 30, 79–90 (2005).Article 

Google Scholar 
Zhao, H. et al. Ecophysiological influences of prescribed burning on wetland plants: a case study in Sanjiang Plain wetlands, northeast China. Fresenius Environ. Bull 20, 2932–2938 (2011).CAS 

Google Scholar 
Pick, J. L., Nakagawa, S. & Noble, D. W. Reproducible, flexible and high‐throughput data extraction from primary literature: The metaDigitise r package. Methods in Ecol. Evol. 10, 426–431 (2019).Article 

Google Scholar 
Team, R. C. R: A Language and Environment for Statistical Computing. https://www.R-project.org/ (2020).Higgins, J. P. et al. Cochrane Handbook for Systematic Reviews of Interventions. (John Wiley & Sons, 2019).Lüdecke, D., Lüdecke, M. D. & David, B. W. Package ‘esc’. https://strengejacke.github.io/esc (2017).Schwarzer, G. meta: An R package for meta-analysis. R news 7, 40–45 (2007).
Google Scholar 
Borenstein, M., Hedges, L. V., Higgins, J. P. & Rothstein, H. R. A basic introduction to fixed‐effect and random‐effects models for meta‐analysis. Res. Synth. Methods 1, 97–111 (2010).Article 

Google Scholar 
Harrer, M., Cuijpers, P., Furukawa, T. A. & Ebert, D. D. Doing Meta-Analysis with R: a Hands-on Guide. (Chapman and Hall, 2019).Wilke, C. O., Wickham, H. & Wilke, M. C. O. Package ‘cowplot’. Streamlined Plot Theme and Plot Annotations for ‘ggplot2 (Cowplot, 2019).Fill, J. M., Davis, C. N. & Crandall, R. M. Climate change lengthens southeastern USA lightning‐ignited fire seasons. Glob. Change Biol. 25, 3562–3569 (2019).Article 

Google Scholar 
Halofsky, J. E., Peterson, D. L. & Harvey, B. J. Changing wildfire, changing forests: the effects of climate change on fire regimes and vegetation in the Pacific Northwest, USA. Fire Ecol. 16, 4 (2020).Article 

Google Scholar 
Kraaij, T., Cowling, R. M., van Wilgen, B. W., Rikhotso, D. R. & Difford, M. Vegetation responses to season of fire in an aseasonal, fire-prone fynbos shrubland. PeerJ 5, e3591 (2017).Article 

Google Scholar 
Peel, M. C., Finlayson, B. L. & McMahon, T. A. Updated world map of the Köppen-Geiger climate classification. Hydrol. Earth Syst. Sci. 11, 1633–1644 (2007).Article 

Google Scholar 
Murphy, B. P. et al. Fire regimes of Australia: a pyrogeographic model system. J. Biogeogr. 40, 1048–1058 (2013).Article 

Google Scholar 
McColl-Gausden, S. C., Bennett, L. T., Duff, T. J., Cawson, J. G. & Penman, T. D. Climatic and edaphic gradients predict variation in wildland fuel hazard in south-eastern Australia. Ecography 43, 443–455 (2020).Article 

Google Scholar 
Pausas, J. G. & Keeley, J. E. Evolutionary ecology of resprouting and seeding in fire‐prone ecosystems. New Phytol. 204, 55–65 (2014).Article 

Google Scholar 
Lamont, B. B., Maitre, D. C. L., Cowling, R. M. & Enright, N. J. Canopy seed storage in woody plants. Bot. Rev. 57, 277–317 (1991).Article 

Google Scholar 
Tangney, R. et al. Data supporting: Success of post-fire plant recovery strategies varies with shifting fire seasonality. Zenodo https://doi.org/10.5061/dryad.7sqv9s4t5 (2022).Rothstein, H. R., Sutton, A. J. & Borenstein, M. Publication Bias in Meta-Analysis: Prevention, Assessment and Adjustments (John Wiley & Sons, 2006).Head, M. L., Holman, L., Lanfear, R., Kahn, A. T. & Jennions, M. D. The extent and consequences of P-hacking in science. PLoS Biol. 13, e1002106 (2015).Article 
CAS 

Google Scholar 
Egger, M., Smith, G. D., Schneider, M. & Minder, C. Bias in meta-analysis detected by a simple, graphical test. BMJ 315, 629–634 (1997).CAS 
Article 

Google Scholar 
Simonsohn, U., Nelson, L. D. & Simmons, J. P. P-curve: a key to the file-drawer. J. Exp. Psychol.Gen. 143, 534 (2014).Article 

Google Scholar  More

In the present work, we extend the base model for the spatial mosquito population dynamics24 to include wild male mosquitoes and genetically modified male mosquitoes. Thus, five populations will be considered: the aquatic mosquito population, including larvae and pupae, the egg mosquito population, the reproductive female mosquito population, the wild male mosquito population, and the genetically modified male population. Similar approaches can be found in the literature25,26.In the following system, we represent mosquito population densities (mosquitoes per m(^2)) by: E – in the egg phase, A – in the aquatic phase, F – female in the reproductive phase, M – wild males, and G – genetically modified male mosquitoes. Due to the very high resistance of the egg phase (up to 450 days27) and as we are interested in an urban spatial macro-scale modeling, we do not consider the mortality in the egg phase. The model is described by the following system of partial differential equations:$$begin{aligned} {left{ begin{array}{ll} partial _t E &{} = alpha beta F M -e E, \ partial _t A &{} = e left( 1 – dfrac{A}{k} right) E -(eta _a+{mu _a})A, \ partial _t F &{} = nabla cdot (D_m nabla F) -mu _f F + reta _{a} A, \ partial _t M &{} = nabla cdot (D_m nabla M) -mu _m M + (1-r)eta _{a} A, \ partial _t G &{} = nabla cdot (D_g nabla G) -mu _{g}G + l, end{array}right. } end{aligned}$$
(1)
where ( alpha ) represents the proportion of wild male mosquitoes to the total number of male mosquitoes (wild males + genetically modified males); (beta ) represents the expected quantity of eggs from the successful encounter between wild females and males; e is the egg hatching rate; k is the carrying capacity of the aquatic phase; ( eta _a ) is the emergence rate for mosquitoes from the aquatic phase to the female or male phases; ( mu _a), (mu _f), (mu _m), and (mu _{g}) are the mortality rates of mosquitoes in the aquatic phase, females, males, and genetically modified males, respectively; r is the proportion of females to males (typically (r=0.5)); (l=l(x,y,t)) is the function representing the number of genetically modified mosquitoes released in a unit of time at any point of the domain; (D_m) is the diffusion coefficient of wild mobiles females and males; (D_g) is the diffusion coefficient of genetically modified males. The proposed model (1) can naturally deal with heterogeneous parameters, such as mortality, diffusion, and carrying capacity coefficients. Thus it is possible to model the influence of rain, wind, and human action. In the context of this work, we are considering that the city neighborhood is divided into two environments: houses and streets. Due to lack of data, we restrict the investigated heterogeneity only to the carrying capacity coefficient.The proposed model can be regarded as an extension of other “economic” models20,24 in the effort to qualitatively reproduce the complex phenomena by using as few parameters as possible. Following this idea, the carrying capacity was neglected in the egg phase because of the skip oviposition phenomenon28 i.e., the female lays the number of eggs that the place holds, without more space, she migrates to other environments to finish laying the eggs. We also do not consider this coefficient in the winged phase as limitations in the winged phase were not reported in any study. On the other hand, we consider it in the aquatic phases (larvae and pupae), where it is effective29.The term ( alpha ), which multiplies the probability of encounters between male and female, represents the impact of the insertion of genetically modified males in the mosquito population to the immobile phase and is defined as$$begin{aligned} alpha = left{ begin{array}{cc} 1, &{} text{ if } M=G= 0, \ dfrac{M}{M + G}, &{} text{ otherwise }. end{array} right. end{aligned}$$
(2)
Similar modeling approach can be found in the literature16. As the release rate of genetically modified males increases, the alpha value decreases, and, consequently, the probability of encounter between females and wild males also decreases. Thus, there is a greater probability of encounter between genetically modified males and females. This approach presents an advantage, when compared to the models found in the literature25, as System  (1) does not present singularities at the equilibrium states, allowing mathematical analysis and numerical simulations. From the biological point of view, the increment of male wild mosquitoes over some critical value does not affect the egg deposition. At first glance, the term FM can lead to a misunderstanding that such property is not satisfied in the presented model. However, in Section “Equilibrium points considering the application of genetically modified male mosquitoes,” we argue that both male and female populations possess mathematical attractor equilibria, blocking the wild male population from growing beyond this value.Finally, any acceptable population model should be invariant in the definition domain, meaning its solution does not present senseless values. Setting the variable domain as$$begin{aligned} 0 le E(x,y,t)< infty ,;; 0 le A(x,y,t) le k, ;; 0 le F(x,y,t)< infty ,;; 0 le M(x,y,t)< infty ,;; 0 le G(x,y,t) < infty , end{aligned}$$ (3) we can verify that it is invariant under the time evolution by the System (1). To prove this statement, it is sufficient to verify that the vector field defined by the right side of (1) points into the domain when (E, A, F, M, G) approaches the domain boundary. When E approaches zero, the right side of the first equation in (1) is not negative. When A approaches zero, the right side of the second equation in (1) is not negative. When A approaches k (bottom), the first term on the right side of the second equation in (1) tends to zero, while the second term remains negative. Since the term ( nabla cdot (D_m nabla F) ) cannot change the F sign, when F approaches zero, the right side of the third equation in (1) is not negative Since the term ( nabla cdot (D_m nabla M) ) cannot change the M sign, when M approaches zero, the right side of the fourth equation in (1) is not negative. Since the term ( nabla cdot (D_g nabla G) ) cannot change the G sign, when G approaches zero, the right side of the fifth equation in (1) is not negative. In the rest of this section, let us explain how to estimate one-by-one all the parameters used in this model from experimental data available in the literature. It is a challenging task as, typically, the development of the Ae. aegypti mosquito depends on food variation30, temperature variations14,15 and rainfall31. This data is not available in the literature in the organized and systematic form. Because of that, we assume the environment is under optimal conditions of temperature, availability of food, and humidity.How to estimate emergence rate ((eta _a)) The emergence rate describes the rate at which the aquatic phase of the mosquito emerges into the adult phases. In the present model, for simplicity, it was considered that no mosquito from the crossing between genetically modified males and females reaches adulthood. Thus, the emergence rate is calculated on the crossing between females and wild males. Under optimal conditions and feeding distribution, based on the literature30, the emergence rate is 0.5596 (text{ day}^{-1}).How to estimate diffusion coefficients ((D_m,D_g)) The diffusion coefficient is one of the most important parameters describing the mosquitoes’ movement. We use the methodology proposed in the previous work24 to obtain the diffusion coefficient of adult mosquitoes (females and males) and genetically modified males.The estimate is done by assuming that all mosquitoes are released at (0, 0), and their movement is described by the corresponding equation in (1) neglecting other terms than diffusion. The population starts spreading in all directions. We define the spreading distance R(t) as the radius of the region centered in (0, 0) where (90%) of the initial mosquitoes population density is present. In Silva et al.24 it is shown that$$begin{aligned} R(t) = sqrt{4Dt} ;text {erf}^{-1}(0.9). end{aligned}$$ (4) Now corresponding diffusion coefficient is estimated by using the average flight distance of the mosquitoes and the characteristic time related to their life expectancy. Under favorable weather conditions, the average lifetime flight distance of females and males is approximately32,33 65 m, while the same for GM males is34 67.3 m. Based on the literature, we consider that the characteristic time for wild females and males32 is 7 days, and the same for genetically modified males is34 2.17 days. Using (4) we estimate the values for (D_m) and (D_g) summarized in Table 1. It would be natural to consider that the mosquitoes’ movement changes in different environments. Unfortunately, we were unable to find the corresponding experimental data, and because of that, we considered that (D_m) and (D_g) are the same in streets and house blocks.How to estimate mortality rates ((mu _a), (mu _f), (mu _m), (mu _{g}))The mortality coefficient represents an average quantity of mosquitoes in the corresponding phase dying each day. As mentioned before, we disregard the mortality rate in the egg phase, as it is negligible due to its great durability27, it does not affect the numerical results, and it complicates analytical estimates. Thus, the aquatic phase mortality rate coefficient is equal to the same for larvae’s coefficient, which is approximately29 (mu _a = 0.025) (1/day).There is no solid agreement on the mortality rate of male and female wild mosquitoes in the literature. Although some results29,30 suggest they are similar, we follow these authors and consider them equal. Considering both natural death and accidental ones, approximately (10%) of females and male mosquitoes in the adult phase die at each day35. Under optimal conditions, the mortality coefficient can be estimated from this data by using the proposed model (1) by neglecting diffusion and emergence terms in the corresponding equation; details can be found in the previous work24. The resulting parameter values are summarized in Table 1.It would be natural to consider that the mosquitoes mortality rate depends on the environment. Unfortunately, we were unable to find the corresponding experimental data, and because of that, we considered that (mu _a), (mu _f), (mu _m), and (mu _{g}) are the same in streets and house blocks.How to estimate the expected egg number ((beta ))This coefficient represents the average quantity of eggs a wild female lays per day, assuming a successful meeting with a wild male. Considering the number of times a female lays eggs in its lifetime36, the average quantity of eggs per lay and the mosquito’s life expectancy, under favorable conditions, this coefficient is estimated as (beta = 34).How to estimate the hatching rate (e)This coefficient determines the average number of eggs hatching in one day. Experimental data37 suggest that, under optimal humidity conditions, the mean value of the hatch rate coefficient is 0.24 given a temperature of 28 ((^{circ })C), which is considered ideal for mosquito development. This is the value used in the present work.How to estimate carrying capacity coefficient (k)The carrying capacity k represents the space limitation of one phase due to situations present in the environment37,38, such as competition for food among the larvae39. In general, it depends on external factors such as food availability, climate, terrain properties, making direct estimation almost impossible. In the Analytical results section, we show how to estimate this coefficient for each grid block. When considering spatial population dynamics in a heterogeneous environment, carrying capacity is one of the most influential parameters as it varies significantly. For example, house block offer more food and a shelter against natural predators resulting to a larger carrying capacity when compared with street environment. Following the literature32 we assume that the 80% of the mosquito’s breeding places are in houses resulting in the relation (k_h=5k_s), where (k_h) and (k_s) are the carrying capacities of the house blocks and in the streets.Genetically modified mosquitoes release rate (l)Function l(x, y, t) determines how many genetically modified mosquitoes are released in the location (x, y) at time t.In a normal situation, the sex ratio between males and females is 1 : 1. The increment of this proportion favoring GM males increases the probability of females to mate with these mosquitoes. As reported in the literature12,30 the initial launch size is 11 times larger than the adult female population, and it is done in some spots in the city. In this work, we analyze different release strategies maintaining the (11times 1) proportion in some scenarios.Table 1 All parameter values are directly taken or estimated from the literature as explained in section Modeling.Full size table More

DesignPlant richness. Sixteen locally frequent native plant species in the barren mountain areas (around Taizhou University, Zhejiang, China) invaded by the exotic plant Symphyotrichum subulatum60 were selected as the native species pool. These species were chosen because they spanned the dicotyledon plant taxonomy (including 7 Orders, 10 Families, and 14 Genus, in the Class Magnoliopsida), differed widely in their functional traits (related to height, life form, dominance in local communities, and leaf habit) (Supplementary Table 3), and were occasionally found to be associated with the invasive species Symphyotrichum subulatum60 in the local secondary-succession communities. With this species pool, we were able to imitate the locally natural, spatially stochastic, compositionally ruderal, and functionally varied plant community61, which is a typical attribute of the secondary-succession communities in the local barren mountains invaded by the exotic plant Symphyotrichum subulatum. Based on this native species pool, monocultures of each species (16 total), and random mixtures of 2, 4 or 8 species (with 10, 10, or 9 distinct assemblages, respectively) were designed, creating a complete set (Fig. 1d) of 45 different plant assemblages (pots) in total. Each plant assemblage was replicated 6 times, for a total of 270 pots. To eliminate the non-random effects during the 1-year development of the 270 pots, their distributions were randomized, such that not all replicates of an assemblage were next to each other (Fig. 1d–f).DroughtAfter 1-year development of the native plant assemblages, three drought treatments (non-, moderate-, and intensive-drought) were manipulated by adjusting irrigation using automatic drip irrigation systems, with 100%, 50%, and 25% of the equivalent to the amount received in the areas where native species were collected, respectively. Two random complete sets were selected for each drought treatment, each complete set being composed of 45 different plant assemblages (Fig. 1d–f).Exotic plant invasionNine months after drought treatment, the two complete sets (Fig. 1d) of each drought treatment were randomly exposed (invasion) or not exposed to (non-invasion) the invasive species Symphyotrichum subulatum (Michx.) G. L. Nesom (Fig. 1e, f). S. subulatum, an annual herbaceous plant native to North America, is a common invasive species in the subtropical and tropical regions of China18,60, and tends to interact with the native species via, for example, competing for space and resources62,63, enriching for pathogens or herbivores, and changing soil faunal, bacterial or fungal microbiomes18,64,65.ExperimentThe experiment based on the design mentioned above was conducted at Taizhou University, Zhejiang province, China (28.66°N, 121.39°E). The seeds of the 16 native plant species (Supplementary Table 3) and the soil were collected from nearby mountain areas (Wugui, 28.65°N, 121.38°E; Baiyun, 28.67°N, 121.42°E; Beigu, 28.86°N, 121.11°E). The seed-mixtures were obtained by mixing seeds of the 16 species pro rata, in proportion to germination rates. The soil (fine-loamy, mixed, semiative, mosic, Humic Hapludults) was sieved to pass a 2-mm mesh, and thoroughly mixed. 270 plastic pots (72 cm length × 64 cm width × 42 cm depth) were prepared, and each was filled with a 27-cm soil layer, followed by a 10-cm mixture of soil and vermiculite-compost to provide water-, air- and fertility-support for germination, seedling establishment, and plant growth (Supplementary Table 4).Native plant assemblagesAll the 270 pots were placed inside a plastic shelter, which allowed for both air ventilation and protection from rain. Each pot was sown with a seed-mixture of ca. 800 seeds. One month after germination, for each pot, the undesired seedlings were removed manually according to the plant richness design (Fig. 1d–f), and thus 32 vigorous seedlings (with the same number of seedlings per species, e.g., 4 seedlings for each species of the 8-species mixtures) were spatial-evenly retained. In this manner, the plant richness was manipulated for each plant assemblage. During the development of the 270 plant assemblages, the soil volumetric water content was controlled at ca. 20%, which was similar to that of the nearby mountainous soil, using the automatic drip irrigation systems. Weeds and undesired species were removed monthly (Fig. 1f).Drought treatmentAfter 1-year development of native plant assemblages, the drought treatments (non-, moderate-, and intensive-drought) were manipulated according to the experimental design mentioned above (Fig. 1d, e). Two complete sets (Fig. 1d) of different plant assemblages (2 × 45 pots) were selected for each drought treatment. Every other week, 40 pots each drought treatment were randomly selected for measuring soil water content and soil temperature at the depth of 0–20 cm, using the ProCheck analyzer (Decagon, Pullman, Washington, USA), and irrigation was adjusted accordingly using automatic drip irrigation systems. The irrigation for non-, moderate-, or intensive-drought was adjusted to accomplish an irrigation level amounts to 100%, 50%, or 25% that of the mountain areas where seeds were collected. Because of the distinct seasonal temperature and evaporation conditions, the irrigation frequencies were approximately daily in May-September, every other day in March–April and October–December, and weekly in January–February. With this manipulation, the volumetric soil water contents of non-, moderate-, and intensive-drought were controlled within ranges of 13.8–23.4%, 6.8–13.7%, and 1.4–7.4%, respectively, throughout the manipulation of drought treatment (Fig. 1e, f). Eight months after drought introduction, fresh litter was collected form the two replicate pots of each drought treatment, and then oven-dried at 40 °C, cut into ca. 2-cm pieces, and filled into litterbags (2-g litter in each litterbag).Invasion treatmentNine months after drought introduction, one complete set (45 pots) of the plant assemblages (Fig. 1d) from each drought treatment, was chosen and exposed to invasion disturbance by sowing 50 seeds of S. subulatum in each pot, and the other was specified as the non-invasion treatment (Fig. 1e, f). The prepared litterbags were embedded under the litter-layer of each pot (5 litterbags in each pot), correspondingly.SamplingSix months after invasion introduction, one litterbag was collected for litter-fauna extraction. Nine months after invasion, five soil cores (20-cm depth) were collected with augers (6.4 cm in diameter) and mixed for extraction of soil-fauna, and measurement of soil property and enzyme activity (Fig. 1f). The aboveground biomass of both native and invasive plants in each pot was harvested, sorted to species, oven-dried to a constant mass at 80 °C, and weighed. The belowground plant biomass was also sampled, sorted to native and invasive groups, oven-dried, and weighed (Fig. 1f).Plant, litter-, and soil-faunal communitiesPlant communitySince exotic plant invasion was treated as a disturbance factor, the biomass of the invasive species S. subulatum was not included for analyses concerning plant community and ecosystem (multi)functionality. The aboveground biomasses of native plant species in each of the 270 pots were collected for plant community analysis.Litter- and soil-faunal communitiesOne litterbag or fifty grams of mixed-soil samples were used for litter- or soil-fauna extraction using a Tullgren funnel apparatus (dry funnel method)66. The obtained microarthropods were stored in 70% alcohol, identified with double-tube anatomical lens, and classified to Family level. For both litter and soil samples, the numbers (abundances) of all faunal taxa were counted for litter/soil-faunal community analysis.Phylogenetic information of plant, litter-, and soil-faunal communitiesSimilar procedures were used to construct the plant and faunal phylogenetic trees. First, protein sequences of 12 faunal mitochondrial coding genes and 16 plant plastid coding genes (Supplementary Data 1) were obtained by searching plant or faunal taxonomies from NCBI protein database (https://www.ncbi.nlm.nih.gov/protein/) with Edirect software (https://www.ncbi.nlm.nih.gov/books/NBK179288/). All available sequences at plant species level or faunal Family level were fetched. If unavailable, the missing sequences were sampled from plant genus or faunal Order level. Sequoiadendron giganteum and Echinococcus were specified as out-group references for plant and faunal trees, respectively. Then, the sequences of each plant or faunal taxon were clustered at 97% or 90% identity independently, and the centroids were used as representative markers. The markers were aligned with MUSCLE67, followed by concatenation. Finally, using MEGA X68, the maximum likelihood trees were constructed based on BioNJ initial trees69 and 500 bootstrap checking nodal support. The parameters for plant tree construction were specified as follow: 70% partial deletion (with 4824 positions retained) and the best-fit substitution model JTT + G + I + F70,71; parameters for faunal tree: 90% partial deletion (2778 positions) and LG + G + I + F model71,72. The Linux codes for processing the protein sequences were submitted to GitHub (https://github.com/YuanGe-Lab/JZW_2022/tree/main/linux)The plant and faunal taxonomies, representative markers, and marker accessions are provided as Supplementary Data 1.Ecosystem function-related variablesA total of 14 individual function-related variables were collected. These variables belonged to three functional groups: (1) biomass production, including aboveground and belowground biomass of native plants, light interception efficiency, litter-fauna abundance, and soil-fauna abundance; (2) soil properties, including contents of soil organic carbon, soil nitrogen, soil phosphorus, and GRSP (relating to soil physical properties and stocks of carbon and nutrient73); and (3) processes, including rate of litter decomposition, and activities of β-glucosidase, protease, nitrate reductase and dehydrogenase.Light interception efficiency, the fraction of incident photosynthetically active radiation (PAR) intercepted by each plant community canopy, was determined between 12:00 and 14:00 on clear days using LI-191R line PAR sensors (LI-COR Inc., NE, USA), and the mean of 4 measurements (monthly from May to August the third year; Fig. 1f) was used. Total soil organic carbon and nitrogen were measured with an elemental analyzer (vario Max; Elementar, Germany). Total soil phosphorus was determined using the molybdenum blue method with a UV–visible spectrophotometer (Shimadzu, Kyoto, Japan). GRSP was determined using the method described by Shen et al.18. Litter decomposition rate was assessed by embedding litterbags and fitting litter mass loss against decomposition time (Fig. 1f). Enzyme activities were analyzed by the spectrophotometric method using the substrates, p-Nitrophenyl-β-d-glucopyranoside (pNPG; for β-glucosidase), caseinate (protease), nitrate (nitrate reductase) and triphenyltetrazolium chloride (TTC; dehydrogenase)18.Quantifying community stability and multifunctional stabilityCommunity data was comprised of native plant biomasses or faunal abundances, and the associated phylogenetic information. Multifunctionality data was comprised of 14 function-related variables, each variable (V) being transformed (V’) using the formula ({V}^{{prime} }=frac{V-{{{{{rm{min }}}}}}left(Vright)}{{{{{{rm{sd}}}}}}left(Vright)}) to guarantee even contribution to global variance. We calculated community similarity (1 minus Weighted-UniFrac distance) and multifunctional similarity (1 minus Bray–Curtis distance), based on the community data and the multifunctionality data, respectively. The specific subsets of each symmetric similarity matrix were used to assess three different aspects of stability: (1) Invariability (against stochastic fluctuations), reflected as the pairwise similarities (1476 pairs) within treatment groups, at same plant richness*drought*invasion condition; (2) Drought resistance, the similarities (2148 pairs) between drought (moderate- and intensive-drought) and non-drought treatments, at same plant richness*invasion condition; and (3) Invasion resistance, the similarities (n = 1611 pairs) between invasion and non-invasion treatments, at same plant richness*drought condition (Supplementary Fig. 1).We also assessed the three aspects of stability of each individual function in a similar way, but by calculating the similarity using the formula ({{{{{{{mathrm{SIM}}}}}}}}_{{ij}}=1-frac{|{V}_{i}-{V}_{j}|}{{V}_{i}+{V}_{j}}) (Vi and Vj are ith and jth elements in a function vector; SIMij is the similarity between Vi and Vj).Statistics and reproducibilityPERMANOVA (10,000 randomizations) was conducted to test the influences of the manipulated factors on ecosystem multifunctionality or communities of plant, litter- and soil-fauna, using “vegan::adonis” in R74. Mantel test (10,000 randomizations; Spearman’s R) was conducted to test the community-community or the community-multifunctionality relationships, using “vegan::mantel” in R74.As each similarity-pair of each aspect of community or multifunctional stability mentioned above was in strict correspondence to single level of each manipulated factor (plant richness, drought, and invasion) (Supplementary Fig.  1), the direct/indirect effects of treatments on the community or multifunctional stability can be assessed using SEM. To test direct and indirect effects (by modulating community stability) of the manipulated factors on multifunctional stability, we built three SEMs (Fig. 1a–c) based on three different aspects of stability (i.e., invariability, drought resistance, and invasion resistance) under the conditions of corresponding parings of manipulated factors (Supplementary Fig. 1), with the LAVAAN package75. The standardized paths (direct effects) in SEMs can be conceived as the partial correlations after teasing all side effects away. Bootstrapping with 10,000 randomizations was conducted to generate the unbiased mean effect. The significance of effect was tested using a Mantel-like permutation (10,000 randomizations) test76, where the null hypotheses (H0) were that the independent factors plant richness, drought, and invasion, had no direct/indirect effects (effect = 0) on multifunctional stability. Based on H0, permutation procedure was conducted by permuting the index of dependent factors (both columns and rows of a symmetric matrix; Supplementary Fig. 1) simultaneously to gain null models and null effects. p-values (probability of H0 acceptance) were calculated as the percentage of observed positive (or negative) effect that was greater (or less) than the null effects. We also assessed the direct and indirect effects of factors on the stability of each individual function based on the same SEMs, to consolidate our findings on multifunctional stability. The R codes and examples solving the permutation test for the significance of effects derived from SEMs that based on multidimensional similarity (or distance) were submitted to GitHub (https://github.com/YuanGe-Lab/JZW_2022/tree/main/R). All the analyses were conducted using R (https://www.r-project.org).Reporting summaryFurther information on research design is available in the Nature Research Reporting Summary linked to this article. More

Drainage area and branching ratio: a matter of scaleGeomorphological and ecological viewpoints on river networks generally differ owing to discordant definitions of the fundamental unit (the node) used to analyze them. From a geomorphological perspective, the determination of a river network entails the definition of an observational scale. Real river networks can be extracted from digital elevation models (DEMs) via algorithms for flow direction determination such as D8 (i.e., each pixel drains towards the lowest of its 8 nearest neighbors53). After the outlet location has been specified (and hence the upstream area A spanned by the river network), the first observational scale required is thus the pixel length l of the DEM, which defines the extent of a network node. A second scale is then needed to distinguish the portion of the drainage network effectively belonging to the channel network. The simplest but still widely used method53 defines channels as those pixels whose drainage area exceeds a threshold value AT. Hydrologically based criteria to determine the appropriate value for AT exist54; however, for the sake of simplicity, we here consider AT as a free parameter.BBTs and RBNs are random constructs, and as such they do not satisfy the optimality criterion of minimizing total energy expenditure, which is the fundamental physical process shaping fluvial landscapes. Furthermore, neither of these networks is a spanning tree, which is a key attribute of real fluvial landforms10: in fact, in both BBTs and RBNs, the extent of the drained domain is not defined. As a result, the drainage area at an arbitrary network node cannot in principle be attributed, unless by using the number of upstream nodes as a proxy. This has practical implications from an ecological viewpoint because drainage area is the master variable controlling several attributes of a river, such as width, depth, discharge, or slope3,55, which in turn impact habitat characteristics and the ecology of organisms therein56.In BBTs and RBNs, branching probability p has been defined35,38,45,46,47 as the probability that a network node is branching, i.e. connected to two upstream nodes. As such, the branching probability of a realized river network (be it a real river or a synthetic construct) could be evaluated as the ratio between the number of links NL constituting a network and the total number of network nodes N; if a unit distance between two adjacent nodes is assumed, the denominator equals the total network length. We note that the former definition of branching probability only holds in the context of the generation of a synthetic random network; it is in fact improper to refer to a “probability” when analyzing the properties of a realized river network. We clarify this aspect by introducing the concept of branching ratio pr for the latter definition (pr = NL/N). Moreover, in the case of BBTs, p and pr do not coincide (see Methods). Importantly, p and pr have no parallel in the literature on fluvial forms, nor do they refer to any of the well-studied measures of rivers’ fractal character.The choice of different observational scales for the same drainage network results in different values of NL and N, and hence of pr. Remarkably, the very same drainage network can result in river networks that virtually assume any value of pr (ranging from 0 to 1) and N (up to the upper bound A) depending on the choice of AT and A (the latter corresponding to a given l value when measured in the number of pixels; Fig. 1d–i); networks with low AT/A ratios result in high N (Fig. 2a), while networks with low AT result in high pr (Fig. 2b). Furthermore, pr does not identify the inherent (i.e., scale-independent) branching character of a given river network in relation to other river networks. In fact, by extracting different river networks at various scales (i.e., various AT values) and assessing the rivers’ rank in terms of pr, one observes that rivers that look more “branching” (i.e., have higher pr) than others for a given AT value can become less “branching” for a different AT value (Fig. 3). We therefore conclude that branching probability is a non-descriptive property of a river network, which by no means describes its inherent branching character, and depends on the observational scale.Fig. 2: Variation of N and pr as a function of observational scales for OCNs and real river networks.a Expected value of number of network nodes N as a function of threshold area AT and total drained area A (from Eq. (1)); the white dots indicate the values of AT and A used to generate the OCNs used in this analysis. b Expected value of branching ratio pr as a function of AT and A (from Eq. (1)); symbols as in a.Full size imageFig. 3: Values of branching ratio as a function of AT for the 50 real river networks analyzed in this study.a Natural values of pr in logarithmic scale. b z-normalized branching ratios (i.e., for each AT value, values of pr are normalized so that they have null mean and unit standard deviation), which better shows how rivers rank differently in terms of pr for different observation scales (i.e., AT). Lines connect dots relative to the same river. For visual purposes, rivers that rank first, second, second-to-last or last in at least one of the AT groups are displayed in colors; the other rivers are displayed in grey.Full size imageScaling is also crucial when looking at river networks from an ecological perspective. In this case, the relevant scale determining the dimension l of a node is the extent of habitat within which individuals (due to e.g. physical constrains) can be assigned to a single population57,58; the riverine connectivity and ensuing dispersal among these populations give rise to a metapopulation at the river network level. The specific spatial scale largely depends on the targeted species (e.g. being larger for fish than for aquatic insects), and it is conceivably much larger than (or, at least, it has no reason to be equal to) the pixel size of the DEM on which the river network is extracted. Since the evaluation of pr depends on the number of nodes N, which, in turn, is defined based on the scale length l, the resulting pr of a river network under this perspective would depend on the characteristics of the target taxa, which is inconsistent with the alleged role of pr as a scale-invariant property of river networks.Note also that using the ecological definition of l (i.e., spatial range of a local population) to discretize a real river network into N nodes, and from there calculate the branching ratio pr = NL/N, is problematic. Indeed, this would imply an elongation of all links shorter than l (which constitute a non-negligible fraction of the total links, under the assumption of exponential distribution of link lengths51), hence preventing a correct estimation of the connectivity patterns (i.e., distances between nodes) and the resulting ecological metrics of the river network (see section Ecological implications).From an ecological perspective, it could be reasonable to consider AT as a parameter expressing how a particular taxon perceives the suitable landscape, rather than a value to be determined from geomorphological arguments: for instance, large fishes inhabit wide and deep river reaches, and do not access small headwaters56. In this case, imposing a large AT would result in a coarser, less branching network constituted by few main channels (Fig. 1f, i), which could mimic the potentially available habitat for such species. Conversely, aquatic insects inhabit also small headwaters17,59, therefore their perceived landscape would resemble the finely resolved networks of Fig. 1d, g, characterized by low AT and higher (apparent) pr.Topology and scaling of river networks and random analoguesTo verify the topological (i.e., Horton’s laws on bifurcation and length ratios) and scaling (i.e., probability distribution of drainage areas) relationships of the different network types, we extracted from DEMs 50 real river networks encompassing a wide range of drainage areas (Fig. 4), and we generated 50 OCNs, 50 RBNs and 50 BBTs of comparable size (see Methods).Fig. 4: Location of real river basins used in the analysis.River basins are shown in dark grey; countries in light grey. Rivers’ numbering is sorted in ascending order according to drainage area values.Full size imageTypical values3,7,60 for the bifurcation ratio RB lie between 3 and 5, while length ratios (RL) range between 1.5 and 3.5. As expected, we observed that the real rivers and OCNs used in our analysis have RB and RL values within the aforementioned ranges (Fig. 5a, b). The same is true for RBNs, while the RB and RL values found for BBTs are lower than the typical ranges. This finding holds regardless of the scale (subsumed by AT) at which real river networks and OCNs are extracted (Supplementary Figs. 1 and 2). Remarkably, BBTs fail to satisfy Horton’s laws despite the statistical inevitability of such laws for any network argued by ref. 61. To this regard, we note that the networks analyzed by ref. 61 did not include constructs where all paths from the source nodes to the outlet have the same length, which is the defining feature of BBTs (Fig. 1a).Fig. 5: Comparison of topological and scaling properties of the different networks.a Scaling of number of network links Nω as a function of stream order ω for the various network types (rivers and OCNs obtained with AT = 20 pixels; RBNs and BBTs derived accordingly – see Methods). b Mean link length Lω (in units of l) as a function of ω. Networks used are as in panel a. c Scaling of drainage areas: probability P[A ≥ a] to randomly sample a node with drainage area A ≥ a as a function of a. The displayed trend lines are fitted on the ensemble values for the 50 network replicates, by excluding nodes with drainage area larger than 2000 pixels (cutoff value marked with a black solid line). The scaling coefficients β reported correspond to the slopes of the fitted trend lines. Extended details on all panels are provided in the Supplementary Methods.Full size imageWhile the power-law scaling of areas in OCNs (Fig. 5c) has an exponent β ≈ 0.45 that closely resembles the one found for the real rivers (β ≈ 0.46) and within the typically observed range8,10 β = 0.43 ± 0.02, drainage areas of RBNs scale as a power law with an exponent β ≈ 0.51, which departs from the observed range. Conversely, BBTs do not show any power-law scaling of areas. Scaling exponents of drainage areas fitted separately for each real river network yielded values in the range 0.36÷0.57 (Supplementary Table 1). In particular, we observed that these values tend to the expected range β = 0.43 ± 0.02 for increasing values of A, expressed in number of pixels (Supplementary Fig. 3), hence implying that highly resolved catchments are required in order to properly estimate β. Interestingly, the observed values of Horton ratios and scaling exponent β for RBNs are compatible with the values RB = 4, RL = 2, β = 0.5 predicted for Shreve’s random topology model3,60,62, which is actually equivalent to a RBN with infinite links.Ecological implicationsWe compared the different network types via two metrics that express the ecological value of a landscape for a metapopulation: the coefficient of variation of a metapopulation CVM and the metapopulation capacity λM. The coefficient of variation of a metapopulation63 is a measure of metapopulation stability (a metapopulation being more stable the lower CVM is), while the metapopulation capacity42,64 expresses the potential for a metapopulation to persist in the long run (persistence being more likely the higher λM is). Both measures are among the most universal metrics describing dynamics of spatially fragmented populations24,40. In order to assess the impact of the two landscape features mostly affecting metapopulation dynamics, i.e. spatial connectivity and spatial distribution of habitat patches, we calculated these metrics for the four network types under two different scenarios: uniform (CVM,U, λM,U) and non-uniform (CVM,H, λM,H) spatial distribution of habitat patch sizes. In the first scenario, CVM,U and λM,U assess stability and persistence (respectively) of a metapopulation solely based on pairwise distances between network nodes; in the second scenario, CVM,H and λM,H depend on the interplay between pairwise distances and spatially heterogeneous habitat availability (namely, downstream nodes being larger than upstream ones).We found that the values of CVM (be it derived with uniform (CVM,U) or nonuniform (CVM,H) distributions of patch sizes) obtained for OCNs match strikingly well those of real rivers (Fig. 6). These CVM values are consistently lower than those found for RBNs, while values of CVM for BBTs are even higher. Notably, this result holds for different values of AT (and hence different pr values) at which real rivers and OCNs are extracted (Fig. 6a–c; g–i), and for values of mean dispersal distance α (see Methods) spanning multiple orders of magnitude (Supplementary Figs. 4–7).Fig. 6: Comparison of values of metapopulation metrics across river network types and observational scales (AT).a–c CVM,U. d–f λM,U. g–i CVM,H. j–l λM,H. Boxplot elements are as follows: center line, median; notches, (pm 1.58cdot {{{{{{{rm{IQR}}}}}}}}/sqrt{50}), where IQR is the interquartile range; box limits, upper and lower quartiles; whiskers, extending up to the most extreme data points that are within ±1.5 ⋅ IQR; circles, outliers. Metapopulation metric values were obtained by setting α = 100 l. Note that in Eq. (1), given A = 40, 000, AT = 20 results in E[N] ≈ 4574, E[pr] ≈ 0.228; AT = 100 yields E[N] ≈ 2231, E[pr] ≈ 0.098; AT = 500 results in E[N] ≈ 1088, E[pr] ≈ 0.042.Full size imageFor a constant α value, the CVM of real rivers, OCNs and RBNs decreases as the resolution at which the network is extracted increases (i.e., AT decreases; see Fig. 6 and Supplementary Figs. 4–7). This is expected63, since a decrease in AT corresponds to an increase in N (Fig. 2a), leading to a decrease in CVM. Indeed, a larger ecosystem, constituted of more patches, has the potential to include a larger (and more diverse) number of subpopulations, which increases stability at a metapopulation level through statistical averaging–a phenomenon widely known as the portfolio effect65. We also found that BBT networks do not generally follow the above-described pattern of decreasing CVM with increasing N; rather, the CVM of BBTs increases with N when the mean dispersal distance α is set to intermediate to high values (Fig. 6 and Supplementary Figs. 5–7), and only when α is very low (e.g. α = 10 l as in Supplementary Fig. 4) and a uniform patch-size distribution is assumed does CVM,U follow the expected decreasing trend with increasing N.However, we need to warn against the conclusion that river networks with higher values of pr (and hence lower AT, see Fig. 2b) are inherently associated with higher metapopulation stability. Indeed, our result was obtained by changing the scale at which we observed the same river networks, and not by increasing the river networks’ size. If the number of network nodes (and, consequently, the branching ratio pr) is determined by the scale at which the landscape is observed, one cannot directly assume that any of such nodes is a node (or patch) in the ecological sense, i.e. the geographical span of a local population: the extent of such patches should be determined based on the mobility characteristics of the focus species, and should be independent of the scale at which the river network is observed. In contrast, we note that, if different river networks spanning different catchment areas (say, in km2) are compared, all of them extracted from the same DEM (same l and same AT in km2), then the larger river network will appear more branching (i.e., have larger pr). Indeed, by selecting catchments with larger A (in km2) for fixed l and AT (in km2), one moves towards the top-left corner of Fig. 2a, b (i.e., perpendicular to the level curves AT/A). The apparent higher “branchiness” of the river network with larger A will result in lower values of CVM; however, the higher metapopulation stability of the larger network will not be due to its (alleged) inherent more branching character, but only dictated by its larger habitat availability.We observed that metapopulation capacity λM values of OCNs (be it evaluated under uniform (λM,U) or non-uniform (λM,H) patch-size distribution assumption) are the closest to those of real rivers, while RBNs (and even more so BBTs) generally overestimate λM with respect to real rivers and OCNs (Fig. 6d–f; j–l). This result holds irrespective of the choice of AT and for intermediate to high values of α (Supplementary Figs. 5–7). When the mean dispersal distance is instead set to very low values (α = 10 l – Supplementary Fig. 4) and the river network is extracted at a high resolution (i.e., low AT), the metapopulation capacity of OCNs under assumption of uniform patch-size distribution (λM,U) is underestimated with respect to that of real rivers. A likely explanation for this apparent mismatch is that, for low values of AT, the number of nodes N tends to be somewhat higher for the extracted river networks used in this analysis than for OCNs (Supplementary Fig. 8), and the effect of the different dimensionality of real rivers and OCNs in the metapopulation capacity estimation tends to be more evident as the mean dispersal distance decreases. Interestingly, such mismatch is absent when a non-uniform patch size distribution is assumed, as λM,H values for OCNs match those for real rivers regardless of the mean dispersal distance value and the river network resolution (Fig. 6; Supplementary Figs. 4–7).The OCN construct encapsulates both random and deterministic processes, the former related to the stochastic nature of the OCN generation algorithm, and the latter pertaining to the minimization of total energy expenditure that characterizes OCN configurations. As such, OCNs reproduce the aggregation patterns of real river networks. From an ecological viewpoint, this implies that both pairwise distances between nodes and the distribution of patch sizes (expressed as a function of drainage areas, or of a proxy thereof such as the number of nodes upstream) are much closer to those of real networks than is the case for fully random synthetic networks as BBTs and RBNs. In particular, BBTs and (to a lesser extent) RBNs tend to underestimate pairwise distances with respect to real rivers and OCNs, as documented by a comparison of mean pairwise distances across network types (Supplementary Fig. 9a–c). Our analysis shows that the connectivity structure of these random networks (subsumed by the matrix of pairwise distances) is too compact with respect to that of real rivers, which leads to an overestimation of the role of dispersal in increasing the ability of a metapopulation to persist in the long run, but also an increased likelihood of synchrony among the different local populations, which results in higher instability.Comparison of patch size distributions among the network types expressed in terms of CVM,0 (i.e., the portion of CVM,H that uniquely depends on the distribution of patch sizes and not on pairwise distances) shows that, while for coarsely resolved networks (AT = 500) no clear differences in CVM,0 emerged, for highly resolved networks (AT = 20) BBTs heavily underestimate the CVM,0 of real rivers and OCNs, while RBNs slightly overestimate it (Supplementary Fig. 9d–f). As a result of the interplay of differences in distance matrices and patch size distributions, BBTs and (to a lesser extent) RBNs generally tend to overestimate the coefficient of variation of a metapopulation and the metapopulation capacity of real rivers and OCNs in both scenarios of uniform and non-uniform patch size distribution. The only exception to this trend occurs for the metapopulation capacity λM,H of very large BBTs (corresponding to AT = 20) in the case of very high dispersal distances (α = 1000 l – Supplementary Fig. 7): here, the patch-size effect (i.e., underestimation of CVM,0) predominates over the distance effect (i.e., overestimation of mean dij), resulting in an underestimation of λM,H with respect to real rivers and OCNs.Our results were derived under a number of simplifying assumptions. In particular, we acknowledge that, while the distance matrix of a landscape and the distribution of patch sizes have in general important implications for metapopulation dynamics, other factors not considered here, such as Euclidean between-patch distance48, fat-tailed dispersal kernel66 and density-dependent dispersal67 could also play a relevant role in this respect. However, it needs to be noted that, especially with regards to the assessment of the Moran effect in metapopulation synchrony (i.e., increased synchrony in local fluvial populations that are geographically close but not flow-connected48), the use of OCNs allows integration of Euclidean distances in a metapopulation model, while this is not possible for RBNs and BBTs, where Euclidean distances are not defined. Moreover, if a larger degree of realism is required for a specific ecological modelling study, such as heterogeneity in abiotic factors (e.g. water temperature or flow rates), the use of OCNs as model landscapes allows a direct integration of these variables, as they can conveniently be expressed as functions of drainage area3,55. In contrast, this is not possible for RBNs or BBTs, because only OCNs verify the scaling of areas (Fig. 5c), while RBNs and BBTs lack a proper definition of drainage areas.Our comparison of synthetic and real river networks showed that riverine metapopulations are more stable and less invasible than what would be predicted by random network analogues. Conversely, the use of OCNs as model landscapes allows capturing not only the scaling features of real rivers, but also drawing ecological conclusions that are in line with those that could be observed in real river networks. We thus support the use of OCNs as analogues of real river networks in theoretical and applied ecological modelling studies. While we found that BBTs are highly inaccurate in reproducing ecological metrics of real river networks and should be therefore discarded altogether in future modelling applications, RBNs show a certain degree of similarity with OCNs and real river networks in this respect; moreover, RBNs (as is the case for any random tree61) satisfy Horton’s laws on bifurcation and length ratios. A relevant advantage of RBNs over OCNs is that their generation algorithm is at least one order of magnitude faster49. Therefore, we acknowledge that RBNs could be considered as a suitable surrogate for real river networks as null models in cases where a large number of network replicates is required. To this end, we encourage researchers exploiting synthetic river networks (whether they be OCNs or RBNs) to always clarify the observational scales (that is, total area drained, size of a node, area drained by a headwater) subsumed by the synthetic network and which give rise to a certain complexity measure (i.e., branching ratio). Only in such a way could the predictions from these studies be compared with real river networks.In conclusion, our results advocate a tighter integration between physical (geomorphology, hydrology) and biological (ecology) disciplines in the study of freshwater ecosystems, and particularly in the perspective of a mechanistic understanding of drivers of persistence and loss of biodiversity. More

Harvey, P. H. & Pagel, M. D. The comparative method in evolutionary biology. Vol. 239 (Oxford University Press, 1991).Oyston, J. W., Hughes, M., Wagner, P. J., Gerber, S. & Wills, M. A. What limits the morphological disparity of clades? Interface Focus 5, 0042 (2015).Article 

Google Scholar 
Jetz, W., Thomas, G. H., Joy, J. B., Hartmann, K. & Mooers, A. O. The global diversity of birds in space and time. Nature 491, 444–448 (2012).CAS 
PubMed 
Article 

Google Scholar 
Webb, C. O. Exploring the phylogenetic structure of ecological communities: an example for rain forest trees. Am. Naturalist 156, 145–155 (2000).Article 

Google Scholar 
Purvis, A., Gittleman, J. L. & Brooks, T. Phylogeny and conservation. (Cambridge University Press, 2005).Page, R. D. M. Parallel phylogenies: reconstructing the history of host-parasite assemblages. Cladistics 10, 155–173 (1994).Article 

Google Scholar 
Weaver, S. C. & Vasilakis, N. Molecular evolution of dengue viruses: contributions of phylogenetics to understanding the history and epidemiology of the preeminent arboviral disease. Infect., Genet. Evolution 9, 523–540 (2009).CAS 
Article 

Google Scholar 
Tassy, P. Trees before and after Darwin. J. Zool. Syst. Evolut. Res. 49, 89–101 (2011).Article 

Google Scholar 
Heather, J. M. & Chain, B. The sequence of sequencers: The history of sequencing DNA. Genomics 107, 1–8 (2016).CAS 
PubMed 
Article 

Google Scholar 
Pyron, R. A. Post-molecular systematics and the future of phylogenetics. Trends Ecol. Evolution 30, 384–389 (2015).Article 

Google Scholar 
Sansom, R. S. & Wills, M. A. Differences between hard and soft phylogenetic data. Proc. R. Soc. B: Biol. Sci. 284, 20172150 (2017).Article 

Google Scholar 
Scotland, R. W., Olmstead, R. G. & Bennett, J. R. Phylogeny reconstruction: the role of morphology. Syst. Biol. 52, 539–548 (2003).PubMed 
Article 

Google Scholar 
Regier, J. C. et al. Arthropod relationships revealed by phylogenomic analysis of nuclear protein-coding sequences. Nature 463, 1079–1083 (2010).CAS 
PubMed 
Article 

Google Scholar 
Callender-Crowe, L. M. & Sansom, R. S. Osteological characters of birds and reptiles are more congruent with molecular phylogenies than soft characters are. Zool. J. Linn. Soc. 194, 1–13 (2022).Article 

Google Scholar 
Wahlberg, N. et al. Synergistic effects of combining morphological and molecular data in resolving the phylogeny of butterflies and skippers. Proc. R. Soc. B: Biol. Sci. 272, 1577–1586 (2005).CAS 
Article 

Google Scholar 
He, L. et al. A molecular phylogeny of selligueoid ferns (Polypodiaceae): Implications for a natural delimitation despite homoplasy and rapid radiation. Taxon 67, 237–249 (2018).Article 

Google Scholar 
Fernández, R., Edgecombe, G. D. & Giribet, G. Phylogenomics illuminates the backbone of the Myriapoda Tree of Life and reconciles morphological and molecular phylogenies. Sci. Rep. 8, 1–7 (2018).
Google Scholar 
Eme, L., Spang, A., Lombard, J., Stairs, C. W. & Ettema, T. J. G. Archaea and the origin of eukaryotes. Nat. Rev. Microbiol. 15, 711–723 (2017).CAS 
PubMed 
Article 

Google Scholar 
Asher, R. J., Bennett, N. & Lehmann, T. The new framework for understanding placental mammal evolution. BioEssays 31, 853–864 (2009).CAS 
PubMed 
Article 

Google Scholar 
Shoshani, J. & McKenna, M. C. Higher taxonomic relationships among extant mammals based on morphology, with selected comparisons of results from molecular data. Mol. Phylogenetics Evolution 9, 572–584 (1998).CAS 
Article 

Google Scholar 
Beck, R. M. D. & Baillie, C. Improvements in the fossil record may largely resolve current conflicts between morphological and molecular estimates of mammal phylogeny. Proc. R. Soc. B: Biol. Sci. 285, 20181632 (2018).Article 

Google Scholar 
Zou, Z. T. & Zhang, J. Z. Morphological and molecular convergences in mammalian phylogenetics. Nat. Commun. 7, 1–9 (2016).
Google Scholar 
Hillis, D. M. Molecular versus morphological approaches to systematics. Annu. Rev. Ecol. Syst. 18, 23–42 (1987).Article 

Google Scholar 
Thompson, N. Alfred Russell Wallace Contributions to the theory of Natural Selection, 1870, and Charles Darwin and Alfred Wallace, ‘On the Tendency of Species to form Varieties’ (Papers presented to the Linnean Society 30th June 1858). (Routledge, 2004).Croizat, L. Panbiogeography; or an introductory synthesis of zoogeography, phytogeography, and geology, with notes on evolution, systematics, ecology, anthropology, etc., Vol. 1, 2a & 2b (Published by the author, Caracas., 1958).Means, J. C. & Marek, P. E. Is geography an accurate predictor of evolutionary history in the millipede family Xystodesmidae? PeerJ 5, e3854 (2017).PubMed 
PubMed Central 
Article 

Google Scholar 
Wills, M. A., Barrett, P. M. & Heathcote, J. F. The modified gap excess ratio (GER*) and the stratigraphic congruence of dinosaur phylogenies. Syst. Biol. 57, 891–904 (2008).PubMed 
Article 

Google Scholar 
Fisher, D. C. Stratocladistics: integrating temporal data and character data in phylogenetic inference. Annu. Rev. Ecol., Evolution Syst. 39, 365–385 (2008).Article 

Google Scholar 
Lazarus, D. B. & Prothero, D. R. The role of stratigraphic and morphologic data in phylogeny. J. Paleontol. 58, 163–172 (1984).
Google Scholar 
Camerini, J. R. Evolution, biogeography, and maps: an early history of Wallace’s Line. Isis 84, 700–727 (1993).CAS 
PubMed 
Article 

Google Scholar 
Upchurch, P., Hunn, C. A. & Norman, D. B. An analysis of dinosaurian biogeography: evidence for the existence of vicariance and dispersal patterns caused by geological events. Proc. R. Soc. B: Biol. Sci. 269, 613–621 (2002).Article 

Google Scholar 
Ferreira, G. S., Bronzati, M., Langer, M. C. & Sterli, J. Phylogeny, biogeography and diversification patterns of side-necked turtles (Testudines: Pleurodira). R. Soc. Open Sci. 5, 171773 (2018).PubMed 
PubMed Central 
Article 

Google Scholar 
Ronquist, F. & Sanmartín, I. Phylogenetic methods in biogeography. Annu. Rev. Ecol., Evolution, Syst. 42, 441–464 (2011).Article 

Google Scholar 
IUCN. The IUCN Red List of Threatened Species. Version 2019-2., https://www.iucnredlist.org (2019).GBIF.org. GBIF Home Page, https://www.gbif.org/ (2019).Uetz, P., Freed, P., Aguilar, R. & Hošek, J. The reptile database., http://www.reptiledatabase.org (2019).Archie, J. W. Homoplasy excess ratios: new indices for measuring levels of homoplasy in phylogenetic systematics and a critique of the consistency index. Syst. Zool. 38, 253–269 (1989).Article 

Google Scholar 
Wilkinson, M. On phylogenetic relationships within Dendrotriton (Amphibia: Caudata: Plethodontidae) is there sufficient evidence? Herpetological J. 7, 55–65 (1997).
Google Scholar 
O’Connor, A. & Wills, M. A. Measuring stratigraphic congruence across trees, higher taxa, and time. Syst. Biol. 65, 792–811 (2016).PubMed 
PubMed Central 
Article 

Google Scholar 
Colless, D. H. Review of phylogenetics: the theory and practice of phylogenetic systematics. Syst. Zool. 31, 100–104 (1982).Article 

Google Scholar 
Lartillot, N. & Philippe, H. Improvement of molecular phylogenetic inference and the phylogeny of Bilateria. Philos. Trans. R. Soc. B: Biol. Sci. 363, 1463–1472 (2008).Article 

Google Scholar 
Sansom, R. S., Choate, P. G., Keating, J. N. & Randle, E. Parsimony, not Bayesian analysis, recovers more stratigraphically congruent phylogenetic trees. Biol. Lett. 14, 20180263 (2018).PubMed 
PubMed Central 
Article 

Google Scholar 
Rosa, B. B., Melo, G. A. & Barbeitos, M. S. Homoplasy-based partitioning outperforms alternatives in Bayesian analysis of discrete morphological data. Syst. Biol. 68, 657–671 (2019).PubMed 
Article 

Google Scholar 
Lucena, D. A. & Almeida, E. A. Morphology and Bayesian tip-dating recover deep Cretaceous-age divergences among major chrysidid lineages (Hymenoptera: Chrysididae). Zool. J. Linn. Soc. 194, 36–79 (2022).Article 

Google Scholar 
O’Reilly, J. E. et al. Bayesian methods outperform parsimony but at the expense of precision in the estimation of phylogeny from discrete morphological data. Biol. Lett. 12, 20160081 (2016).PubMed 
PubMed Central 
Article 

Google Scholar 
Smith, M. R. Bayesian and parsimony approaches reconstruct informative trees from simulated morphological datasets. Biol. Lett. 15, 20180632 (2019).PubMed 
PubMed Central 
Article 

Google Scholar 
Wiens, J. The role of morphological data in phylogeny reconstruction. Syst. Biol. 53, 653–661 (2004).PubMed 
Article 

Google Scholar 
O’Leary, M. A. & Kaufman, S. G. MorphoBank 3.0: Web application for morphological phylogenetics and taxonomy., http://www.morphobank.org (2012).de Queiroz, A. & Gatesy, J. The supermatrix approach to systematics. Trends Ecol. Evolution 22, 34–41 (2007).Article 

Google Scholar 
Wilkinson, M. A comparison of two methods of character construction. Cladistics 11, 297–308 (1995).Article 

Google Scholar 
Brazeau, M. D. Problematic character coding methods in morphology and their effects. Biol. J. Linn. Soc. 104, 489–498 (2011).Article 

Google Scholar 
Drummond, A. J., Ho, S. Y. W., Phillips, M. J. & Rambaut, A. Relaxed phylogenetics and dating with confidence. PLoS Biol. 4, e88 (2006).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
O’Reilly, J. E., Puttick, M. N., Pisani, D. & Donoghue, P. C. Probabilistic methods surpass parsimony when assessing clade support in phylogenetic analyses of discrete morphological data. Palaeontology 61, 105–118 (2018).PubMed 
Article 

Google Scholar 
Keating, J. N., Sansom, R. S., Sutton, M. D., Knight, C. G. & Garwood, R. J. Morphological phylogenetics evaluated using novel evolutionary simulations. Syst. Biol. 69, 897–912 (2020).PubMed 
PubMed Central 
Article 

Google Scholar 
Makarenkov, V. et al. Weighted bootstrapping: a correction method for assessing the robustness of phylogenetic trees. BMC Evolut. Biol. 10, 1–16 (2010).Article 
CAS 

Google Scholar 
Stayton, C. T. The definition, recognition, and interpretation of convergent evolution, and two new measures for quantifying and assessing the significance of convergence. Evolution 69, 2140–2153 (2015).PubMed 
Article 

Google Scholar 
Sattler, R. Homology – a continuing challenge. Syst. Bot. 9, 382–394 (1984).Article 

Google Scholar 
Jenner, R. A. & Schram, F. R. The grand game of metazoan phylogeny: rules and strategies. Biol. Rev. 74, 121–142 (1999).Article 

Google Scholar 
Pisani, D. & Wilkinson, M. Matrix representation with parsimony, taxonomic congruence, and total evidence. Syst. Biol. 51, 151–155 (2002).PubMed 
Article 

Google Scholar 
Arcila, D. et al. Testing the utility of alternative metrics of branch support to address the ancient evolutionary radiation of tunas, stromateoids, and allies (Teleostei: Pelagiaria). Syst. Biol. 70, 1123–1144 (2021).PubMed 
Article 

Google Scholar 
Felsenstein, J. Phylogenies and the comparative method. Am. Naturalist 125, 1–15 (1985).Article 

Google Scholar 
Bremer, K. Branch support and tree stability. Cladistics 10, 295–304 (1994).Article 

Google Scholar 
Johnson, W. E. et al. The late Miocene radiation of modern Felidae: a genetic assessment. Science 311, 73–77 (2006).CAS 
PubMed 
Article 

Google Scholar 
Van der Made, J. Biogeography and climatic change as a context to human dispersal out of Africa and within Eurasia. Quat. Sci. Rev. 30, 1353–1367 (2011).Article 

Google Scholar 
May, F., Rosenbaum, B., Schurr, F. M. & Chase, J. M. The geometry of habitat fragmentation: Effects of species distribution patterns on extinction risk due to habitat conversion. Ecol. Evolution 9, 2775–2790 (2019).Article 

Google Scholar 
Swofford, D. L. et al. Bias in phylogenetic estimation and its relevance to the choice between parsimony and likelihood methods. Syst. Biol. 50, 525–539 (2001).CAS 
PubMed 
Article 

Google Scholar 
Jaeger, J. J. & Martin, M. African marsupials – vicariance or dispersion? Nature 312, 379–379 (1984).Article 

Google Scholar 
Smith, B. T. et al. The drivers of tropical speciation. Nature 515, 406–409 (2014).CAS 
PubMed 
Article 

Google Scholar 
Simkanin, C. et al. Exploring potential establishment of marine rafting species after transoceanic long-distance dispersal. Glob. Ecol. Biogeogr. 28, 588–600 (2019).Article 

Google Scholar 
Raxworthy, C. J., Forstner, M. R. J. & Nussbaum, R. A. Chameleon radiation by oceanic dispersal. Nature 415, 784–787 (2002).CAS 
PubMed 
Article 

Google Scholar 
Stehli, F. G. & Webb, S. D. The great American biotic interchange., Vol. 4 (Springer Science & Business Media, 2013).Ronquist, F. Dispersal-vicariance analysis: A new approach to the quantification of historical biogeography. Syst. Biol. 46, 195–203 (1997).Article 

Google Scholar 
Ricklefs, R. E. & Bermingham, E. The concept of the taxon cycle in biogeography. Glob. Ecol. Biogeogr. 11, 353–361 (2002).Article 

Google Scholar 
Ma, H. An analysis of the equilibrium of migration models for biogeography-based optimization. Inf. Sci. 180, 3444–3464 (2010).Article 

Google Scholar 
Yiming, L., Niemelä, J. & Dianmo, L. Nested distribution of amphibians in the Zhoushan archipelago, China: can selective extinction cause nested subsets of species? Oecologia 113, 557–564 (1998).CAS 
PubMed 
Article 

Google Scholar 
Crisci, J. V., Katinas, L. & Posadas, P. Historical Biogeography: An Introduction. (Harvard University Press, 2003).Chen, R. et al. Adaptive innovation of green plants by horizontal gene transfer. Biotechnol. Adv. 46, 107671 (2021).CAS 
PubMed 
Article 

Google Scholar 
Schönknecht, G., Weber, A. P. & Lercher, M. J. Horizontal gene acquisitions by eukaryotes as drivers of adaptive evolution. BioEssays 36, 9–20 (2014).PubMed 
Article 
CAS 

Google Scholar 
Smith, A. B. Echinoderm phylogeny: morphology and molecules approach accord. Trends Ecol. Evolution 7, 224–229 (1992).CAS 
Article 

Google Scholar 
Bateman, R. M., Hilton, J. & Rudall, P. J. Morphological and molecular phylogenetic context of the angiosperms: contrasting the ‘top-down’ and ‘bottom-up’ approaches used to infer the likely characteristics of the first flowers. J. Exp. Bot. 57, 3471–3503 (2006).CAS 
PubMed 
Article 

Google Scholar 
Morris, J. L. et al. The timescale of early land plant evolution. Proc. Natl Acad. Sci. 115, E2274–E2283 (2018).CAS 
PubMed 
PubMed Central 

Google Scholar 
Richter, S. The Tetraconata concept: hexapod-crustacean relationships and the phylogeny of Crustacea. Org. Diversity Evolution 2, 217–237 (2002).Article 

Google Scholar 
Dunn, C. W. et al. Broad phylogenomic sampling improves resolution of the animal tree of life. Nature 452, 745–749 (2008).CAS 
PubMed 
Article 

Google Scholar 
Caravas, J. & Friedrich, M. Of mites and millipedes: recent progress in resolving the base of the arthropod tree. BioEssays 32, 488–495 (2010).CAS 
PubMed 
Article 

Google Scholar 
Howard, R. J. et al. The Ediacaran origin of Ecdysozoa: integrating fossil and phylogenomic data. J. Geol. Soc. https://doi.org/10.1144/jgs2021-107 (2022).Newman, M. E. J. A model of mass extinction. J. Theor. Biol. 189, 235–252 (1997).CAS 
PubMed 
Article 

Google Scholar 
Cobbett, A., Wilkinson, M. & Wills, M. A. Fossils impact as hard as living taxa in parsimony analyses of morphology. Syst. Biol. 56, 753–766 (2007).PubMed 
Article 

Google Scholar 
Ruta, M., Krieger, J., Angielczyk, K. & Wills, M. A. The evolution of the tetrapod humerus: morphometrics, disparity, and evolutionary rates. Earth Environ. Sci. Trans. R. Soc. Edinb. 109, 351–369 (2018).
Google Scholar 
Puttick, M. N., Thomas, G. H. & Benton, M. J. High rates of evolution preceded the origins of birds. Evolution 68, 1497–1510 (2014).PubMed 
PubMed Central 
Article 

Google Scholar 
Sansom, R. S. & Wills, M. A. Fossilization causes organisms to appear erroneously primitive by distorting evolutionary trees. Sci. Rep. 3, 1–5 (2013).Article 

Google Scholar 
Brinkworth, A., Sansom, R. & Wills, M. A. Phylogenetic incongruence and homoplasy in the appendages and bodies of arthropods: why broad character sampling is best. Zool. J. Linn. Soc. 187, 100–116 (2019).Article 

Google Scholar 
Brown, J. W. & Smith, S. A. The past sure is tense: on interpreting phylogenetic divergence time estimates. Syst. Biol. 67, 340–353 (2018).PubMed 
Article 

Google Scholar 
Barba-Montoya, J., Dos Reis, M. & Yang, Z. H. Comparison of different strategies for using fossil calibrations to generate the time prior in Bayesian molecular clock dating. Mol. Phylogenetics Evolution 114, 386–400 (2017).CAS 
Article 

Google Scholar 
Sanderson, M. J. & Donoghue, M. J. Patterns of variation in levels of homoplasy. Evolution 43, 1781–1795 (1989).PubMed 
Article 

Google Scholar 
Alroy, J. Fossilworks: Gateway to the Paleobiology Database, http://fossilworks.org (2019).Benton, M. J. The Fossil Record 2. (Chapman & Hall, 1993).Cohen, K. M., Harper, D. A. T. & Gibbard, P. L. ICS International Chronostratigraphic Chart 2021/02, http://www.stratigraphy.org/ (2021).Gradstein, F. & Ogg, J. Geologic time scale 2004–why, how, and where next! Lethaia 37, 175–181 (2004).Article 

Google Scholar 
Rohde, R. A. The GeoWhen Database, (2005).O’Leary, M. A. et al. The placental mammal ancestor and the post–K-Pg radiation of placentals. Science 339, 662–667 (2013).PubMed 
Article 
CAS 

Google Scholar 
Kluge, A. G. A concern for evidence and a phylogenetic hypothesis of relationships among Epicrates (Boidae, Serpentes). Syst. Biol. 38, 7–25 (1989).Article 

Google Scholar 
Tolson, P. J. Phylogenetics of the boid snake genus Epicrates and Caribbean vicariance theory. Occasional Pap. Mus. Zool., Univ. Mich. 715, 1–68 (1987).
Google Scholar 
Clopper, C. J. & Pearson, E. S. The use of confidence or fiducial limits illustrated in the case of the binomial. Biometrika 26, 404–413 (1934).Article 

Google Scholar  More

Costanza, R. et al. Nature 387, 253–260 (1997).Article 

Google Scholar 
Toman, M. Ecol. Econom. 25, 57–60 (1998).Article 

Google Scholar 
Barbier, E. B., Burgess, J. C. & Folke, C. (eds) Paradise Lost? The Ecological Economics of Biodiversity (Routledge, 1994).
Google Scholar 
Daily, G. C. (ed.) Nature’s Services: Societal Dependence on Natural Ecosystems (Island, 1997).
Google Scholar 
Watson, R. T. et al. (eds) Global Biodiversity Assessment: Summary for Policy-Makers (Cambridge Univ. Press, 1995).
Google Scholar 
Reid, W. V. et al. Ecosystems and Human Well-Being: Synthesis (A Report of the Millennium Ecosystem Assessment) (Island, 2005).
Google Scholar 
Chichilnisky, G. & Heal, G. Nature 391, 629–630 (1998).Article 

Google Scholar 
Umaña Quesada, A. in Green Growth That Works: Natural Capital Policy and Finance Mechanisms from Around the World (eds Mandle, L. et al.) Ch. 13, 195–212 (Island, 2019).
Google Scholar 
Ouyang, Z. et al. in Green Growth That Works: Natural Capital Policy and Finance Mechanisms from Around the World (eds Mandle, L. et al.) Ch. 12, 177–194 (Island, 2019).
Google Scholar 
Salzman, J., Bennett, G., Carroll, N., Goldstein, A. & Jenkins, M. Nature Sustain. 1, 136–144 (2018).Article 

Google Scholar 
Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services. The Global Assessment Report on Biodiversity and Ecosystem Services: Summary for Policymakers (eds Díaz, S. et al.) (IPBES, 2019).
Google Scholar 
Chaplin-Kramer, R. et al. Science 366, 255–258 (2019).PubMed 
Article 

Google Scholar 
Hamel, P. et al. npj Nature Urban Sustain. 1, 25 (2021).Article 

Google Scholar 
Mandle, L., Ouyang, Z., Salzman, J. & Daily, G. C. (eds) Green Growth That Works: Natural Capital Policy and Finance Mechanisms from Around the World (Island, 2019).
Google Scholar 
Dasgupta, P. The Economics of Biodiversity: The Dasgupta Review. Abridged Version (HM Treasury, 2021).
Google Scholar 
Arkema, K. K. et al. Proc. Natl Acad. Sci. USA 112, 7390–7395 (2015).PubMed 
Article 

Google Scholar 
Ouyang, Z. et al. Proc. Natl Acad. Sci. USA 117, 14593–14601 (2020).PubMed 
Article 

Google Scholar  More