Ecology

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

Woese, C. R., Kandler, O. & Wheelis, M. L. Proc. Natl Acad. Sci. USA 87, 4576–4579 (1990).CAS 
Article 

Google Scholar 
Sapp, J. Microbiol. Mol. Biol. Rev. 69, 292–305 (2005).CAS 
Article 

Google Scholar 
Chatton, É. Ann. Sci. Nat. Zool. 8, 1–84 (1925).
Google Scholar 
Soyer-Gobillard, M.-O. & Schrevel, J. The Discoveries and Artistic Talents of Édouard Chatton and André Lwoff, Famous Biologists (Cambridge Scholars Publishing, 2020).Macfarlane, J. M. The Causes and Course of Organic Evolution: A Study in Bioenergics (Macmillan, 1918).Haeckel, E. Systematische Phylogenie. Erster Theil (Verlag von Georg Reimer, 1894).Stanier, R. Y., Douderoff, M. & Adelberg, E. The Microbial World 2nd edn (Prentice Hall, 1963).Williams, T. A., Cox, C. J., Foster, P. G., Szöllősi, G. J. & Embley, T. M. Nat. Ecol. Evol. 4, 138–147 (2020).Article 

Google Scholar 
Steckbeck, W. Science 98, 487–488 (1943).CAS 
Article 

Google Scholar 
Creese, M. R. S. & Creese, T. M. Ladies in the Laboratory III (Scarecrow Press, 2010). More

Mariotte, P. et al. Plant–soil feedback: Bridging natural and agricultural sciences. Trends Ecol. Evol. 33, 129–142 (2018).PubMed 
Article 

Google Scholar 
Daryanto, S., Fu, B., Wang, L., Jacinthe, P. A. & Zhao, W. Quantitative synthesis on the ecosystem services of cover crops. Earth-Sci. Rev. 185, 357–373 (2018).ADS 
CAS 
Article 

Google Scholar 
Shackelford, G. E., Kelsey, R. & Dicks, L. V. Effects of cover crops on multiple ecosystem services: Ten meta-analyses of data from arable farmland in California and the Mediterranean. Land Use Policy 88, 104204 (2019).Article 

Google Scholar 
McDaniel, M. D., Tiemann, L. K. & Grandy, A. S. Does agricultural crop diversity enhance soil microbial biomass and organic matter dynamics? A meta-analysis. Ecol. Appl. 24, 560–570 (2014).CAS 
PubMed 
Article 

Google Scholar 
Wittwer, R. A., Dorn, B., Jossi, W. & van der Heijden, M. G. A. A. Cover crops support ecological intensification of arable cropping systems. Sci. Rep. 7, 41911 (2017).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Chahal, I. & Van Eerd, L. L. Cover crops increase tomato productivity and reduce nitrogen losses in a temperate humid climate. Nutr. Cycl. Agroecosyst. 119, 195–211 (2021).CAS 
Article 

Google Scholar 
Belfry, K. D., Trueman, C., Vyn, R. J., Loewen, S. A. & Van Eerd, L. L. Winter cover crops on processing tomato yield, quality, pest pressure, nitrogen availability, and profit margins. PLoS ONE 12, 1–17 (2017).Article 
CAS 

Google Scholar 
Wall, L. G. et al. Changes of paradigms in agriculture soil microbiology and new challenges in microbial ecology. Acta Oecologica 95, 68–73 (2019).ADS 
Article 

Google Scholar 
Schmidt, R., Gravuer, K., Bossange, A. V., Mitchell, J. & Scow, K. Long-term use of cover crops and no-till shift soil microbial community life strategies in agricultural soil. PLoS ONE 13, 1–19 (2018).
Google Scholar 
Schmidt, R., Mitchell, J. & Scow, K. Cover cropping and no-till increase diversity and symbiotroph:saprotroph ratios of soil fungal communities. Soil Biol. Biochem. 129, 99–109 (2019).CAS 
Article 

Google Scholar 
Ali, A. et al. Hiseq base molecular characterization of soil microbial community, diversity structure, and predictive functional profiling in continuous cucumber planted soil affected by diverse cropping systems in an intensive greenhouse region of Northern China. Int. J. Mol. Sci. 20, 2619 (2019).CAS 
PubMed Central 
Article 

Google Scholar 
Kim, N., Zabaloy, M. C., Guan, K. & Villamil, M. B. Do cover crops benefit soil microbiome? A meta-analysis of current research. Soil Biol. Biochem. 142, 107701 (2020).CAS 
Article 

Google Scholar 
Vukicevich, E., Lowery, T., Bowen, P., Úrbez-Torres, J. R. & Hart, M. Cover crops to increase soil microbial diversity and mitigate decline in perennial agriculture. A review. Agron. Sustain. Dev. 36, 1–14 (2016).CAS 
Article 

Google Scholar 
Nevins, C. J., Nakatsu, C. & Armstrong, S. Characterization of microbial community response to cover crop residue decomposition. Soil Biol. Biochem. 127, 39–49 (2018).CAS 
Article 

Google Scholar 
Peralta, A. L., Sun, Y., McDaniel, M. D. & Lennon, J. T. Crop rotational diversity increases disease suppressive capacity of soil microbiomes. Ecosphere 9, e02235 (2018).Article 

Google Scholar 
Cloutier, M. L. et al. Fungal community shifts in soils with varied cover crop treatments and edaphic properties. Sci. Rep. 10, 1–15 (2020).Article 
CAS 

Google Scholar 
Finney, D. M., Buyer, J. S. & Kaye, J. P. Living cover crops have immediate impacts on soil microbial community structure and function. J. Soil Water Conserv. 72, 361–373 (2017).Article 

Google Scholar 
Calderón, F. J., Nielsen, D., Acosta-Martínez, V., Vigil, M. F. & Lyon, D. Cover crop and irrigation effects on soil microbial communities and enzymes in semiarid agroecosystems of the central great plains of North America. Pedosphere 26, 192–205 (2016).Article 
CAS 

Google Scholar 
Romdhane, S. et al. Cover crop management practices rather than composition of cover crop mixtures affect bacterial communities in no-till agroecosystems. Front. Microbiol. 10, 1–11 (2019).Article 

Google Scholar 
Blanco-Canqui, H. & Lal, R. Crop residue removal impacts on soil productivity and environmental quality. CRC. Crit. Rev. Plant Sci. 28, 139–163 (2009).CAS 
Article 

Google Scholar 
Turmel, M. S., Speratti, A., Baudron, F., Verhulst, N. & Govaerts, B. Crop residue management and soil health: A systems analysis. Agric. Syst. 134, 6–16 (2015).Article 

Google Scholar 
Yang, Q., Wang, X. & Shen, Y. Comparison of soil microbial community catabolic diversity between rhizosphere and bulk soil induced by tillage or residue retention. J. Soil Sci. Plant Nutr. https://doi.org/10.4067/S0718-95162013005000017 (2013).Article 

Google Scholar 
Tang, H. et al. Tillage and crop residue incorporation effects on soil bacterial diversity in the double-cropping paddy field of southern China. Arch. Agron. Soil Sci. 67, 435–446 (2021).CAS 
Article 

Google Scholar 
Zhang, Y. et al. Long-term harvest residue retention could decrease soil bacterial diversities probably due to favouring oligotrophic lineages. Microb. Ecol. 76, 771–781 (2018).CAS 
PubMed 
Article 

Google Scholar 
Zhang, C. et al. Straw retention efficiently improves fungal communities and functions in the fallow ecosystem. BMC Microbiol. 21, 52 (2021).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Chahal, I. & Van Eerd, L. L. Cover crop and crop residue removal effects on temporal dynamics of soil carbon and nitrogen in a temperate, humid climate. PLoS ONE 15, e0235665 (2020).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Chahal, I. & Van Eerd, L. L. Evaluation of commercial soil health tests using a medium-term cover crop experiment in a humid, temperate climate. Plant Soil 427, 351–367 (2018).CAS 
Article 

Google Scholar 
Ruis, S. J. & Blanco-Canqui, H. Cover crops could offset crop residue removal effects on soil carbon and other properties: A review. Agron. J. 109, 1785–1805 (2017).CAS 
Article 

Google Scholar 
Zhao, M. et al. Intercropping affects genetic potential for inorganic nitrogen cycling by root-associated microorganisms in Medicago sativa and Dactylis glomerata. Appl. Soil Ecol. 119, 260–266 (2017).ADS 
Article 

Google Scholar 
Wardle, D. A. et al. Ecological linkages between aboveground and belowground biota. Science (80-). 304, 1629–1633 (2004).ADS 
CAS 
Article 

Google Scholar 
Xiong, C. et al. Host selection shapes crop microbiome assembly and network complexity. New Phytol. 229, 1091–1104 (2021).CAS 
PubMed 
Article 

Google Scholar 
McDaniel, M. D., Grandy, A. S., Tiemann, L. K. & Weintraub, M. N. Eleven years of crop diversification alters decomposition dynamics of litter mixtures incubated with soil. Ecosphere 7, e01426 (2016).Article 

Google Scholar 
Buyer, J. S., Teasdale, J. R., Roberts, D. P., Zasada, I. A. & Maul, J. E. Factors affecting soil microbial community structure in tomato cropping systems. Soil Biol. Biochem. 42, 831–841 (2010).CAS 
Article 

Google Scholar 
Fernandez-Gnecco, G. et al. Microbial community analysis of soils under different soybean cropping regimes in the Argentinean south-eastern Humid Pampas. FEMS Microbiol. Ecol. 97, 1–14 (2021).Article 
CAS 

Google Scholar 
Semenov, M. V., Krasnov, G. S., Semenov, V. M. & van Bruggen, A. H. C. Long-term fertilization rather than plant species shapes rhizosphere and bulk soil prokaryotic communities in agroecosystems. Appl. Soil Ecol. 154, 103641 (2020).Article 

Google Scholar 
White, C. M. & Weil, R. R. Forage radish cover crops increase soil test phosphorus surrounding radish taproot holes. Soil Sci. Soc. Am. J. 75, 121–130 (2011).ADS 
CAS 
Article 

Google Scholar 
Schulz, M., Marocco, A., Tabaglio, V., Macias, F. A. & Molinillo, J. M. G. Benzoxazinoids in rye allelopathy—From discovery to application in sustainable weed control and organic farming. J. Chem. Ecol. 39, 154–174 (2013).CAS 
PubMed 
Article 

Google Scholar 
Cheng, F. & Cheng, Z. Research progress on the use of plant allelopathy in agriculture and the physiological and ecological mechanisms of allelopathy. Front. Plant Sci. 6, 1020 (2015).PubMed 
PubMed Central 

Google Scholar 
Thapa, V. R., Ghimire, R., Acosta-Martínez, V., Marsalis, M. A. & Schipanski, M. E. Cover crop biomass and species composition affect soil microbial community structure and enzyme activities in semiarid cropping systems. Appl. Soil Ecol. 157, 103735 (2021).Article 

Google Scholar 
Drost, S. M., Rutgers, M., Wouterse, M., de Boer, W. & Bodelier, P. L. E. Decomposition of mixtures of cover crop residues increases microbial functional diversity. Geoderma 361, 114060 (2020).ADS 
CAS 
Article 

Google Scholar 
Di Rauso Simeone, G., Müller, M., Felgentreu, C. & Glaser, B. Soil microbial biomass and community composition as affected by cover crop diversity in a short-term field experiment on a podzolized Stagnosol-Cambisol. J. Plant Nutr. Soil Sci. 183, 539–549 (2020).Article 
CAS 

Google Scholar 
Maul, J. E. et al. Microbial community structure and abundance in the rhizosphere and bulk soil of a tomato cropping system that includes cover crops. Appl. Soil Ecol. 77, 42–50 (2014).Article 

Google Scholar 
Huang, J. et al. Allocation and turnover of rhizodeposited carbon in different soil microbial groups. Soil Biol. Biochem. 150, 107973 (2020).CAS 
Article 

Google Scholar 
Strickland, M. S. & Rousk, J. Considering fungal:bacterial dominance in soils—Methods, controls, and ecosystem implications. Soil Biol. Biochem. 42, 1385–1395 (2010).CAS 
Article 

Google Scholar 
Leff, J. W. et al. Predicting the structure of soil communities from plant community taxonomy, phylogeny, and traits. ISME J. 12, 1794–1805 (2018).PubMed 
PubMed Central 
Article 

Google Scholar 
Milcu, A. et al. Functionally and phylogenetically diverse plant communities key to soil biota. Ecology 94, 1878–1885 (2013).PubMed 
Article 

Google Scholar 
Lozupone, C. A., Hamady, M., Kelley, S. T. & Knight, R. Quantitative and qualitative β diversity measures lead to different insights into factors that structure microbial communities. Appl. Environ. Microbiol. 73, 1576–1585 (2007).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Lay, C.-Y., Hamel, C. & St-Arnaud, M. Taxonomy and pathogenicity of Olpidium brassicae and its allied species. Fungal Biol. 122, 837–846 (2018).PubMed 
Article 

Google Scholar 
Liu, L., Zhu, K., Wurzburger, N. & Zhang, J. Relationships between plant diversity and soil microbial diversity vary across taxonomic groups and spatial scales. Ecosphere 11, e02999 (2020).
Google Scholar 
Hartwright, L. M., Hunter, P. J. & Walsh, J. A. A comparison of Olpidium isolates from a range of host plants using internal transcribed spacer sequence analysis and host range studies. Fungal Biol. 114, 26–33 (2010).CAS 
PubMed 
Article 

Google Scholar 
Barel, J. M. et al. Winter cover crop legacy effects on litter decomposition act through litter quality and microbial community changes. J. Appl. Ecol. 56, 132–143 (2019).CAS 
Article 

Google Scholar 
Austin, E. E., Wickings, K., McDaniel, M. D., Robertson, G. P. & Grandy, A. S. Cover crop root contributions to soil carbon in a no-till corn bioenergy cropping system. GCB Bioenergy 9, 1252–1263 (2017).CAS 
Article 

Google Scholar 
Bai, Z., Liang, C., Bodé, S., Huygens, D. & Boeckx, P. Phospholipid 13C stable isotopic probing during decomposition of wheat residues. Appl. Soil Ecol. 98, 65–74 (2016).Article 

Google Scholar 
Põlme, S. et al. FungalTraits: A user-friendly traits database of fungi and fungus-like stramenopiles. Fungal Divers. 105, 1–16 (2020).Article 

Google Scholar 
Pepe, A., Giovannetti, M. & Sbrana, C. Lifespan and functionality of mycorrhizal fungal mycelium are uncoupled from host plant lifespan. Sci. Rep. 8, 1–10 (2018).
Google Scholar 
Frey, S. D. Mycorrhizal fungi as mediators of soil organic matter dynamics. Annu. Rev. Ecol. Evol. Syst. 50, 237–259 (2019).Article 

Google Scholar 
Saleem, M., Hu, J. & Jousset, A. More than the sum of its parts: Microbiome biodiversity as a driver of plant growth and soil health. Annu. Rev. Ecol. Evol. Syst. 50, 145–168 (2019).Article 

Google Scholar 
Wei, Z. et al. Initial soil microbiome composition and functioning predetermine future plant health. Sci. Adv. 5, 1–12 (2019).
Google Scholar 
Ozimek, E. & Hanaka, A. Mortierella species as the plant growth-promoting fungi present in the agricultural soils. Agriculture 11, 7 (2020).Article 
CAS 

Google Scholar 
Li, F. et al. Mortierella elongata’s roles in organic agriculture and crop growth promotion in a mineral soil. L. Degrad. Dev. 29, 1642–1651 (2018).Article 

Google Scholar 
Sansinenea, E. Bacillus spp.: As plant growth-promoting bacteria. in Secondary Metabolites of Plant Growth Promoting Rhizomicroorganisms: Discovery and Applications 225–237 (Springer, 2019). https://doi.org/10.1007/978-981-13-5862-3_11.Palaniyandi, S. A., Yang, S. H., Zhang, L. & Suh, J.-W. Effects of actinobacteria on plant disease suppression and growth promotion. Appl. Microbiol. Biotechnol. 97, 9621–9636 (2013).CAS 
PubMed 
Article 

Google Scholar 
Jung, M.-Y. et al. Ammonia-oxidizing archaea possess a wide range of cellular ammonia affinities. ISME J. 16, 272–283 (2022).CAS 
PubMed 
Article 

Google Scholar 
Zhong, Y. et al. Microbial community assembly and metabolic function during wheat straw decomposition under different nitrogen fertilization treatments. Biol. Fertil. Soils 56, 697–710 (2020).CAS 
Article 

Google Scholar 
Liu, X. et al. Decomposing cover crops modify root-associated microbiome composition and disease tolerance of cash crop seedlings. Soil Biol. Biochem. 160, 108343 (2021).CAS 
Article 

Google Scholar 
Larkin, R. P., Griffin, T. S. & Honeycutt, C. W. Rotation and cover crop effects on soilborne potato diseases, tuber yield, and soil microbial communities. Plant Dis. 94, 1491–1502 (2010).PubMed 
Article 

Google Scholar 
van der Putten, W. H., Bradford, M. A., Brinkman, E. P., van de Voorde, T. F. J. & Veen, G. F. Where, when and how plant–soil feedback matters in a changing world. Funct. Ecol. 30, 1109–1121 (2016).Article 

Google Scholar 
Menalled, U. D., Seipel, T. & Menalled, F. D. Farming system effects on biologically mediated plant–soil feedbacks. Renew. Agric. Food Syst. 36, 1–7 (2021).Article 

Google Scholar 
Fierer, N. & Jackson, J. Assessment of soil microbial community structure by use of taxon-specific quantitative PCR assays. Appl. Environ. Microbiol. 71, 4117 (2005).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Vainio, E. J. & Hantula, J. Direct analysis of wood-inhabiting fungi using denaturing gradient gel electrophoresis of amplified ribosomal DNA. Mycol. Res. 104, 927–936 (2000).CAS 
Article 

Google Scholar 
Caporaso, J. G. et al. Global patterns of 16S rRNA diversity at a depth of millions of sequences per sample. Proc. Natl. Acad. Sci. 108, 4516–4522 (2011).ADS 
CAS 
PubMed 
Article 

Google Scholar 
White, T. J., Bruns, T., Lee, S. & Taylor, J. Amplification and direct sequencing of fungal ribosomal RNA genes for phylogenetics. In PCR Protocols: A Guide to Methods and Applications (eds Innis, M. A. et al.) 315–322 (Academic Press, 1990).

Google Scholar 
Bolyen, E. et al. Reproducible, interactive, scalable and extensible microbiome data science using QIIME 2. Nat. Biotechnol. 37, 852–857 (2019).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Callahan, B. J. et al. DADA2: High-resolution sample inference from Illumina amplicon data. Nat. Methods 13, 581–583 (2016).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Rivers, A. R., Weber, K. C., Gardner, T. G., Liu, S. & Armstrong, S. D. ITSxpress: Software to rapidly trim internally transcribed spacer sequences with quality scores for marker gene analysis. F1000Research 7, 1418 (2018).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Katoh, K., Misawa, K., Kuma, K. & Miyata, T. MAFFT: A novel method for rapid multiple sequence alignment based on fast Fourier transform. Nucleic Acids Res. 30, 3059–3066 (2002).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Price, M. N., Dehal, P. S. & Arkin, A. P. FastTree 2—Approximately maximum-likelihood trees for large alignments. PLoS ONE 5, e9490 (2010).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Bokulich, N. A. et al. Optimizing taxonomic classification of marker-gene amplicon sequences with QIIME 2’s q2-feature-classifier plugin. Microbiome 6, 90 (2018).PubMed 
PubMed Central 
Article 

Google Scholar 
Quast, C. et al. The SILVA ribosomal RNA gene database project: Improved data processing and web-based tools. Nucleic Acids Res. 41, D590–D596 (2013).CAS 
PubMed 
Article 

Google Scholar 
Abarenkov, K. et al. UNITE QIIME release for Fungi. https://doi.org/10.15156/bio/786385 (2020).R Core Team. R: A Language and Environment for Statistical Computing. (2020).Oksanen, J. et al. vegan: Community Ecology Package. (2020).Anderson, M. J., Gorley, R. N. & Clarke, K. R. PERMANOVA+ for PRIMER: Guide to Software and Statistical Methods. (PRIMER-E, 2008).Anderson, M. J. & Willis, T. J. Canonical analysis of principal coordinates: A useful method of constrained ordination for ecology. Ecology 84, 511–525 (2003).Article 

Google Scholar 
Cáceres, M. D. & Legendre, P. Associations between species and groups of sites: Indices and statistical inference. Ecology 90, 3566–3574 (2009).PubMed 
Article 

Google Scholar 
Fernandes, A. D. et al. Unifying the analysis of high-throughput sequencing datasets: Characterizing RNA-seq, 16S rRNA gene sequencing and selective growth experiments by compositional data analysis. Microbiome 2, 15 (2014).PubMed 
PubMed Central 
Article 

Google Scholar  More

Marzluff, J. M. Worldwide urbanization and its effects on birds. In Avian Ecology and Conservation in an Urbanizing World (eds Marzluff, J. et al.) 19–47 (Springer, Boston, 2001).Chapter 

Google Scholar 
McKinney, M. L. Effects of urbanization on species richness: A review of plants and animals. Urban Ecosyst. 11, 161–176 (2008).Article 

Google Scholar 
Luniak, M. Synurbization–adaptation of animal wildlife to urban development in Proceedings 4th international urban wildlife symposium (eds. Shaw, W., Harris, L.,Vandruff, L.) 50–55 (University of Arizona, Tucson, ARI, 2004).Isaksson, C. Impact of urbanization on birds in Bird Species how they arise, modify and vanish (ed. Tietze D. T.) 235–257 (Springer, 2018).Minias, P. Successful colonization of a novel urban environment is associated with an urban behavioural syndrome in a reed-nesting waterbird. Ethology 121, 1178–1190 (2015).Article 

Google Scholar 
Møller, A. P. et al. Urban habitats and feeders both contribute to flight initiation distance reduction in birds. Behav. Ecol. 26, 861–865 (2015).Article 

Google Scholar 
Jokimäki, J. & Suhonen, J. Distribution and habitat selection of wintering birds in urban environments. Landsc. Urban Plan. 39, 253–263 (1998).Article 

Google Scholar 
Francis, R. A. & Chadwick, M. A. What makes a species synurbic?. Appl. Geogr. 32, 514–521 (2012).Article 

Google Scholar 
Møller, A. P. et al. High urban population density of birds reflects their timing of urbanization. Oecologia 170, 867–875 (2012).PubMed 
Article 
ADS 

Google Scholar 
Tella, J. L. et al. Offspring body condition and immunocompetence are negatively affected by high breeding densities in a colonial seabird: A multiscale approach. Proc. R. Soc. B 268, 1455–1461 (2001).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Savoca, M. S., Bonter, D. N., Zuckerberg, B., Dickinson, J. L. & Ellis, J. C. Nesting density is an important factor affecting chick growth and survival in the Herring Gull. Condor 113, 565–571 (2011).Article 

Google Scholar 
Minias, P., Włodarczyk, R. & Janiszewski, T. Opposing selective pressures may act on the colony size in a waterbird species. Evol. Ecol. 29, 283–297 (2015).Article 

Google Scholar 
Kamiński, M. et al. Density-dependence of nestling immune function and physiological condition in semi-precocial colonial bird: A cross-fostering experiment. Front. Zool. 18, 7 (2021).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Ward, P. & Zahavi, A. The importance of certain assemblages of birds as “information-centres” for food-finding. Ibis 115, 517–534 (1973).Article 

Google Scholar 
Danchin, E. & Wagner, R. H. The evolution of coloniality: The emergence of new perspectives. Trends Ecol. Evol. 12, 342–347 (1997).CAS 
PubMed 
Article 

Google Scholar 
Brown, C. R. & Brown, M. B. Coloniality in the Cliff Swallow: The Effect of Group Size on Social Behavior (University of Chicago Press, 1996).
Google Scholar 
Evans, J. C., Votier, S. C. & Dall, S. R. Information use in colonial living. Biol. Rev. 91, 658–672 (2016).PubMed 
Article 

Google Scholar 
Brown, C. R. & Brown, M. B. Avian coloniality. In Current Ornithology (eds Brown, C. R. & Brown, M. B.) 1–82 (Springer, Boston, 2001).
Google Scholar 
Coulson, J. C., Duncan, N. & Thomas, C. Changes in the breeding biology of the herring gull (Larus argentatus) induced by reduction in the size and density of the colony. J. Anim. Ecol. 51, 739–756 (1982).Article 

Google Scholar 
Ots, I. & Horak, P. Great tits Parus major trade health for reproduction. Proc. R. Soc. B. 263, 1443–1447 (1996).CAS 
PubMed 
Article 
ADS 

Google Scholar 
Richner, H. & Tripet, F. Ectoparasitism and the trade-off between current and future reproduction. Oikos 86, 535–538 (1999).Article 

Google Scholar 
Fokkema, R. W., Ubels, R. & Tinbergen, J. M. Great tits trade off future competitive advantage for current reproduction. Behav. Ecol. 27, 1656–1664 (2016).
Google Scholar 
Horak, P. & Leberton, J. D. Survival of adult Great Tits Parus major in relation to sex and habitat; a comparison of urban and rural populations. Ibis 140, 205–209 (1998).Article 

Google Scholar 
Stracey, C. M. & Robinson, S. K. Are urban habitats ecological traps for a native songbird? Season-long productivity, apparent survival, and site fidelity in urban and rural habitats. J. Avian Biol. 43, 50–60 (2012).Article 

Google Scholar 
Sepp, T., McGraw, K. J., Kaasik, A. & Giraudeau, M. A review of urban impacts on avian life-history evolution: Does city living lead to slower pace of life?. Glob. Change Biol. 24, 1452–1469 (2018).Article 
ADS 

Google Scholar 
Phillips, J. N., Gentry, K. E., Luther, D. A. & Derryberry, E. P. Surviving in the city: Higher apparent survival for urban birds but worse condition on noisy territories. Ecosphere 9, e02440 (2018).Article 

Google Scholar 
Johnston, R. F. & Janiga, M. Feral Pigeons (Oxford University Press on Demand, 1995).
Google Scholar 
Giunchi, D., Mucci, N., Bigi, D., Mengoni, C. & Baldaccini, N. E. Feral pigeon populations: Their gene pool and links with local domestic breeds. Zoology 142, 125817 (2020).PubMed 
Article 

Google Scholar 
Sol, D. Artificial selection, naturalization, and fitness: Darwin’s pigeons revisited. Biol. J. Linn. Soc. 93, 657–665 (2008).Article 

Google Scholar 
Giunchi, D., Albores-Barajas, Y. V., Baldaccini, N. E., Vanni, L. & Soldatini, C. Feral pigeons: Problems, dynamics and control methods. In Integrated Pest Management and Pest Control. Current and Future Tactics (eds Soloneski, S. & Larramendy, M.) 215–240 (InTechOpen, London, 2012).
Google Scholar 
Senar, J. C., Navalpotro, H., Pascual, J. & Montalvo, T. Nicarbazin has no effect on reducing feral pigeon populations in Barcelona. Pest Manag. Sci. 77, 131–137 (2021).CAS 
PubMed 
Article 

Google Scholar 
Rose, E., Nagel, P. & Haag-Wackernagel, D. Spatio-temporal use of the urban habitat by feral pigeons (Columba livia). Behav. Ecol. Sociobiol. 60, 242–254 (2006).Article 

Google Scholar 
Corbel, H. et al. Stress response varies with plumage colour and local habitat in feral pigeons. J. Ornithol. 157, 825–837 (2016).Article 

Google Scholar 
Møller, A. P., Merino, S., Brown, C. R. & Robertson, R. J. Immune defense and host sociality: A comparative study of swallows and martins. Am. Nat. 158, 136–145 (2001).PubMed 
Article 

Google Scholar 
Drzewińska-Chańko, J. et al. Immunocompetent birds choose larger breeding colonies. J. Anim. Ecol. https://doi.org/10.1111/1365-2656.13540 (2021).Article 
PubMed 

Google Scholar 
Saino, N., Suffritti, C., Martinelli, R., Rubolini, D. & Møller, A. P. Immune response covaries with corticosterone plasma levels under experimentally stressful conditions in nestling barn swallows (Hirundo rustica). Behav. Ecol. 14, 318–325 (2003).Article 

Google Scholar 
Goutte, A. et al. Long-term survival effect of corticosterone manipulation in black-legged kittiwakes. Gen. Comp. Endocrinol. 167, 246–251 (2010).CAS 
PubMed 
Article 

Google Scholar 
Møller, A. P., Christe, P., Erritzøe, J. & Mavarez, J. Condition, disease and immune defence. Oikos 83, 301–306 (1998).Article 

Google Scholar 
Navarro, C., Marzal, A., De Lope, F. & Møller, A. P. Dynamics of an immune response in house sparrows Passer domesticus in relation to time of day, body condition and blood parasite infection. Oikos 101, 291–298 (2003).Article 

Google Scholar 
Toïgo, C., Gaillard, J. M., Van Laere, G., Hewison, M. & Morellet, N. How does environmental variation influence body mass, body size, and body condition? Roe deer as a case study. Ecography 29, 301–308 (2006).Article 

Google Scholar 
Jacquin, L. et al. A potential role for parasites in the maintenance of color polymorphism in urban birds. Oecologia 173, 1089–1099 (2013).CAS 
PubMed 
Article 
ADS 

Google Scholar 
Meillère, A., Brischoux, F., Parenteau, C. & Angelier, F. Influence of urbanization on body size, condition, and physiology in an urban exploiter: A multi-component approach. PLoS ONE https://doi.org/10.1371/journal.pone.0135685 (2015).Article 
PubMed 
PubMed Central 

Google Scholar 
Peig, J. & Green, A. J. New perspectives for estimating body condition from mass/length data: The scaled mass index as an alternative method. Oikos 118, 1883–1891 (2009).Article 

Google Scholar 
Jacquin, L. et al. Melanin-based coloration is related to parasite intensity and cellular immune response in an urban free living bird: The feral pigeon Columba livia. J. Avian Biol. 42, 11–15 (2011).Article 

Google Scholar 
Liker, A., Papp, Z., Bókony, V. & Lendvai, A. Z. Lean birds in the city: Body size and condition of house sparrows along the urbanization gradient. J. Anim. Ecol. 77, 789–795 (2008).CAS 
PubMed 
Article 

Google Scholar 
Audet, J. N., Ducatez, S. & Lefebvre, L. The town bird and the country bird: Problem solving and immunocompetence vary with urbanization. Behav. Ecol. 27, 637–644 (2016).Article 

Google Scholar 
Kurucz, K., Purger, J. J. & Batáry, P. Urbanization shapes bird communities and nest survival, but not their food quantity. Glob. Ecol. Conserv. 26, e01475 (2021).Article 

Google Scholar 
Partecke, J., Schwabl, I. & Gwinner, E. Stress and the city: Urbanization and its effects on the stress physiology in European blackbirds. Ecology 87, 1945–1952 (2006).PubMed 
Article 

Google Scholar 
Bailly, J. et al. Negative impact of urban habitat on immunity in the great tit Parus major. Oecologia 182, 1053–1062 (2016).PubMed 
Article 
ADS 

Google Scholar 
Glądalski, M. et al. Differences in use of bryophyte species in tit nests between two contrasting habitats: An urban park and a forest. Eur. Zool. J. 88, 807–815 (2021).Article 

Google Scholar 
Tella, J. L., Scheuerlein, A. & Ricklefs, R. E. Is cell–mediated immunity related to the evolution of life-history strategies in birds?. Proc. R. Soc. B 269, 1059–1066 (2002).PubMed 
PubMed Central 
Article 

Google Scholar 
Brown, C. R. & Brown, M. B. Empirical measurement of parasite transmission between groups in a colonial bird. Ecology 85, 1619–1626 (2004).Article 

Google Scholar 
O’Brien, V. A. & Brown, C. R. Group size and nest spacing affect Buggy Creek virus (Togaviridae: Alphavirus) infection in nestling house sparrows. PLoS ONE 6, e25521 (2011).CAS 
PubMed 
PubMed Central 
Article 
ADS 

Google Scholar 
Wilcoxen, T. E. et al. Effects of bird-feeding activities on the health of wild birds. Conserv. Physiol. 3, cov058 (2015).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Moyers, S. C., Adelman, J. S., Farine, D. R., Thomason, C. A. & Hawley, D. M. Feeder density enhances house finch disease transmission in experimental epidemics. Philos. Trans. R. Soc. B 373, 20170090 (2018).Article 
CAS 

Google Scholar 
Møller, A. P. Successful city dwellers: A comparative study of the ecological characteristics of urban birds in the Western Palearctic. Oecologia 159, 849–858 (2009).PubMed 
Article 
ADS 

Google Scholar 
Watson, H., Videvall, E., Andersson, M. N. & Isaksson, C. Transcriptome analysis of a wild bird reveals physiological responses to the urban environment. Sci. Rep. 7, 44180 (2017).CAS 
PubMed 
PubMed Central 
Article 
ADS 

Google Scholar 
Hasselquist, D. & Nilsson, J. Å. Physiological mechanisms mediating costs of immune responses: What can we learn from studies of birds?. Anim. Behav. 83, 1303–1312 (2012).Article 

Google Scholar 
Biard, C., Monceau, K., Motreuil, S. & Moreau, J. Interpreting immunological indices: The importance of taking parasite community into account. An example in blackbirds Turdus merula. Methods Ecol. Evol. 6, 960–972 (2015).Article 

Google Scholar 
Leclaire, S., Czirják, G. Á., Hammouda, A. & Gasparini, J. Feather bacterial load shapes the trade-off between preening and immunity in pigeons. BMC Evol. Biol. 15, 60 (2015).PubMed 
PubMed Central 
Article 

Google Scholar 
Vinkler, M., Adelman, J. S. & Ardia, D. R. Evolutionary and ecological immunology. In Avian Immunology 3rd edn (eds Kaspers, B. et al.) 519–558 (Academic Press, London, 2021).
Google Scholar 
Davis, A. K., Maney, D. L. & Maerz, J. C. The use of leukocyte profiles to measure stress in vertebrates: A review for ecologists. Funct. Ecol. 22, 760–772 (2008).Article 

Google Scholar 
Indykiewicz, P., Podlaszczuk, P., Kamiński, M., Włodarczyk, R. & Minias, P. Central–periphery gradient of individual quality within a colony of Black-headed Gulls. Ibis 161, 744–758 (2019).Article 

Google Scholar 
Vleck, C. M., Vertalino, N., Vleck, D. & Bucher, T. L. Stress, corticosterone, and heterophil to lymphocyte ratios in free-living Adélie penguins. Condor 102, 392–400 (2000).Article 

Google Scholar 
Davis, A. K., Cook, K. C. & Altizer, S. Leukocyte profiles in wild house finches with and without mycoplasmal conjunctivitis, a recently emerged bacterial disease. EcoHealth 1, 362–373 (2004).Article 

Google Scholar 
Lobato, E., Moreno, J., Merino, S., Sanz, J. J. & Arriero, E. Haematological variables are good predictors of recruitment in nestling pied flycatchers (Ficedula hypoleuca). Ecoscience 12, 27–34 (2005).Article 

Google Scholar 
Bobby Fokidis, H., Greiner, E. C. & Deviche, P. Interspecific variation in avian blood parasites and haematology associated with urbanization in a desert habitat. J. Avian Biol. 39, 300–310 (2008).Article 

Google Scholar 
Padgett, D. A. & Glaser, R. How stress influences the immune response. Trends Immunol. 24, 444–448 (2003).CAS 
PubMed 
Article 

Google Scholar 
Dimitrov, S. et al. Cortisol and epinephrine control opposing circadian rhythms in T cell subsets. Blood 113, 5134–5143 (2009).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Ilmonen, P., Hasselquist, D., Langefors, Å. & Wiehn, J. Stress, immunocompetence and leukocyte profiles of pied flycatchers in relation to brood size manipulation. Oecologia 136, 148–154 (2003).PubMed 
Article 
ADS 

Google Scholar 
Minias, P., Gach, K., Włodarczyk, R. & Janiszewski, T. Colony size affects nestling immune function: A cross-fostering experiment in a colonial waterbird. Oecologia 190, 333–341 (2019).PubMed 
PubMed Central 
Article 
ADS 

Google Scholar 
Cyr, N. E., Earle, K., Tam, C. & Romero, L. M. The effect of chronic psychological stress on corticosterone, plasma metabolites, and immune responsiveness in European starlings. Gen. Comp. Endocrinol. 154, 59–66 (2007).CAS 
PubMed 
Article 

Google Scholar 
Schoech, S. J., Bowman, R. & Reynolds, S. J. Food supplementation and possible mechanisms underlying early breeding in the Florida Scrub-Jay (Aphelocoma coerulescens). Horm. Behav. 46, 565–573 (2004).CAS 
PubMed 
Article 

Google Scholar 
Ibáñez-Álamo, J. D. et al. Physiological stress does not increase with urbanization in European blackbirds: Evidence from hormonal, immunological and cellular indicators. Sci. Total Environ. 721, 137332 (2020).PubMed 
Article 
ADS 
CAS 

Google Scholar 
Bonier, F. Hormones in the city: Endocrine ecology of urban birds. Horm. Behav. 61, 763–772 (2012).CAS 
PubMed 
Article 

Google Scholar 
Valdebenito, J. O. et al. Seasonal variation in sex-specific immunity in wild birds. Sci. Rep. 11, 1349 (2021).CAS 
PubMed 
PubMed Central 
Article 
ADS 

Google Scholar 
Hetmański, T. Timing of breeding in the Feral Pigeon Columba livia f. domestica in Słupsk (NW Poland). Acta Ornithol. 39, 105–110 (2004).Article 

Google Scholar 
Dijkstra, C. et al. An adaptive annual rhythm in the sex of first pigeon eggs. Behav. Ecol. Sociobiol. 64, 1393–1402 (2010).PubMed 
PubMed Central 
Article 

Google Scholar 
Swanson, D. L. Seasonal variation of vascular oxygen transport in the dark-eyed junco. Condor 92, 62–66 (1990).Article 

Google Scholar 
Niedojadlo, J., Bury, A., Cichoń, M., Sadowska, E. T. & Bauchinger, U. Lower haematocrit, haemoglobin and red blood cell number in zebra finches acclimated to cold compared to thermoneutral temperature. J. Avian Biol. 49, e01596 (2018).Article 

Google Scholar 
Roulin, A. Condition-dependence, pleiotropy and the handicap principle of sexual selection in melanin-based colouration. Biol. Rev. 91, 328–348 (2016).PubMed 
Article 

Google Scholar 
Statistics Poland. https://stat.gov.pl/en/ (2021).Sol, D. & Senar, J. C. Urban pigeon populations: Stability, home range, and the effect of removing individuals. Can. J. Zool. 73, 1154–1160 (1995).Article 

Google Scholar 
Minias, P. Reproduction and survival in the city: Which fitness components drive urban colonization in a reed-nesting waterbird?. Curr. Zool. 62, 79–87 (2016).PubMed 
PubMed Central 
Article 

Google Scholar 
Meissner, W. & Fischer, I. Sexing of common gull, Larus canus, using linear measurements. Folia Zool. 66, 183–188 (2017).Article 

Google Scholar 
Haag-Wackernagel, D., Heeb, P. & Leiss, A. Phenotype-dependent selection of juvenile urban feral pigeons Columba livia. Bird Study 53, 163–170 (2006).Article 

Google Scholar 
Harter, T. S., Reichert, M., Brauner, C. J. & Milsom, W. K. Validation of the i-STAT and HemoCue systems for the analysis of blood parameters in the bar-headed goose, Anser indicus. Conserv. Physiol. 3, cov021 (2015).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Minias, P. The use of haemoglobin concentrations to assess physiological condition in birds: A review. Conserv. Physiol. 3, cov007 (2015).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
Martin, L. B. et al. Phytohemagglutinin-induced skin swelling in birds: Histological support for a classic immunoecological technique. Funct. Ecol. 20, 290–299 (2006).Article 

Google Scholar 
Brown, G. P., Shilton, C. M. & Shine, R. Measuring amphibian immunocompetence: Validation of the phytohemagglutinin skin-swelling assay in the cane toad, Rhinella marina. Methods Ecol. Evol. 2, 341–348 (2011).Article 

Google Scholar 
Kennedy, M. W. & Nager, R. G. The perils and prospects of using phytohaemagglutinin in evolutionary ecology. Trends Ecol. Evol. 21, 653–655 (2006).PubMed 
Article 

Google Scholar 
Vinkler, M., Bainová, H. & Albrecht, T. Functional analysis of the skin-swelling response to phytohaemagglutinin. Funct. Ecol. 24, 1081–1086 (2010).Article 

Google Scholar 
Turmelle, A. S., Ellison, J. A., Mendonça, M. T. & McCracken, G. F. Histological assessment of cellular immune response to the phytohemagglutinin skin test in Brazilian free-tailed bats (Tadarida brasiliensis). J. Comp. Physiol. B 180, 1155–1164 (2010).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Griffiths, R., Double, M. C., Orr, K. & Dawson, R. J. A DNA test to sex most birds. Mol. Ecol. 7, 1071–1075 (1998).CAS 
PubMed 
Article 

Google Scholar 
Çakmak, E., Akın Pekşen, Ç. & Bilgin, C. C. Comparison of three different primer sets for sexing birds. J. Vet. Diagn. Investig. 29, 59–63 (2017).Article 
CAS 

Google Scholar 
Kaiser, H. F. The application of electronic computers to factor analysis. Educ. Psychol. Meas. 20, 141–151 (1960).Article 

Google Scholar 
Bates, D., Mächler, M., Bolker, B. & Walker, S. Fitting linear mixed-effects models using lme4. J. Stat. Softw. 67, 1–48 (2015).Article 

Google Scholar 
Kuznetsova, A., Brockhoff, P. B. & Christensen, R. H. B. lmerTest package: Tests in linear mixed effects models. J. Stat. Softw. 82, 1–26 (2017).Article 

Google Scholar 
Jaeger, B. C., Edwards, L. J., Das, K. & Sen, P. K. An R 2 statistic for fixed effects in the generalized linear mixed model. J. Appl. Stat. 44, 1086–1105 (2017).MathSciNet 
MATH 
Article 

Google Scholar 
Johnson, P. C. Extension of Nakagawa & Schielzeth’s R2GLMM to random slopes models. Methods Ecol. Evol. 5, 944–946 (2014).PubMed 
PubMed Central 
Article 

Google Scholar 
Bartoń, K. MuMIn: Multi-model inference. R package ver. 1.43.17. CRAN: The Comprehensive R Archive Network, Berkeley, CA, USA. https://CRAN.R-project.org/package=MuMIn (2020).Wickham, H. ggplot2: Elegant Graphics for Data Analysis (Springer, 2016).MATH 
Book 

Google Scholar 
Kahle, D. J. & Wickham, H. ggmap: Spatial visualization with ggplot2. R J. 5, 144–161 (2013).Article 

Google Scholar  More

Data UsedWe used cadastral maps from 1858 to reconstruct the LULC structure of Aksu and Kestel for the mid-nineteenth century. General Staff of the Ottoman Army produced these maps in 1:10,000 scale. These maps were one of the earliest attempts of creating cadastral maps in the Ottoman Empire. The images of historical maps scanned at 1270 dpi resolutions are provided by the Turkish Presidency State Archives of the Republic of Turkey – Department of Ottoman Archives (Map collection, HRT.h, 561–567). Individual tiles of cadastral maps are of a 1:2,000 scale. However, these maps are kept separated from their accompanying cadastral registers or documentation regarding their production process in the archives. There are no studies on the production of these cadastral maps, but few studies used them until now35,36.The LULC structures of Aksu and Kestel for the mid-twentieth century were generated using aerial photographs from June 23, 1955, with a scale of 1:30,000. These aerial photographs were captured by the US Navy Photographic Squadron VJ-62 (established on April 10, 1952, re-designated to VAP-62 on July 1956, and disestablished on October 15, 1969). The squadron conducted mapping and special photographic projects worldwide37. Lastly, the VHR satellite images of WorldView-3 (WV-3) with 0.3 m of spatial resolution, collected on September 6, 2020, were used to show the up-to-date LULC patterns of Aksu and Kestel.MethodologyFigure 2 shows the flowchart of steps followed in this study to detect the LULC changes. The workflow includes three phases: preprocessing, LULC mapping, and statistical analysis of LULC changes.Figure 2Flowchart of the processing steps for the LULC change analysis for Kestel.Full size imageData preprocessingOrthorectification is the first important step in ensuring accurate spatial positioning among the multi-temporal and multi-source images, eliminating geometric distortions, and defining all data sets on a common projection system. To align the multi-modal geospatial datasets, we first performed the orthorectification of the satellite images and then we used the orthorectified satellite images as reference for the georeferencing of the cadastral maps and aerial photographs.Satellite imagery orthorectificationWe first pan-sharpened the WV-3 images by applying the PANSHARP2 algorithm38 to fuse the panchromatic (PAN) image of 0.3 m spatial resolution with four multispectral bands (R, G, B, and near-infrared (NIR)) of 1.2 m. We then geometrically corrected the pan-sharpened (PSP) WV-3 imageries using an ALOS Global Digital Surface Model with a horizontal resolution of approximately 30 m (ALOS World 3D – 30 m), rational polynomial coefficients (RPC) file, and additional five ground control points (GCPs) for the refinement. As a geometric model, we used the RPC model with zero-order polynomial adjustment39, and orthorectified images were registered to the Universal Traverse Mercator (UTM) Zone 35 N as the reference coordinate system.Georeferencing of scanned cadastral maps and aerial photographsWe used orthorectified WV-3 imageries as a reference for the geometric correction of the historical cadastral maps and the aerial photographs. We selected the spline (triangulation) transformation, a rubber sheeting method, useful for local accuracy and requiring a minimum of 10 control points, as the transformation method to determine the correct map coordinate location for each cell in the historical maps and aerial photographs. The spline transformation provides superior accuracies for the geometric correction of the historical geospatial data40,41.The lack of topographic properties of planimetric features in the historical cadastral maps and the inherent distortions of the aerial photographs due to terrain and camera tilts causes difficulties in precise georeferencing of these data sets. It increases the number of required ground control points (GCPs) for optimal image rectification. Adequate and homogenously distributed GCPs, identified from cadastral maps and aerial photographs, can ensure precise spatial alignment among different geospatial data. The best locations for GCPs were border intersections of agricultural fields and roads. As for artificial objects, places of worship and schools, which are highly probable that have remained unchanged, are other optimal locations for GCPs to perform the accurate geometric correction. Figure 3 displays samples of GCPs selected from cadastral maps and aerial photographs. We obtained 2.11 m or better overall RMSE (Root Mean Square Error) values for the geometric correction of the historical maps and aerial photographs.Figure 3Examples of GCPs selection (red crosses in blue circles) on (a), (c) Cadastral maps and their counterparts on (b), (d) Aerial photographs.Full size imageLULC classification schemeWe defined our classification scheme by analyzing the LULC classes distinguished in all three datasets (i.e., cadastral maps, aerial photographs, and WV-3 imageries). We used the classification scheme shown in Table 1. We also present codes and names of the land cover (LC) classes derived from Corine LC nomenclature42.Table 1 Classification scheme of the study.Full size tableThe legends provided on the historical cadastral maps of Aksu and Kestel delineate 15 LULC categories, which are: (1) buildings, (2) home gardens, (3) roads, (4) arable land, (5) gardens, (6) mulberry groves, (7) chestnut groves, (8) olive groves, (9) vegetable gardens, (10) forest, (11) courtyards, (12) vineyards, (13) arable fields, (14) cemeteries, (15) watercourses. Categorizing the land cover types of cadastral maps is limited with the classes available in the map legend. The legend of cadastral maps categorizes the forested area in one class named “forest”. Therefore, it was not possible to use third-level LC sub-categories in our classification schema for forest area which is separating forested areas into three subclasses (3.1.1, 3.1.2, and 3.1.3) according to the type of tree cover. Although some of the third-level LC sub-categories could be extracted from the cadastral map legend, it was not possible to extract all third level agricultural classes from single-band monochromatic aerial photographs. Although the spatial extent of fruit trees as a permanent crop could be determined from aerial photographs, it was not possible to classify these trees into different fruit types (e.g. 2.2.1 Vineyards, 2.2.2 Fruit trees and berry plantations, 2.2.3 Olive groves). Limitation on the number of forest classes is due to the historical cadastral map content; whereas limitation on the number of agricultural classes is mainly offset by the aerial photographs which have only one spectral band and we did not have any field survey or ancillary geographical data that could help the specific identification of fruit trees.Our primary focus is to find out agricultural land abandonment, deforestation/afforestation, urbanization, and rural depopulation within the historical periods. Therefore, most of the second level LULC classes are sufficient for our purpose. LULC changes within the third class level such as the conversion of third level agriculture classes among each other were not aimed to analyze in this research. Our datasets allow us to use Level 3 CORINE classes for the artificial surfaces. These classes are useful to determine residential area implications of rural depopulation or urbanization, one of the focused transformation types for our analysis.We re-organized and categorized the LULC types used in the cadastral maps, with minimum possible manipulation (only for 2.4 and 3.2 LC classes) according to the classification scheme, as shown in Table 2.Table 2 Correspondence between Corine Land Cover and historical cadastral maps nomenclature.Full size tableLULC mappingAfter aligning all geospatial data, we used the georeferenced cadastral maps, aerial photographs, and satellite images for the LULC mapping. We set the spatial extent of the selected regions based on boundaries digitized from the cadastral maps of 1858. Then we detected historical LULC changes within these extents for all geospatial datasets covering 1900 ha and 830 ha of the Aksu and Kestel regions, respectively. Figures 4 and 5 show the selected extents from the historical maps, aerial photographs, and satellite images of the Kestel and Aksu sites, respectively.Figure 4Geospatial dataset for the Kestel study region. (a) 1858 Cadastral map, (b) 1955 aerial photo, and (c) 2020 WV-3 satellite image (finer details shown in the inset images highlighted by Blue boxes).Full size imageFigure 5Geospatial dataset for the Aksu study region. (a) 1858 Cadastral map, (b) 1955 aerial photo, and (c) 2020 WV-3 satellite image (finer details shown in the inset images highlighted by red boxes).Full size imageDigitization of cadastral maps-1858 LULC mapsWe visually interpreted and manually digitized the geographic features on the historical maps and created vector data for each class. The road features in cadastral maps lack topological properties. They also include spatial errors possibly generated in the process of surveying and map production. Therefore, we cross-checked digitized road segments by visual inspection of the road data of the aerial photographs from 1955. We then further modified road polygons to represent accurate road widths. Afterward, we categorized vectorized features of the cadastral maps into the LULC classes defined in Table 1. Finally, we created the vectorized 1858 LULC map. Figure 6 presents the vectorized 1858 cadastral maps of Aksu and Kestel.Figure 6Vectorized cadastral maps of (a) Kestel and (b) Aksu with Red and green lines showing the vector boundaries.Full size imageObject-based image analysis of aerial photographs-1955 LULC mapsAt the second stage of LULC mapping, we performed the segmentation and classification of the aerial photographs using an object-based approach for generating the 1955 LULC map. The object-based image analysis (OBIA) approach in LULC mapping provides advantages over the traditional per-pixel techniques such as higher classification accuracy, depicting more accurate LULC change, and differentiating extra LULC classes33,43,44. We used the eCognition® software (Trimble Germany GmbH, Munich) to implement an object-based image analysis (OBIA). The OBIA approach contains two phases including the segmentation and classification phases that are performed to locate meaningful objects in an image and categorize the created objects, respectively.Multiple ancillary datasets have been used to support different phases of OBIA. The Open Street Map (OSM) vector data, an open-source geospatial dataset (http://www.openstreetmap.org/), has been utilized as ancillary vector data in OBIA to improve the classification of the remotely sensed images. Sertel et al. (2018) used OSM as a thematic layer for road extraction7. Since there are several limitations in extracting the roads from aerial imagery, the OSM road network data could be useful. A majority of unpaved roads in single-band aerial photographs can easily be misclassified as homogeneous areas of arable lands. Precise detection of the roads from monoband aerial photographs without multi-spectral information is difficult. Therefore, we overlaid the OSM road network data with the aerial photographs to extract the revised aerial road vectors through visual interpretation and manual digitization.We segmented the 1955 aerial photographs with the integration of 1858 LULC map produced from cadastral maps. We implemented the multi-resolution segmentation algorithm. In this segmentation method, a parameter called scale determines the size of resulting objects, and the shape and compactness parameters determine the boundaries of objects. The segmentation process of the aerial photographs was performed at multiple stages with various scale, shape, and compactness parameter values. At the initial stage, we segmented the regions according to the 1858 LULC map and we utilized large-scale parameters. The scale parameter was set to 100 and the shape parameter and the compactness were set as 0.7 and 0.3, respectively. At this stage, we focused on interpreting the objects that have not changed between 1858 and 1955. We classified the segments using the thematic layer attribute (LULC classes defined by the cadastral maps) with the highest coverage. Image objects in which the land surface has changed during 1858–1955 period were detected by visual interpretation and unclassified for further segmentation. We followed this approach to reduce the manual effort. We defined unchanged objects between 1858 and 1955 and assigned the same classes of 1858 LULC map to the objects in 1955 aerial photographs. We then segmented the remaining segments, the last time into smaller objects with the scale parameter set as 25, the shape parameter set as 0.2, and the compactness set as 0.8.We classified the remaining unclassified objects through the development of rulesets. An object can be described by several possible features as explanatory variables which are provided by eCognition. In the classification ruleset, different features and parameters can be defined to describe and extract object classes of interest and thresholds for each feature can be defined by the trial-and-error method. We tested sets of variables for the classification of the monoband aerial photographs. Object features such as the mean value of the monoband, texture after Haralick, distance to neighbor objects, shape features (e.g., rectangular fit and asymmetry), and extent features (e.g., area and length/width) were the most useful alternatives. The classification process of the parcels of the aerial photographs with LULC change started with the classification of roads constructed between 1858 and 1955 by utilizing the aerial road map. The watercourse class was the most difficult to classify since shrubs or trees mostly covered the watercourses. These areas were misclassified as forest or agricultural land. Therefore, experts in historical map reading with local geographical information performed the detection and classification of the water course class and interpreted by the cadastral map (1858) and the google map (2020). After roads and watercourses, we classified forest and agricultural lands using the optimal thresholds for the brightness feature. We calculated the thresholds using the single band of the aerial photograph combined with the area and rectangular fit features. The heterogeneous agricultural areas class principally occupied by agriculture with significant areas of natural grass and trees within the same object are separated from the arable lands using the standard deviation of the digital number (DN) values of the aerial photographs. The texture feature helped classify the permanent crops. The brightness, shape, asymmetry, and distance to road class features were the best-performing ones for classifying the remaining artificial surfaces. The manual interpretation was performed for the classification of sub-classes of artificial surface class, including the continuous/discontinuous urban fabric, industrial, commercial, and transport units, mine, dump and construction sites, and artificial, non-agricultural vegetated areas. Since these land use classes contain one or more land cover and land use categories (e.g., artificial non-agriculture land or industrial or commercial units), finding the optimal threshold and exact feature for distinguishing the subclasses of artificial surfaces is difficult. Especially in the case of using the single-band aerial photographs, manual interpretation was required.Object-based image analysis of satellite images-2020 LULC mapsWe segmented WV-3 satellite images using multi-resolution segmentation algorithm and ancillary geographic data. Similar to the aerial road map, the road network of the study region in 2020, named, WV-3 road map, was extracted by overlaying the OSM road data with the WV-3 satellite image. In the segmentation process of the WV-3 image, we used the vector boundaries of the classified aerial photograph (the 1955 LULC map) and the WV-3 road map as ancillary thematic layers. We opted for the same segmentation and classification approach used for the aerial photographs for the WV-3 image.Firstly, we segmented the satellite image into spectrally homogeneous objects using vector data of the 1955 LULC map by applying large-scale parameters. We implemented scale parameter values of 300, 200, 100, and 50 to find the optimal scale to classify objects that have not changed between 1955 and 2020. The best multi-resolution segmentation configuration was the scale of 100 and the shape and compactness parameters of 0.3 and 0.7, respectively. We classified the segments using the thematic layer attribute (LULC classes defined by the aerial maps) with the highest coverage. Segments with LULC change, e.g. the image objects in which the land surface has changed during 1955–2020 period were detected by visual interpretation and unclassified for further segmentation. As a result, we excluded the objects which were remained unchanged during 1955–2020 by assigning the prepared labels which were allocated in the previous step during the classification of 1955 aerial photographs. We then segmented the remaining objects into smaller objects to identify the changed areas in detail. At this step, the scale, shape, and compactness parameters were set as 25, 0.2, and 0.8, respectively.Except for the additional sets of variables utilized to classify the WV-3 images, we applied the rule-set developed for the classification of the aerial photograph for the classification of the remaining objects of 2020 satellite images. The additional sets of variables include the mean of G, B, R, and NIR and two spectral indices, the Normalized Difference Water Index (NDWI), and the Normalized Difference Vegetation Index (NDVI). NDVI was calculated as the normalized difference of reflectance values in the red and NIR bands; whereas , NDWI was determined as the normalized difference of reflectance values of the green and NIR bands. Through the logical conditions, objects having specified values of NDVI and NDWI can be assigned to vegetation and water classes, respectively. The use of NDVI facilitated the delineation of terrains covered by vegetation and the NDWI improved the extraction of water bodies due to its ability to separate water and non-water objects. We separated different sub-classes of agricultural areas and forests by using optimal thresholds for NDVI which were defined by a trial and error method. Also we utilized assigning the optimal threshold to NDWI to separate water bodies from other land covers. In addition, the mean blue band layer was useful in classifying the artificial surfaces. We assessed the accuracy of each classification using error matrices (overall, user’s and producer’s accuracies, and Kappa statistics)45,46.Estimating LULC changes and LULC conversionsAfter the production of LULC maps of Aksu and Kestel for 1858, 1955, and 2020, the vector data of the LULC maps were used to quantify the LULC conversions for two different periods which are 1858–1955 and 1955–2020. To compare the LULC maps of study areas between two different dates of each study period, we provided detailed “from-to” LULC change information by calculating the LULC change transition matrix computed using overlay functions in ArcGIS.We overlaid LULC maps of 1858 and 1955 and intersected the vector boundaries of the 1858 and 1955 LULC maps to determine the conversion types of LULC classes (from which class to which class). Similarly, to quantify the LULC changes between 1955 and 2020, we overlaid the 1955 and 2020 LULC maps. Then we created transition matrices and performed statistical analysis utilizing the matrices. Finally, we discussed the main LULC change types and the driving factors of the changes in the selected study areas. More