Philippa Kaur

Bulgarelli, D., Schlaeppi, K., Spaepen, S., Ver Loren van Themaat, E. & Schulze-Lefert, P. Structure and functions of the bacterial microbiota of plants. Annu. Rev. Plant Biol. 64, 807–838 (2013).CAS 
PubMed 

Google Scholar 
Leach, J. E., Triplett, L. R., Argueso, C. T. & Trivedi, P. Communication in the phytobiome. Cell 169, 587–596 (2017).CAS 
PubMed 

Google Scholar 
Vorholt, J. A., Vogel, C., Carlstrom, C. I. & Muller, D. B. Establishing causality: opportunities of synthetic communities for plant microbiome research. Cell Host Microbe 22, 142–155 (2017).CAS 
PubMed 

Google Scholar 
Jiang, Y. et al. Plant cultivars imprint the rhizosphere bacterial community composition and association networks. Soil Biol. Biochem. 109, 145–155 (2017).CAS 

Google Scholar 
Garbeva, P., van Elsas, J. D. & van Veen, J. A. Rhizosphere microbial community and its response to plant species and soil history. Plant Soil 302, 19–32 (2008).CAS 

Google Scholar 
Li, Y. et al. Rhizobacterial communities of five co-occurring desert halophytes. PeerJ 6, e5508 (2018).PubMed 
PubMed Central 

Google Scholar 
Matthews, A., Pierce, S., Hipperson, H. & Raymond, B. Rhizobacterial community assembly patterns vary between crop species. Front. Microbiol. 10, 581 (2019).PubMed 
PubMed Central 

Google Scholar 
Perez-Jaramillo, J. E., Mendes, R. & Raaijmakers, J. M. Impact of plant domestication on rhizosphere microbiome assembly and functions. Plant Mol. Biol. 90, 635–644 (2016).CAS 
PubMed 

Google Scholar 
Xu, J. et al. The structure and function of the global citrus rhizosphere microbiome. Nat. Commun. 9, 4894 (2018).ADS 
PubMed 
PubMed Central 

Google Scholar 
Trivedi, P., Leach, J. E., Tringe, S. G., Sa, T. & Singh, B. K. Plant-microbiome interactions: from community assembly to plant health. Nat. Rev. Microbiol. 18, 607–621 (2020).CAS 
PubMed 

Google Scholar 
Howard, M. M., Munoz, C. A., Kao-Kniffin, J. & Kessler, A. Soil microbiomes from fallow fields have species-specific effects on crop growth and pest resistance. Front. Plant Sci. 11, 1171 (2020).PubMed 
PubMed Central 

Google Scholar 
Yan, Y., Kuramae, E. E., de Hollander, M., Klinkhamer, P. G. & van Veen, J. A. Functional traits dominate the diversity-related selection of bacterial communities in the rhizosphere. ISME J. 11, 56–66 (2017).PubMed 

Google Scholar 
Bakker, P. A., Berendsen, R. L., Doornbos, R. F., Wintermans, P. C. & Pieterse, C. M. The rhizosphere revisited: root microbiomics. Front. Plant Sci. 4, 165 (2013).PubMed 
PubMed Central 

Google Scholar 
Lundberg, D. S. et al. Defining the core Arabidopsis thaliana root microbiome. Nature 488, 86–90 (2012).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Busby, P. E. et al. Research priorities for harnessing plant microbiomes in sustainable agriculture. PLoS Biol. 15, e2001793 (2017).PubMed 
PubMed Central 

Google Scholar 
de Vries, F. T., Griffiths, R. I., Knight, C. G., Nicolitch, O. & Williams, A. Harnessing rhizosphere microbiomes for drought-resilient crop production. Science 368, 270–274 (2020).ADS 
PubMed 

Google Scholar 
Hamonts, K. et al. Field study reveals core plant microbiota and relative importance of their drivers. Environ. Microbiol. 20, 124–140 (2018).CAS 
PubMed 

Google Scholar 
Xu, Q. et al. Long-term chemical-only fertilization induces a diversity decline and deep selection on the soil bacteria. mSystems 5, e00337–20 (2020).CAS 
PubMed 
PubMed Central 

Google Scholar 
Richter, A., Schöning, I., Kahl, T., Bauhus, J. & Ruess, L. Regional environmental conditions shape microbial community structure stronger than local forest management intensity. Ecol. Manag. 409, 250–259 (2018).
Google Scholar 
Bais, H. P., Weir, T. L., Perry, L. G., Gilroy, S. & Vivanco, J. M. The role of root exudates in rhizosphere interactions with plants and other organisms. Annu. Rev. Plant Biol. 57, 233–266 (2006).CAS 
PubMed 

Google Scholar 
Wallenstein, M. D. Managing and manipulating the rhizosphere microbiome for plant health: a systems approach. Rhizosphere 3, 230–232 (2017).
Google Scholar 
Kuzyakov, Y. & Xu, X. Competition between roots and microorganisms for nitrogen: mechanisms and ecological relevance. N. Phytol. 198, 656–669 (2013).CAS 

Google Scholar 
Roller, B. R., Stoddard, S. F. & Schmidt, T. M. Exploiting rRNA operon copy number to investigate bacterial reproductive strategies. Nat. Microbiol. 1, 16160 (2016).CAS 
PubMed 
PubMed Central 

Google Scholar 
Wu, L. et al. Microbial functional trait of rRNA operon copy numbers increases with organic levels in anaerobic digesters. ISME J. 11, 2874–2878 (2017).CAS 
PubMed 
PubMed Central 

Google Scholar 
Nuccio, E. E. et al. Niche differentiation is spatially and temporally regulated in the rhizosphere. ISME J. 14, 999–1014 (2020).CAS 
PubMed 
PubMed Central 

Google Scholar 
Fan, K., Weisenhorn, P., Gilbert, J. A. & Chu, H. Wheat rhizosphere harbors a less complex and more stable microbial co-occurrence pattern than bulk soil. Soil Biol. Biochem. 125, 251–260 (2018).CAS 

Google Scholar 
Fan, K. et al. Rhizosphere-associated bacterial network structure and spatial distribution differ significantly from bulk soil in wheat crop fields. Soil Biol. Biochem. 113, 275–284 (2017).CAS 

Google Scholar 
Peiffer, J. A. et al. Diversity and heritability of the maize rhizosphere microbiome under field conditions. Proc. Natl Acad. Sci. USA 110, 6548–6553 (2013).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Baudoin, E., Benizri, E. & Guckert, A. Impact of artificial root exudates on the bacterial community structure in bulk soil and maize rhizosphere. Soil Biol. Biochem. 35, 1183–1192 (2003).CAS 

Google Scholar 
Kuzyakov, Y. & Razavi, B. S. Rhizosphere size and shape: temporal dynamics and spatial stationarity. Soil Biol. Biochem. 135, 343–360 (2019).CAS 

Google Scholar 
Ren, Y. et al. Functional compensation dominates the assembly of plant rhizospheric bacterial community. Soil Biol. Biochem. 150, 107968 (2020).CAS 

Google Scholar 
Chen, Y. et al. Organic amendments shift the phosphorus-correlated microbial co-occurrence pattern in the peanut rhizosphere network during long-term fertilization regimes. Appl. Soil Ecol. 124, 229–239 (2018).ADS 

Google Scholar 
Atulba, S. L. et al. Evaluation of rice root oxidizing potential using digital image analysis. J. Korean Soc. Appl. Bi 58, 463–471 (2015).CAS 

Google Scholar 
Schmidt, H., Eickhorst, T. & Tippkötter, R. Monitoring of root growth and redox conditions in paddy soil rhizotrons by redox electrodes and image analysis. Plant Soil 341, 221–232 (2011).CAS 

Google Scholar 
Pausch, J., Zhu, B., Kuzyakov, Y. & Cheng, W. Plant inter-species effects on rhizosphere priming of soil organic matter decomposition. Soil Biol. Biochem. 57, 91–99 (2013).CAS 

Google Scholar 
Finn, D., Kopittke, P. M., Dennis, P. G. & Dalal, R. C. Microbial energy and matter transformation in agricultural soils. Soil Biol. Biochem. 111, 176–192 (2017).CAS 

Google Scholar 
Jones, R. T. et al. A comprehensive survey of soil acidobacterial diversity using pyrosequencing and clone library analyses. ISME J. 3, 442–453 (2009).CAS 
PubMed 

Google Scholar 
Zhao, S. et al. Biogeographical distribution of bacterial communities in saline agricultural soil. Geoderma 361, 114095 (2020).ADS 
CAS 

Google Scholar 
Eiler, A., Heinrich, F. & Bertilsson, S. Coherent dynamics and association networks among lake bacterioplankton taxa. ISME J. 6, 330–342 (2012).CAS 
PubMed 

Google Scholar 
Zhou, J. et al. Generation of arbitrary two-point correlated directed networks with given modularity. Phys. Lett. A 374, 3129–3135 (2010).ADS 
CAS 
MATH 

Google Scholar 
Herron, P. M., Gage, D. J., Arango Pinedo, C., Haider, Z. K. & Cardon, Z. G. Better to light a candle than curse the darkness: illuminating spatial localization and temporal dynamics of rapid microbial growth in the rhizosphere. Front. Plant Sci. 4, 323 (2013).PubMed 
PubMed Central 

Google Scholar 
Blagodatskaya, E., Blagodatsky, S., Anderson, T. H. & Kuzyakov, Y. Microbial growth and carbon use efficiency in the rhizosphere and root-free soil. PLoS ONE 9, e93282 (2014).ADS 
PubMed 
PubMed Central 

Google Scholar 
Berendsen, R. L., Pieterse, C. M. & Bakker, P. A. The rhizosphere microbiome and plant health. Trends Plant Sci. 17, 478–486 (2012).CAS 
PubMed 

Google Scholar 
Mendes, L. W., Kuramae, E. E., Navarrete, A. A., van Veen, J. A. & Tsai, S. M. Taxonomical and functional microbial community selection in soybean rhizosphere. ISME J. 8, 1577–1587 (2014).CAS 
PubMed 
PubMed Central 

Google Scholar 
Hinsinger, P. Bioavailability of soil inorganic P in the rhizosphere as affected by root-induced chemical changes: a review. Plant Soil 237, 173–195 (2001).CAS 

Google Scholar 
Kuzyakov, Y. & Blagodatskaya, E. Microbial hotspots and hot moments in soil: Concept & review. Soil Biol. Biochem. 83, 184–199 (2015).CAS 

Google Scholar 
Loeppmann, S., Blagodatskaya, E., Pausch, J. & Kuzyakov, Y. Substrate quality affects kinetics and catalytic efficiency of exo-enzymes in rhizosphere and detritusphere. Soil Biol. Biochem. 92, 111–118 (2016).CAS 

Google Scholar 
Ma, X. et al. Spatial patterns of enzyme activities in the rhizosphere: Effects of root hairs and root radius. Soil Biol. Biochem. 118, 69–78 (2018).CAS 

Google Scholar 
Kroener, E., Zarebanadkouki, M., Kaestner, A. & Carminati, A. Nonequilibrium water dynamics in the rhizosphere: How mucilage affects water flow in soils. Water Resour. Res. 50, 6479–6495 (2014).ADS 

Google Scholar 
Carminati, A. Rhizosphere wettability decreases with root age: a problem or a strategy to increase water uptake of young roots? Front. Plant Sci. 4, 298 (2013).PubMed 
PubMed Central 

Google Scholar 
Holz, M., Zarebanadkouki, M., Kaestner, A., Kuzyakov, Y. & Carminati, A. Rhizodeposition under drought is controlled by root growth rate and rhizosphere water content. Plant Soil 423, 429–442 (2018).CAS 

Google Scholar 
Tripathi, B. M. et al. Trends in taxonomic and functional composition of soil microbiome along a precipitation gradient in Israel. Microb. Ecol. 74, 168–176 (2017).PubMed 

Google Scholar 
Harms, A., Brodersen, D. E., Mitarai, N. & Gerdes, K. Toxins, targets, and triggers: an overview of toxin-antitoxin biology. Mol. Cell 70, 768–784 (2018).CAS 
PubMed 

Google Scholar 
Kearns, P. J. & Shade, A. Trait-based patterns of microbial dynamics in dormancy potential and heterotrophic strategy: case studies of resource-based and post-press succession. ISME J. 12, 2575–2581 (2018).CAS 
PubMed 
PubMed Central 

Google Scholar 
Klappenbach, J. A., Dunbar, J. M. & Schmidt, T. M. rRNA operon copy number reflects ecological strategies of bacteria. Appl. Environ. Microb. 66, 1328–1333 (2000).ADS 
CAS 

Google Scholar 
Schoeps, R. et al. Land-use intensity rather than plant functional identity shapes bacterial and fungal rhizosphere communities. Front. Micro. 9, 2711 (2018).
Google Scholar 
Nemergut, D. R. et al. Decreases in average bacterial community rRNA operon copy number during succession. ISME J. 10, 1147–1156 (2016).CAS 
PubMed 

Google Scholar 
Cui, J. et al. Carbon and nitrogen recycling from microbial necromass to cope with C:N stoichiometric imbalance by priming. Soil Biol. Biochem. 142, 107720 (2020).CAS 

Google Scholar 
Blagodatskaya, E. V., Blagodatsky, S. A., Anderson, T. H. & Kuzyakov, Y. Priming effects in chernozem induced by glucose and N in relation to microbial growth strategies. Appl. Soil Ecol. 37, 95–105 (2007).
Google Scholar 
Lecomte, S. M. et al. Diversifying anaerobic respiration strategies to compete in the rhizosphere. Front. Environ. Sci. 6, 139 (2018).
Google Scholar 
Herz, K. et al. Drivers of intraspecific trait variation of grass and forb species in German meadows and pastures. J. Veg. Sci. 28, 705–716 (2017).
Google Scholar 
Ravenek, J. M. et al. Linking root traits and competitive success in grassland species. Plant Soil 407, 39–53 (2016).CAS 

Google Scholar 
Larsen, J., Jaramillo-López, P., Nájera-Rincon, M. & González-Esquivel, C. Biotic interactions in the rhizosphere in relation to plant and soil nutrient dynamics. J. Soil Sci. Plant Nutr. 15, 449–463 (2015).
Google Scholar 
Raaijmakers, J. M., Paulitz, T. C., Steinberg, C., Alabouvette, C. & Moënne-Loccoz, Y. The rhizosphere: a playground and battlefield for soilborne pathogens and beneficial microorganisms. Plant Soil 321, 341–361 (2009).CAS 

Google Scholar 
Ma, H.-K. et al. Steering root microbiomes of a commercial horticultural crop with plant-soil feedbacks. Appl. Soil Ecol. 150, 103468 (2020).
Google Scholar 
Hannula, S. E. et al. Persistence of plant-mediated microbial soil legacy effects in soil and inside roots. Nat. Commun 12, 5686 (2021).ADS 
CAS 
PubMed 
PubMed Central 

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 

Google Scholar 
Hill, T. C., Walsh, K. A., Harris, J. A. & Moffett, B. F. Using ecological diversity measures with bacterial communities. FEMS Microbiol. Ecol. 43, 1–11 (2003).CAS 
PubMed 

Google Scholar 
Lima-Mendez, G. et al. Determinants of community structure in the global plankton interactome. Science 348, 6237 (2015).
Google Scholar 
Noble, W. S. How does multiple testing correction work? Nat. Biotechnol. 27, 1135–1137 (2009).CAS 
PubMed 
PubMed Central 

Google Scholar 
Luo, F., Zhong, J., Yang, Y., Scheuermann, R. H. & Zhou, J. Application of random matrix theory to biological networks. Phys. Lett. A 357, 420–423 (2006).ADS 
CAS 

Google Scholar 
Shannon, P. et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 13, 2498–2504 (2003).CAS 
PubMed 
PubMed Central 

Google Scholar 
Bastian, M., Heymann, S. & Jacomy, M. Gephi: an open source software for exploring and manipulating networks. ICWSM 8, 361–362 (2009).
Google Scholar 
Peng, G. S. & Wu, J. Optimal network topology for structural robustness based on natural connectivity. Phys. A 443, 212–220 (2016).MathSciNet 

Google Scholar 
Ruan, Y., Wang, T., Guo, S., Ling, N. & Shen, Q. Plant grafting shapes complexity and co-occurrence of rhizobacterial assemblages. Microb. Ecol. 80, 643–655 (2020).CAS 
PubMed 

Google Scholar 
Newman, M. E. Modularity and community structure in networks. Proc. Natl Acad. Sci. USA 103, 8577–8582 (2006).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Deng, Y. et al. Molecular ecological network analyses. BMC Bioinforma. 13, 113 (2012).
Google Scholar 
Guimerà, R. & Nunes Amaral, L. A. Functional cartography of complex metabolic networks. Nature 433, 895–900 (2005).ADS 
PubMed 
PubMed Central 

Google Scholar 
Olesen, J. M., Bascompte, J., Dupont, Y. L. & Jordano, P. The modularity of pollination networks. Proc. Natl Acad. Sci. USA 104, 19891 (2007).ADS 
CAS 
PubMed 
PubMed Central 
MATH 

Google Scholar 
Ling, N. et al. Insight into how organic amendments can shape the soil microbiome in long-term field experiments as revealed by network analysis. Soil Biol. Biochem. 99, 137–149 (2016).CAS 

Google Scholar 
Louca, S., Parfrey Laura, W. & Doebeli, M. Decoupling function and taxonomy in the global ocean microbiome. Science 353, 1272–1277 (2016).ADS 
CAS 
PubMed 

Google Scholar 
Douglas, G. M. et al. PICRUSt2 for prediction of metagenome functions. Nat. Biotechnol. 38, 685–688 (2020).CAS 
PubMed 
PubMed Central 

Google Scholar 
Lennon, J. T. & Jones, S. E. Microbial seed banks: the ecological and evolutionary implications of dormancy. Nat. Rev. Microbiol. 9, 119–130 (2011).CAS 
PubMed 

Google Scholar 
Hedges, L. V., Gurevitch, J. & Curtis, P. S. The meta-analysis of response ratios in experimental ecology. Ecology 80, 1150–1156 (1999).
Google Scholar 
Rosenberg, M. S., Adams, D. C. & Gurevitch, J. MetaWin: Statistical software for meta-analysis. Version 2.0. Sinauer (2000).Viechtbauer, W. Conducting meta-analyses in R with the metafor Package. J. Stat. Softw. 36, 1–48 (2010).
Google Scholar 
Egger, M., Smith, G. D., Schneider, M. & Minder, C. Bias in meta-analysis detected by a simple, graphical test. BMJ 315, 629 (1997).CAS 
PubMed 
PubMed Central 

Google Scholar 
Calcagno, V. & de Mazancourt, C. glmulti: an R package for easy automated model selection with (generalized) linear models. J. Stat. Softw. 34, 1–29 (2010).
Google Scholar 
Gloor, G. B., Macklaim, J. M., Pawlowsky-Glahn, V. & Egozcue, J. J. Microbiome datasets are compositional: and this is not optional. Front. Microbiol. 8, 2224 (2017).PubMed 
PubMed Central 

Google Scholar 
Edgar, R. C. MUSCLE: multiple sequence alignment with high accuracy and high throughput. Nucleic Acids Res. 32, 1792–1797 (2004).CAS 
PubMed 
PubMed Central 

Google Scholar 
Letunic, I. & Bork, P. Interactive tree of life (iTOL) v3: an online tool for the display and annotation of phylogenetic and other trees. Nucleic Acids Res. 44, W242–W245 (2016).CAS 
PubMed 
PubMed Central 

Google Scholar  More

Ant experimentsAnt coloniesSeven colonies of the garden ant L. niger collected from the Soka University and a nearby park were used in this study (Extended Data Table S1). They were placed in plastic cases (35 × 25 × 6 cm). Water was provided ad libitum. They were fed a sucrose solution, and were starved for 2–5 days before the start of the experiment. The colonies were queen-less colonies with 200–700 workers. Aqueous sucrose solution was used as a food resource (bait) in the experiments. The laboratory room where the experiments were performed and the ant colonies were kept was maintained at a temperature of 25–27 °C and a humidity of 60–70%. Artificial lights were also installed in this room.ApparatusWe used an apparatus, called “the main apparatus,” with two paths from the nest to the feeding site (length: 30 cm, width: 2 cm, height: 12 cm for the outward path and 15 cm for the return path) (Fig. 1). This apparatus could separate the outward path (bridge) from the inward path (bridge). Here, the outward path refers to that taken by ants from the nest to the feeding site, whereas the inward path refers to the path taken by ants from the feeding site to the nest.Figure 1The main apparatus used in the three experiments (the main experiment and the comparison experiments 1 and 2). Nests are connected to the experimental apparatus by a slope. In the main experiment, on the outward path, there is ant traffic from the nest to the feeding site on a pheromone trail, and on the inward path, there is ant traffic from the feeding site to the nest on a pheromone trail. In the comparison experiment 1, only a pheromone trail is present on both the outward and inward paths. In the comparison experiment 2, no pheromone trail or ant traffic is present on both the outward and inward paths.Full size imageTwo important features of this apparatus were as follows: firstly, it allowed ants to only enter the outward path from the nest. A rat-guard structure at the end of the inward path prevented the ants on the outward path from entering the inward path (Extended Data Fig. S1A). Secondly, we installed a vertical structure at the end of the outward path (height: 4 cm). After climbing the vertical structure, ants were not allowed to return to the outward path (Extended Data Fig. S1B). Moreover, we installed partitions on the feeding site, which also prevented ants from returning to the outward path after reaching the feeding site (Extended Data Fig. S1C). After entering the feeding site, ants had to pass through a narrow gap (width: 0.5 cm) created by the partition. No visual cues were offered as the apparatus was surrounded on all four sides by plastic walls.In this experiment, we made another apparatus for a single ant (target ant), which would be joining the ant trail on the main bridges (Extended Data Fig. S2). This apparatus, called “the confluence device,” was a detachable device that could be connected at right angles to the outward and inward bridges of the main apparatus. To connect this device to the outward bridge, we made the confluence path of this device under the inward bridge of the main apparatus, since the outward bridge was lower than the inward bridge. Thus, we made a slope on the outward confluence path connected to the outward bridge of the main apparatus. Further, because placing the ants directly on the sidewalk sometimes caused them to fall off the sidewalk owing to panic, we constructed a free space and a wall (height: 5 cm) in the middle of the confluence device on which the ants were placed calmly. Owing to this modification, we could let each target ant calm down and then access the main bridge whenever they wanted to. The apparatus used in this experiment was made of white plastic plates.Pheromone trail with ant trafficThis main experiment was limited to once a day for each colony. A sucrose solution was dripped into the feeding site. Target ants, which were walking on a plastic case as foragers, had been moved from their nests to another case immediately before a trail of (nontarget) ants was formed. Thus, dozens of ants were moved in advance to the case to be used as target ants. Subsequently, a trail of (nontarget) ants was formed from the nest to the main apparatus. Considering that it took some time for the ants that had finished foraging and returned to the nest to recruit their mates, the ants were left for approximately 40 min to an hour until a permanent ant trail was formed. It was difficult to form an ant trail immediately after the start of the experiment since no foraging pheromones could be produced in the first foraging trip on the outward path and since experienced foraging ants may make foraging pheromones on the outward path2,21,22. The target ants were allowed to enter bridges of the main apparatus after the establishment of a permanent ant trail. At that time, trails of individual target ants were started. Target ants were allowed to join at right angles to the path on the apparatus, one by one from the confluence device. Individual target ants were allowed to enter the main apparatus at four different points: (1) Left-Left (LL), located at the left side of the center of the outward path. The outward path was on the left side, whereas the inward path was on the right side for the experimenter when seen from the nest. (2) Left–Right (LR), located at the right side of the center of the outward path. (3) and (4) Right-Left (RL) and Right-Right (RR), located at the left and right sides of the center of the inward path, respectively (Fig. 2). We had set these four points to check if target ants tended to turn their body to a certain direction when entering the main bridges, regardless of the movement direction of the other ants. A video camera (Panasonic, AVCHD 30fps) was used to record the migration of ants to the feeding site or nest. Videos were taken from above, and target ants were used only once.Figure 2Four joining points (LL, LR, RL, and RR) and the confluence device (joining device). The confluence (joining) device was connected at right angles to the center of the outward and inward bridges of the main apparatus. Here, the LR version is shown as an example.Full size imageThe goal lines were set at 15 cm from the center of the main paths. We checked the side (nest side or feeding site side) from which a target ant passed the goal line.Pheromone trail with no ant trafficThis comparison experiment 1 was limited to once a day for each colony. Dozens of ants were moved in advance to another case to be used as target ants in a similar manner to the main experiment. The (nontarget) ants were left for about 40 min to an hour until a permanent ant trail was formed. Subsequently, we removed all the ants from the device. Then, target ants were allowed to enter on the side path one by one. In this case, we left the bait in place to control this experiment under the same condition as the main experiment. As the pheromone trail was created on the outward path as well as on the inward path, it was the only decision-making factor for the ants to join at the main path (outward/inward paths). We checked the side (nest side or feeding site side) from which a target ant passed the goal line in a similar manner to the main experiment.No pheromone trail or ant trafficThis comparison experiment 2 was limited to once a day for each colony. Dozens of ants were moved in advance to another case to be used as target ants in a similar manner to the main experiment. This experiment was conducted to investigate ant behavior under the following two conditions: (1) no ant trails and (2) no pheromones trails. The bait was in place in the same manner. We checked the side (nest side or feeding site side) from which a target ant passed the goal line in a similar manner to the main experiment. After each trial (the target ant passed the goal line), we wiped the apparatus with ethanol solution before the next target ant was allowed to enter the main paths.AnalysisThe goal lines were set at 15 cm from the center of the main paths. We checked which goal side the target ants reached the goal line on each trial. A reverse run referred to the goal to the nest on the outward path and the goal to the feeding site on the inward path. A normal run referred to the goal to the feeding site on the outward path and the goal to the nest on the inward path.In some cases of the main experiment, foraging (nontarget) ants that could not reach the feeding site on their outward path or could not return to the nest on their inward path would be against the ant flows. On the outward path, we considered that the ants conducted a “reverse flow” if the position of their heads was on the nest side compared with the position of their stomach. If not, we defined that the ants conducted a “normal flow” (Extended Data Fig. S3). On the inward path, we defined that the ants conducted a “reverse flow” if the position of their head was on the feeding site side compared with the position of their stomach. If not, we defined that the ants conducted a “normal flow” (Extended Data Fig. S3). We focused on the target ants that came in contact with ants with normal flow. Therefore, if an ant with reverse flow was located within 10 cm of the target ant, that trial was excluded from the analysis.Furthermore, we also evaluated if target ants coming in contact with foraging (nontarget) ants immediately after entering the trail would affect the goal choice. Therefore, we conducted an analysis focusing on the contact using the data from the main experiment. We examined whether or not the target ant made contact with other foraging ants until it passed a point 2 cm from the center of the path. As already mentioned, if the target ant came in contact with another ant moving against the normal flow of the ant trail, this contact was excluded from the counts. Moreover, we also excluded cases in which the body of target ants was on a point 2 cm from the center of the path by visual evaluation. Thus, we examined the goal choice of target ants by focusing on whether or not they came in contact with other ants immediately after joining the main bridges.We also conducted a preliminary experiment using a single path apparatus to investigate bi-directional trail behaviour. Please see the Extended Data File S1.Model descriptionThe models were coded using the C programming language. The model description follows the Overview, Design concepts, and Details protocol23,24.PurposeThe purpose of the model was to examine the mechanistic understanding of our findings. We adopted an action of target agents obtained from our ant experiments and compared it with another action of target agents on a trail that was contrary to the fact. To be more precise, target agents were allowed to obey an alignment rule in which they tended to move in the same direction with other agents. We named the former model as the reverse-rule model and the latter model as the alignment-rule model. By doing so, we could find the significance of our findings from ant experiments.Entities, state variables, and scalesWe developed two different models (reverse-rule model and alignment-rule model) that included two types of entities: agents and cells. The agent has the state variable Navigational state, which has two values: Navigational state = {wandering, foraging}. The cell has the state variable Pheromone; this value represents the amount of pheromones in each cell. We used a 2D lattice field and set a straight bridge with 61 cells × 5 cell sizes. We also set goal lines at x-coordinate =  − 30 and 30. If the agents reached coordinates satisfying their x-coordinate =  − 30 or 30, they were removed from the system. If the agents reached y-axis boundaries, their movement direction was restricted. Each trial continued until the target agent reached one of the two goal lines. However, trials were forcibly finished if the target agent never reached any goal line by t = 500-time steps. In total, we conducted 1000 trials.Process overview and schedulingAt the beginning of each trial, an artificial target ant (Navigational state = wandering) was introduced at the center of an artificial simulation field. Foraging agents (Navigational state = foraging) were randomly distributed on the simulation field in advance.Agents on the simulation field selected one direction from two directions (+ x and − x) on each time step and updated their positions. Briefly, an agent at coordinate (x, y) selected one direction from two directions (+ x and − x) and updated its position with one of the three coordinates—(x − 1, y), (x − 1, y + 1), or (x − 1, y − 1)—if it selected the − x direction, or—(x + 1, y), (x + 1, y + 1), or (x + 1, y − 1)—if it selected the + x direction by scanning pheromones on these three coordinates. For example, if an agent at coordinate (0, 2) decided to move in + x direction at one time, the position of this agent was replaced with one of (1, 3), (1, 2) and (1, 1) from (0, 2) by scanning pheromones on these three coordinates. The target agent selected the − x/ + x direction with equal probability on each time step until it met the foragers. In contrast, foraging agents tended to decide to move in the − x direction on each time step with a high probability and therefore they tended to select the − x direction for position updating. Foraging agents deposited pheromones before leaving the current cell (see submodel entitled “Position updating” and submodel entitled “Pheromone updating”). In contrast, the target agents did not deposit pheromones.Using above submodels, artificial ants sometimes met other agents. If the target agent (Navigational state = wandering) met the foragers (Navigational state = foraging), the target agent tended to select one direction from two directions (+ x and − x) on each time step thereafter with a high probability, which was dependent on which direction the met foragers came from. More strictly, in the reverse-rule model, the target agent tended to move in an opposite direction from the foragers if it met the foragers coming from the opposite direction. On the contrary, the target agent in the alignment-rule model tended to move in the same direction with foragers if it met the foragers moving in the same direction (see submodel entitled “The interaction between the target agent and foragers”). For example, in the reverse-rule model, if the target agent at coordinate (x, y), whose previous coordinate was (x − 1, y), met the forager coming from the opposite direction, whose previous coordinate was (x + 1, y), the target agent decided to move in + x direction on each time step thereafter with a high probability until similar events occurred.Design conceptThe mean goal time was the emergent property of the model. Sensing was important as the agents scanned the pheromone concentrations. Stochasticity was used to determine in which direction the agent moved and to select one cell using the pheromone concentrations.InitializationWe set a single agent (target agent) on the coordinate (0, 2) and its Navigational state was set to wandering (Extended Data Fig. S5A). We also set N foraging agents on the bridge whose Navigational state was set to foraging. Therefore, N + 1 agents were on the test field at the beginning of each trial. A target agent was the agent k = 0, whereas foraging agents were agents k = 1, 2, …, N. These foragers were randomly distributed on the bridge. Thus, x(k) (in) {n |− 30 ≤ n ≤ 30, n is an integer} and y(k) (in) {n | 0 ≤ n ≤ 4, n is an integer} for k  > 0.Foraging agents were set to move in the -x direction (Direction(k) for k  > 0 = − x). On the other hand, the target agent randomly chose one direction from two directions at the beginning of each trial (Direction(0) was set to + x or − x with equal probability). Herein, Direction(k) can be − x or + x, which implies bias in the movement direction. The parameter prob(k) indicates the probability of moving in Direction(k). The target agent selected the − x/+ x direction with equal probability on each time step until it met the foragers. Therefore, the parameter prob was set to 0.50 for the target ant (prob(0) = 0.50), whereas prob was set to 0.80 for foraging agents (prob(k) = 0.80 for k  > 0). The amount of pheromones on each cell was set to 1 at the beginning of each trial (pheromone(x, y) = 1) and the pheromone evaporation rate q was set to 0.99.The model descriptions are explained using submodels. A Submodel: the interaction between the target agent and foragers causes differences between two rules (the reverse-rule model and the alignment-rule model).SubmodelsSubmodel: the interaction between the target agent and foragersThe parameters Direction(0) and prob(0) were replaced with new ones whenever the following events occurred.In the reverse-rule model, for any agent k (k  > 0),Herein, (xt(k), yt(k)) indicates the x–y-coordinate for the agent k at time t. Furthermore, (xt(0), yt(0)) = (xt(k), yt(k)) means that the target agent and the agent k occupy the same cell at time t while (xt(0) − xt−1(0)) × (xt(k) − xt−1(k)) =  − 1 indicates that the target agent meets the agent k came from the opposite direction. The target agent replaces Direction(0) with an opposite direction from the forager k (see Extended data Fig. S5B).In the alignment-rule model, for any agent k (k  > 0),(xt(0) − xt−1(0)) × (xt(k) − xt−1(k)) = 1 indicates that the target agent meets the agent k came from the same direction. The target agent replaces Direction(0) with a same direction with the forager k (See Extended Data Fig. S5B).In the reverse-rule model, these events are driven from the experimental observations of real ants. Target ants appear to move against the trail and seem to move straight by contacting those other nestmates that come from the opposite direction. Also, target ants seem to select the reverse goal even if physical contact with ant nestmates does not occur immediately after entering the bridge. So, regarding parameter replacements, we did not consider the position at which the target agent met another agent. Note that foraging agents did not change these parameters until the end of each trial. Further, Direction(0) can be replaced with − x from + x and vice versa whenever the target agent meets foragers that come from the opposite direction.In the alignment-rule model, the target agent tends to move in the same direction with other agents. This is contrary to the experimental observations of real ants.Submodel: position updatingFor all k agents (k = 0–N), the movement direction and position updates are shown as follows (Extended Data Fig. S5C);Here, rnt(k) indicates a random number. Thus, rnt(k) (in) [0.00, 1.00].Prob(0) for the target agent is initially set to 0.50. Therefore, the target agent selects one direction from the two (− x and + x) on each time step randomly before the condition described in submodels—the interaction between the target agent and foragers is satisfied. On the other hand, foraging agents select − x direction with a high probability (= Prob(k)) on each time step. After selecting one direction from two (− x and + x), agents scan three cells in the direction of movement. Using pheromone concentrations on those three cells, they update their positions.If agents reach coordinates satisfying their y-coordinate = 4 or 0, those agents update their position by selecting not three but two coordinates since they are located on the edges of the bridge.Submodel: pheromone updatingForaging agents (k  > 0) deposited pheromones on the current cell when leaving that cell.Then, pheromones are evaporated using the evaporation rate q.For each time iteration, these submodels operated in the following order.STEP 1: The interaction between the target agent and foragers.STEP 2: Position updating.STEP 3: Pheromone updating.AnalysisTo check the accuracy of our model, we counted which goal side the target agent entered the goal line from using the reverse-rule model by setting N = 9. If the target agent passed the goal line at x-coordinate =  − 30 (30), we considered that it reached the normal (reverse) goal. Note that trials in which the target agent never reached any goal lines by t = 500 were excluded from this analysis. Furthermore, to investigate the adaptability of the reverse run mechanism, we examined the time until the target agent reached the goal lines using the reverse-rule model and the alignment-rule model. Herein, we set two different conditions with respect to the number of foraging agents (N = 4 and 9). More

Study subjects and measurement of skin physiological parametersWe analyzed skin microbiome and mycobiome from cheeks and foreheads of healthy younger (19–28 years old, Y-group) and older (60–63 years old, O-group) Korean women who were free from cutaneous disorders (Table 1 and Supplementary Table S1). All 61 subjects had been living in Seoul, Korea, for more than 3 years with normal skin conditions. We preferentially selected those who had sebum secretion greater than 30 arbitrary units and moisture greater than 50 arbitrary units in both groups. Among the measurements of moisture content, pH, sebum content, and transepidermal water loss (TEWL), only sebum and TEWL decreased significantly in the O-group compared to the Y-group in the cheeks (P = 2.25e−06, Wilcoxon rank-sum test; P = 0.019, Welch two-sample t test) and forehead (P = 1.33e−06, Wilcoxon rank-sum test; P = 0.003, Welch two-sample t test). Whereas no significant differences were found in the average values for moisture (cheeks: Y-group, 59.9; O-group, 56.6; forehead: Y-group, 61.1; O-group, 58.7) and pH (cheeks: Y-group, 6.0; O-group, 5.8; forehead: Y-group, 6.0; O-group, 5.6) between the two age groups.Table 1 Characteristics of subjects for aged related skin microbiome and mycobiome study.Full size tableComparisons in cheek and forehead microbiome and mycobiome between the two age groupsWe analyzed bacterial communities from 27 Y-group samples (cheeks, n = 13; forehead, n = 14) and 24 O-group samples (cheeks, n = 12; forehead, n = 12) and fungal communities from 28 Y-group samples (cheeks, n = 15; forehead, n = 13) and 32 O-group samples (cheeks, n = 16; forehead, n = 16), except for samples that were eliminated from the Illumina Mi-Seq sequencing due to low sequence reads (bacteria,  3. 0) (Fig. 4). Pathways belonging to the metabolism category were dominant in each age group. In the cheek of the Y-group, pathways involved in energy metabolism by bacteria, such as glycolysis/gluconeogenesis, citrate cycle, pentose phosphate pathway, fructose and mannose metabolism, galactose metabolism, d-alanine metabolism, and thiamine metabolism, were predominant, whereas in the cheek of the O-group, degradation-related pathways, such as fatty acid degradation, synthesis and degradation of ketone bodies, benzoate degradation, and chloroalkane and chloroalkene degradation, were predominant. In the forehead of the Y-group, glycolysis/gluconeogenesis, pentose phosphate pathway, fructose and mannose metabolism, galactose metabolism, d-glutamine and d-glutamate metabolism, d-alanine metabolism, and thiamine metabolism pathway were significantly more abundant, whereas in the forehead of the O-group, fatty acid degradation, synthesis and degradation of ketone bodies, valine/leucine and isoleucine degradation, and limonene/pinene degradation pathway were significantly more abundant.Figure 4Heat map for significantly different predicted functional pathways on (a) cheeks and (b) foreheads of Korean women by age based on LEfSe analysis (LDA score  > 3.0).Full size imageThe metabolism pathway for biotin, a water-soluble vitamin that is effective for skin health and essential for keratin production15, was more prevalent in the cheek and forehead of the Y-group. Interestingly, the metabolism pathway for lipoic acid, which is known to possess beneficial effects against skin aging and is used widely in cosmetic and dermatological products16,17, was significantly higher in the foreheads of the Y-group. We tracked the specific ASVs possessing these pathways, in both biotin metabolism and lipoic acid metabolism, Cutibacterium sp. (ASV2136 and ASV2130) and Staphylococcus sp. (ASV3008) were predicted to have the top three relative abundances in KOs. The relative abundances in biotin metabolism and lipoic acid metabolism of Cutibacterium sp. (ASV2136) were 24.9% and 26.1%, respectively. The relative abundances for each pathway for Staphylococcus sp. (ASV3008) were 10.2% and 18.7%, and for Cutibacterium sp. (ASV2130), they were 9.3% and 10.0%, respectively. We confirmed these two pathways in the genome of skin bacteria, C. acnes (Supplementary Fig. S2). These additional analyses support the reliability of the function in the skin environment of Cutibacterium. Interestingly, from the LEfSe result, Cutibacterium sp. (ASV2136) had a significantly higher abundance in the cheek and forehead microbiome of the Y-group. The pathway of biosynthesis of lipopolysaccharide, also known as bacterial endotoxins, showed higher abundance in the cheek and forehead microbiome of the O-group. The ASVs that contribute to inferring the LPS biosynthesis pathway were identified as Paraburkholderia sp. (ASV5030) and B. vesicularis (ASV4155). Also, pathways related to antibiotic biosynthesis (biosynthesis of vancomycin group antibiotics) and bacterial motility (bacterial chemotaxis and flagellar assembly; both belonging to the cellular processes category) were prominent in the cheek and forehead of the O-group. PICRUSt2 analysis implied that, regardless of skin site differences, the potential functions of the microbial community that compose the skin microbiome were similar according to age.Network analysis on cheek and forehead microbiome and mycobiomeWe performed SParse InversE Covariance estimation for Ecological Association Inference (SPIEC-EASI) analysis to evaluate the overall network of the skin microbes. The results of network density (D) on 81 cheek and 87 forehead ASVs, calculated using the ratio of the number of edges, showed higher network density in the skin microbiome of the Y-group (D = 0.015 and D = 0.001, in cheek and forehead, respectively) than the O-group (D = 0.007 and D = 0.007, respectively) (Fig. 5). To examine network correlation between bacteria and fungi, network density for Bacteria–Fungi (DBF) was calculated by the actual number of edges and a potential number of edges in a correlation ([bacterial nodes × fungal nodes]/2). We confirmed higher network density in the cheek of the Y-group (DBF = 0.008) than the O-group (DBF = 0) and edges of the major bacterial and fungal taxa, such as Staphylococcus sp. (ASV3008)—M. sympodialis (ASV500) and Roseomonas sp. (ASV4088)—M. restricta (ASV482), were observed in the cheek of the Y-group. In the forehead, edges of Methylobacterium sp. (ASV4314)—M. globosa (ASV454), Methylobacterium sp. (ASV4314)—Zygosaccharomyces rouxii (ASV208), and Venionella sp. (ASV3575)—M. sympodialis (ASV500) were observed in the Y-group, and edges of Cutibacterium sp. (ASV2107)—M. globosa (ASV461), Staphylococcus sp. (ASV3024)—M. arunalokei (ASV446), and Methylobacterium (ASV4314)—M. dermatis (ASV448) were observed in the O-group (DBF = 0.004). We found a network between bacteria and fungi with different kingdom levels in the skin microbiome, and especially, we confirmed that different genus or species level microbe was involved in the microbial network according to skin location and Y-, O-group.Figure 5Network analysis of the ASVs on (a) cheeks and (b) forehead of Korean women. Each node represents the ASV and the size of the node is based on relative abundance of each ASV. Color markings indicate the major taxa except for unidentified bacteria or fungi. Shapes represent the level of kingdom, Bacteria (bold) and Fungi (dotted line). The ASVs were selected for bacterial ASVs found in more than half of all samples on the cheeks and forehead, respectively, and for the fungal ASVs with a relative abundance of more than 0.1% in each of the cheeks and forehead samples. The D value is network density calculated using the ratio of the number of edges.Full size image More

Analytical targetsTwo species of undulation motion rays with different pectoral fin shapes bred in KAIYUKAN were analyzed: sharpnose stingray Dasyatis acutirostra and pitted stingray Dasyatis matsubarai (Fig. 1a,b). Blender 2.7925 was used to construct stingray models from pictures26,27 as accurately as possible; Blender is a free and open-source 3D creation suite used to make realistic characters for movies, etc. Detailed information on how to construct models using Blender is provided in our previous paper28. To focus on the effects of pectoral fin movements, we did not consider the body’s shape as in the previous studies12,29. The height and disk width (WD) of all models were set to 0.01 m and 0.44 m, respectively, considering the previous studies30,31. The disk length of each model was determined from WD, referring to the aspect ratio of the rays’ photographs26,27; the disk length (LD) of D. acutirostra and D. matsubarai are 0.348 m and 0.344 m, respectively.Figure 1Analytical targets and description of motion. (a) Analytical model of D. acutirostra. (b) Analytical model of D. matsubarai. (c) Description of motion, (d) the relationship between any two points on the surface before and after the deformation.Full size imageMotionThe motion was given to satisfy the following equations:$$z = left{ {begin{array}{*{20}l} {A;sin left( {omega left( {t – kTleft( {frac{{angl{text{e}}left( {x_{i} ,y_{i} } right) – 10^{{text{o}}} }}{{All;angl{text{e}}}} – theta } right)} right)h_{1} h_{2} } right.} hfill & {left( {10^{{text{o}}} le angl{text{e}}left( {x_{i} ,y_{i} } right) le 170^{{text{o}}} } right)} hfill \ {A;sin left( {omega left( {t – kTleft( {frac{{350^{{text{o}}} – left( {angl{text{e}}left( {x_{i} ,y_{i} } right) – 10^{{text{o}}} } right)}}{{All;angl{text{e}}}}} right)} right)h_{1} h_{2} } right.} hfill & {left( {190^{{text{o}}} le angl{text{e}}left( {x_{i} ,y_{i} } right) le 350^{{text{o}}} } right)} hfill \ end{array} } right.$$
(1)
$$begin{array}{c}{h}_{1}=a{r}_{i}^{3}+b{r}_{i}^{2}+c{r}_{i}end{array}$$
(2)
$$h_{2} = left{ {begin{array}{*{20}l} {dleft( {angleleft( {x_{i} ,y_{i} } right) – 10^{ circ } } right)^{2} + eleft( {angleleft( {x_{i} ,y_{i} } right) – 10^{ circ } } right)} hfill & {left( {10^{ circ } le angleleft( {x_{i} ,y_{i} } right) le 170^{ circ } } right)} hfill \ {dleft( {350^{ circ } – left( {angleleft( {x_{i} ,y_{i} } right) – 10^{ circ } } right)} right)^{2} + eleft( {350^{ circ } – left( {angleleft( {x_{i} ,y_{i} } right) – 10^{ circ } } right)} right)} hfill & {left( {190^{ circ } le angleleft( {x_{i} ,y_{i} } right) le 350^{ circ } } right)} hfill \ end{array} } right.$$
(3)
$$begin{array}{c}{left({r}_{i}-{r}_{i-1}right)}^{2}+{left({z}_{i}-{z}_{i-1}right)}^{2}={left({r}_{i}^{mathrm{^{prime}}}-{r}_{i-1}^{mathrm{^{prime}}}right)}^{2}+{left({z}_{i}^{mathrm{^{prime}}}-{z}_{i-1}^{mathrm{^{prime}}}right)}^{2}end{array}$$
(4)
$$begin{array}{c}angleleft({x}_{i},{y}_{i}right)= angleleft({x}_{i}^{^{prime}},{y}_{i}^{^{prime}}right).end{array}$$
(5)
Equation (1) represents the amount of movement of the model surface in the z-axis direction, where (A) is the amplitude of the pectoral fin tip, (omega) is the angular velocity, (t) is time, (k) is the wavenumber, (T) is the period, angle(({x}_{i},{y}_{i})) is the angle made by the line connecting the center of rotation and any point (({x}_{i},{y}_{i})) on the model surface with the x-axis, Allangle is the range where the motion is given (160°), and (theta) is the phase difference between the movements of the right and left pectoral fins (Fig. 1c). ({h}_{1}) is the weighting from the center to the radial direction: it is necessary to set the amplitude at the ray’s center to zero and increase the amplitude toward the pectoral fin tip (a = 119.786, b = -7.957, c = 0.498). ({h}_{2}) is the weighting in the circumferential direction: it is necessary to increase the amplitude from the anterior to the tip of the pectoral fin and decrease the amplitude from the tip of the pectoral fin to the posterior (d = − 1.563 × 10–4, e = 0.025). Equation (4) is the condition in which the distance between two neighboring points in the same radial direction is equal before and after the movement (Fig. 1d). (r) is the distance between the center of rotation and any point (({x}_{i},{y}_{i})), defined as (sqrt{{x}_{i}^{2}+{y}_{i}^{2}}). Equation (5) defines angle(({x}_{i},{y}_{i})) as being constant before and after the move (Fig. 1c). The variables after the move are marked with ‘. Variables used in the analysis are A = 0.089 m, T = 0.499, k = 1.270, and ω = 12.599 rad/s. Videos of the created motion from the front and the side are shown in the “Supplement” (Supplement Movies 3, 4).Analytical conditionsAnalysis cases were conducted with eight conditions: two types of pectoral fin shape (Fig. 1a,b) and four types of phase difference (0 (T), 0.25 (T), 0.5 (T), and 0.75 (T)). These conditions were set for investigating the effects of phase differences between left and right pectoral fin movements on swimming and how these effects vary with pectoral fin shape.Numerical methodsA CFD simulation of the ray models in the water flow was performed using OPENFOAM, an open-source finite volume method CFD toolbox32, to calculate the forces acting on the rays in each axial direction and the moment around each axis. The governing equations were the continuity equation and the three-dimensional incompressible Reynolds-averaged Navier–Stokes equation, expressed by:$$begin{array}{*{20}c} {nabla cdot u = 0} \ end{array}$$
(6)
$$begin{array}{*{20}c} {frac{partial u}{{partial {text{t}}}} + nabla cdot left( {uu} right) = – nabla p + nabla cdot left( {vnabla u} right) + nabla cdot left[ {nu left{ {left( {nabla u} right)^{T} – frac{1}{3}nabla cdot uI} right}} right], } \ end{array}$$
(7)
where (u) is the velocity vector, t is the time, p is the static pressure divided by the reference density, (nu) is the kinematic viscosity, and I is the unit tensor. The Reynolds number was defined regarding the previous studies3 as:$$begin{array}{c}{R}_{e}=frac{U{L}_{D}}{nu },end{array}$$
(8)
where (U) (/ms) is the given flow speed3 (1.5 × LD/ms), ({L}_{D}) (m) is the length of the ray models, and (nu) is the kinematic viscosity of water at 20 °C (1.0 × 10–6 m2/s). The Reynolds number in this study is 1.8 × 105; considering this, we used the k–ω shear stress turbulence model33,34. The k–ω shear stress turbulence model is a type of Reynolds-averaged Navier–Stokes equation (RANS) turbulence model that is widely used to calculate for the fish swimming flow35,36,37. The overset grid method was used in this study; it is a generic implementation of overset meshes. For both static and dynamic cases, cell-to-cell mapping between multiple, disconnected mesh regions is employed to generate a composite domain38,39. This method permits complex mesh motions and interactions without the penalties associated with deforming meshes. The process is described in detail by Noack40. The calculation volume was 5.4 WD in length, 5.4 WD in height, and 5.4 WD in width (Fig. 2a,b). A hexahedral volume mesh was created using the snappyHexMesh of OPENFOAM. The fluid region was divided into two parts: the overset region and the background region (Fig. 2a,b). The overset region moves and transforms to match the motion of the ray and was made with fine meshes around the analysis target and coarse meshes in the outlying areas; a one-layer boundary layer mesh was created around the analysis target. The overset region shape is an ellipsoid (Fig. 2a,b). The minimum mesh volume is 7.3 × 10–10 (m3), and the maximum mesh volume is 2.6 × 10–2 (m3). The total number of meshes was 9.0 × 105. At the outlet boundary, the average static relative pressure was set to 0 Pa. The surfaces of the fish model were formed into non-slip surfaces.Figure 2Meshes for CFD simulation and differences in force between different meshes. (a) Meshes at the coronal plane of the whole fluid region. (b) Frontal cross-section of the fluid region at the green line in (a). The red region is the overset region. (c,d) Comparison of the instantaneous drag coefficient and the moment coefficient around the y-axis of D. matsubarai between the coarse, fine, and dense mesh.Full size imageThe drag coefficient ({C}_{D}left(tright)), the lateral force coefficient ({C}_{l}left(tright)), the lift coefficient ({C}_{L}left(tright)), the moment coefficient around the x-axis ({C}_{mx}left(tright)), the moment coefficient around the y-axis ({C}_{my}left(tright)) and the moment coefficient around the z-axis ({C}_{mz}left(tright)) were calculated as:$$begin{array}{c}{C}_{D}left(tright)=frac{Dleft(tright)}{frac{1}{2}rho {U}^{2}{L}_{D}{W}_{D}}end{array}$$
(9)
$$begin{array}{c}{C}_{l}left(tright)=frac{lleft(tright)}{frac{1}{2}rho {U}^{2}{L}_{D}{W}_{D}}end{array}$$
(10)
$$begin{array}{c}{C}_{L}left(tright)=frac{Lleft(tright)}{frac{1}{2}rho {U}^{2}{L}_{D}{W}_{D}}end{array}$$
(11)
$$begin{array}{c}{C}_{mx}left(tright)=frac{{M}_{psi }left(tright)}{frac{1}{2}rho {U}^{2}{L}_{D}^{2}{W}_{D}}end{array}$$
(12)
$$begin{array}{c}{C}_{my}left(tright)=frac{{M}_{phi }left(tright)}{frac{1}{2}rho {U}^{2}{L}_{D}^{2}{W}_{D}}end{array}$$
(13)
$$begin{array}{c}{c}_{mz}left(tright)=frac{{M}_{theta }left(tright)}{frac{1}{2}rho {U}^{2}{L}_{D}^{2}{W}_{D}},end{array}$$
(14)
where (Dleft(tright)) is the calculated drag, (lleft(tright)) is the calculated lateral force, (Lleft(tright)) is the calculated lift, ({M}_{psi }left(tright)) is the calculated moment around the x-axis, ({M}_{phi }left(tright)) is the calculated moment around the y-axis, ({M}_{theta }left(tright)) is the calculated moment around the z-axis, and (rho) (kg/m3) is the density of water at 20 °C (998 kg/m3). As shown in a previous study41. the propulsive efficiency (eta) is defined as the ratio of output power ({P}_{o}) to input power ({P}_{e}) which can be written as:$$begin{array}{c}{P}_{o}left(tright)=frac{1}{T}{int }_{0}^{T}Dleft(tright)Udtend{array}$$
(15)
$$begin{array}{c}{P}_{e}left(tright)=frac{1}{T}{int }_{0}^{T}left[Dleft(tright)dot{x}left(tright)+lleft(tright)dot{y}left(tright)+Lleft(tright)dot{z}left(tright)right]dtend{array}$$
(16)
$$begin{array}{c}eta =frac{{P}_{o}}{{P}_{e}}.end{array}$$
(17)
An in-house program calculated the forces acting on rays in each axial direction and the moment around each axis. The numerical method’s validity and reliability were verified by comparing previous experimental and numerical analytical studies of heaving and pitching on airfoil naca001341. A high degree of similarity to previous studies was confirmed; the mean difference in the propulsive efficiency from the previous study of analysis was 5%, and the difference from the previous study of the experiment was 9%. Detailed information such as mesh, length, and velocity, of this analysis method’s verification is provided in the “Supplement”.A grid sensitivity study was conducted using three meshes: coarse, fine, and dense. The coarse mesh has 8.1 × 105 elements, the fine mesh has 9.0 × 105 elements, and the dense mesh has 9.9 × 105 elements. The analysis was conducted using a condition with no phase difference of D. matsubarai. As shown in Fig. 2c,d, the drag coefficient and the moment coefficient around the y-axis are almost the same when the mesh is fine and when the mesh is dense. The mean drag and propulsive efficiency error of fine and coarse meshes are 2.7% and 3.5%, respectively. The fine mesh was used in all simulation cases considering accuracy. We used the same meshes for all cases. More

Thank you for visiting nature.com. You are using a browser version with limited support for CSS. To obtain
the best experience, we recommend you use a more up to date browser (or turn off compatibility mode in
Internet Explorer). In the meantime, to ensure continued support, we are displaying the site without styles
and JavaScript. More

A workflow for the calculation of CSI is presented in Fig. 1c. For all the analyses, we used the R v. 4.0.3 software environment20 implemented in RStudio v. 1.4.1103. The scripts used for each methodological step are available at the Figshare repository21. After data download from primary sources (PaleoClim and WorldClim), specifically for the CSI-future map set we performed an initial step aimed to obtain individual bioclimatic variables for each future time period for the four SSPs (Fig. 1b). To achieve this, the median values of nine GCMs were calculated in functions compiled in raster R package22 for each individual bioclimatic variable (see a few exceptions of number of GCMs used in Table 2).Table 2 General circulation models (GCM) used to construct the future map sets.Full size tableThe standard deviation (SD) was estimated as a measure of the amount of variation or dispersion along time series, from which the resulting output maps showed the places where climate conditions remained constant or variable across the temporal periods considered (Fig. 1a,b). The SD, as a way to identify stable/unstable climatic areas, was previously used in other climatic or evolutionary studies4,14. To compute the SD output rasters, we applied the mosaic function setting “fun = sd” from raster R package, calculating the SD for each pixel in the 12 time period rasters for CSI-past and five times for CSI-future, independently for each variable. The mosaic function was also used for the range calculation, with “fun = min” and “fun = max” to obtain the minimum and maximum values of input rasters, respectively, with a further step for subtracting maximum to minimum values.Specifically, for CSI-past, as it includes several time periods with sea-level dropping below the present level (T1, T3, T5, T6, T7, T8, T9; Fig. 1a), we applied a mask of the current land surface, i.e. taking the T12 (Anthropocene) as a template. With this additional step, we were able to remove those pixels (grid cells) currently under the sea but that were once emerged. Most of these pixels, however, were only emerged during the LGM (ca. 21 ka), thus having values for bioclimatic variables for just a single time period (instead of the 12 routinely used for the variability estimation). The inclusion of these areas would result in highly climatically stable regions (low SD values; Supplementary Fig. 1), but this would be an obviously biased result. In contrast, we did not remove those areas affected by the sea-level rising periods, as only three periods contained “NoData” values (T2, T4, T10; Fig. 1a). However, to take this fact into consideration, we created a raster file in which these areas submerged during warm periods are indicated (see Supplementary Fig. 1). Finally, for both CSI-past and CSI-future, the resulting SD values were normalized to values between 0 and 1, with 0 representing completely stable areas and 1 the most unstable ones.The next step was focused on the selection of a relatively uncorrelated set of variables for each map set. We used the removeCollinearity function from virtualspecies R package23 that estimates the correlation value among pairs of variables from a given number of random sample points (10,000 in present case) according to a given method (Pearson for the present case) and a threshold of statistic selected (r  > 0.8 as a cut-off value). The function removeCollinearity returns a list of uncorrelated variables according to the settings specified, randomly selecting just one variable from groups of correlated ones (see Table 1 for a complete list of variables used for each map set). As we compiled estimates of variability independently for each variable and map set (e.g. SD bio1 past, SD bio2 past, etc.), each user can define his own CSI, selecting the more interesting variables according to the case of study.The final CSI maps were obtained by summing the SD values of the variables selected and the subsequent outputs normalized (0 to 1) (Figs. 2–4). Histogram plots were represented with ggplot2 R package24 and maps were exported with ArcGIS v.10.2.2 (Esri, Redlands, California, USA 2014). The histograms were computed for these final CSI maps, which represent the frequency and distribution of CSI values. We presented the final CSI maps with two different colour ramp schemes with ArcGIS. The first consisted of defining equal interval breaks from 0 to 1. The second was based on defining 32 categories with different value breaks for past and future map sets according to the value frequency shown by the histogram plot, i.e. the category with the highest CSI values (no. 32) was 0.71–1 in the past map set and 0.356–1 in the future map set.Fig. 2Maps of Climate Stability Index (CSI) values for the past map set from Pliocene (3.3 Ma) to present (1979–2013), at 2.5 arc-min grid resolution. Colours range from blue for low standard deviation (SD) values, which represents areas with low climatic fluctuations (i.e, low values of CSI) during the period Pliocene–present, to red for high SD values, which shows areas where high climatic fluctuations would have taken place (i.e., high values of CSI). On the upper map, the colour ramp shows equal interval breaks. The histogram with frequency and distribution of CSI values is also shown. On the lower map, the colour ramp has been manually adjusted to a defined set of break values (see details in the text).Full size imageFig. 3Maps of Climate Stability Index (CSI) values for the future conditions (Shared Socioeconomic Pathways: SSP1-2.6, SSP2-4.5, SSP3-7.0, and SSP5-8.5) from present (1970–2000) to future (2100), at 2.5 arc-min grid resolution. Colours range from blue for low standard deviation (SD) values, which represents areas with low climatic fluctuations (i.e, low values of CSI) from present to future, to red for high SD values, which shows areas where high climatic fluctuations would have taken place (i.e., high values of CSI). The colour ramp shows equal interval breaks. The histogram with frequency and distribution of CSI values is also shown for each future scenario.Full size imageFig. 4Maps of Climate Stability Index (CSI) values for the future conditions (Shared Socioeconomic Pathways: SSP1-2.6, SSP2-4.5, SSP3-7.0, and SSP5-8.5) from present (1970–2000) to future (2100), at 2.5 arc-min grid resolution. Colours range from blue for low standard deviation (SD) values, which represents areas with low climatic fluctuations (i.e, low values of CSI) from present to future, to red for high SD values, which shows areas where high climatic fluctuations would have taken place (i.e., high values of CSI). The colour ramp has been manually adjusted to a defined set of break values (see details in the text).Full size image More

Study siteAn open field small-scale release of a GM strain of Anopheles mosquitoes was carried-out in July 2019 in the village of Bana in Western Burkina Faso (see Supplementary Fig. 5). The study was granted regulatory authorisation from the National Biosafety Agency (NBA) (order No. 2018-453/MESRSI/SG/ANB of 10 August 2018 authorising the controlled release of genetically modified sterile male mosquitoes) and institutional ethical permission from the Institutional Ethics Committee for Research in Health Sciences: CEIRES (No. A-003/2019-CEIRES granted on January 9th 2019) and a programme of engagement established community acceptance. Details of the extensive stakeholder and communication processes and activities that were conducted in preparation of this release will be published elsewhere. The village of Bana is located in Western Burkina Faso (12°36′00″N, 3°28′59″W), 23 km west of the city of Bobo-Dioulasso.Bana has two main inhabited agglomerations of similar size: Bana Centre (administrative area) and Bana Market (economic area), separated by a 1.5 km unpopulated land band, crossed by a small semi-permanent river and a forest (see Supplementary Fig. 5). In its entirety, the village comprises about 130 compounds for about 759 inhabitants (local census, IRSS 2014). This region is characterised by two seasons: a wet season from June to September and a dry season from November to April. The mean annual rainfall in the village is about 800 mm and the mean temperature is about 27 °C (22–32 °C)52.Study designThe study design followed the format of an MRR experiment with an intensive period of recaptures followed by several months of monitoring to confirm the disappearance of the transgene. Both the period (July) and design (MRR-like experiment) were informed by previous baseline entomological studies and MRR experiments conducted in the same village41,52. Given the low population size expected in July and to avoid over-sampling, a lower recapture effort (reduction of daily swarm sampling number) was implemented than in previous MRR studies performed in the same area.41 The month of July corresponds to the start of the rainy season, when regular rains and cooler weather promote mosquito survival, and the target population of A. coluzzii is at a much lower level than later in the rainy season41,52. In July, plant coverage is still sparse and males tend to seek refuge inside houses and can be captured in good numbers via indoor sampling52.GM sterile strain maintenanceThe mosquito strain used in the experiment was the genetically modified mosquito Anopheles coluzzii sterile male strain referred to as Ac(DSM)2 (for Anopheles coluzzii Dominant Sterile male strain 2). This strain is the product of local introgression (series of backcrosses) of the original Ag(DSM)2 (dominant sterile male on Anopheles gambiae G3 mosquitoes strain 2) with a local A. coluzzii wild-type (WT) colony (female DSM-carrier crossed with male WT)34. The importation of Ag(DSM)2 in Burkina Faso, introgression with local wild type background and maintenance were conducted under regulatory authorisation from the National Biosafety Agency (N°000002/MRSI/SG/ANB of October 21th 2016). The wild-type A. coluzzii strain used for introgression and maintenance of Ac(DSM)2 was colonised in July 2014 from gravid female adults collected in village 7 of the Kou valley (VK7) in western Burkina Faso. Both colonies were maintained in a dedicated ACL2 (Arthropod Containment Level 2) insectary located within the IRSS main campus at Bobo-Dioulasso, Burkina Faso.For general stock-keeping purposes, Ac(DSM)2 was reared in a dedicated and highly secured climate-controlled room at a temperature fixed at 27.4 °C (±0.2, 95% Confidence intervals) and a relative humidity of 76.3% (±3.2, 95% CIs). Rearing rooms have natural light via windows and were supplemented with an artificial lighting regime of LD 12/12 h photoperiod, including dusk (1 h) and dawn (1 h). Larvae were reared in plastic trays (20 × 30 cm) with 1 l of deionized water and fed with an optimised larvae diet regime53. When mosquito larvae reached their level 3 instar (L3) larvae stage they were sorted manually between transgenic and non-transgenic mosquito larvae using a fluorescent stereomicroscope (Olympus SZX7, 2-8 Honduras street, London, United Kingdom) and put in separated trays to continue their development till pupation. At the pupal stage a second round of sorting occurred to separate male and female (sexing) from both strains. The sexing was done under a basic stereomicroscope (Olympus SZX7 basic, 2–8 Honduras street, London, United Kingdom) using a thin soft brush. Pupae from each strain and sex were placed in small plastic cups inside separate fresh adult cages to emerge. Adults were kept in 30 × 30 × 30 cm insect cages (produced locally) and continuously supplied with 10% (w/v) glucose solution (made with deionized water). Each generation, adult female transgenic mosquitoes were mated with male mosquitoes from the wild-type colony and blood-fed with fresh rabbit’s blood, using a membrane feeder (Hemotek® feeder, Hemotek Ltd, Blackburn United Kingdom). Gravid females were allowed to oviposit in plastic Petri dishes containing a wet sponge covered with filter paper. Eggs were collected and hatched in plastic trays. First instar larvae (L1) were then redistributed into several trays to keep similar larvae abundance (about 250 L1 larvae per tray).In accordance with Mendelian inheritance, stock-maintenance crosses between Ac(DSM)2 females and wild type colony males are expected to generate ~50% hemizygous transgenic male and female progeny referred to as Ac(DSM)2 and 50% non-transgenic sibling with a wild-type phenotype referred to as WT-Ac(DSM)2. That the actual phenotypic proportions matched the expected ratio was checked at each generation a part of standard procedures of colony maintenance.Production, sexing, marking and transport of release mosquitoesReleased males were derived from the 41st backcross generation from strain importation. Assuming Mendelian inheritance, the proportion of residual non-local genetic background after so many generations would be negligible (= 0.541).In rearing the release mosquito cohort, some changes were made in the stock-keeping procedure to maximise the fitness of male mosquitoes to be released. These changes aimed to minimise male mosquito handling during the entire process (rearing, sorting, marking and transport). Crucially, no transgenic versus non-transgenic sorting was done at larval stage resulting in a mix of transgenic and non-transgenic sibling males in the release generation. Additionally, to minimise the number of transgenic female mosquitoes released during the study, male versus female sexing was done at both pupae (initial) and adult (complementary) stages leading to a very high sorting accuracy (over 99.5%). Pupae sexing followed the procedure described for stock maintenance. Next, adult sexing focused on removing the few females resulting from errors in pupal sexing. It consisted of removing those rare females from male mosquito cages through inspection by eye of cages and in using a heat source to attract females. Once spotted, these were removed from male release cages using a mouth aspirator.After pupae sexing, male pupae were placed in 25 × 25 × 25 cm emergence cages (made locally and specially designed to fit dimensions of the secured coolboxes used for secure transportation) at a density of ~1400 pupae per cage. Following adult emergence, and over the following days, the cages were inspected by eye daily to check for and remove any females that had not been detected during the pupae sexing process. This procedure led to a total of 15,384 male mosquitoes aged 3–7 days have emerged in 15 cages and ready for marking and release purposes. Screening of ~50 males randomly picked from each emergence cage was conducted in the ACL2 insectary and revealed a slight bias in favour of WT-Ac(DSM)2 sibling males while Ac(DSM)2 male represented 43.3% (39.7–46.9, 95% CIs) of all emerged males. Based on this genotypic ratio, it was estimated that the male release cohort was equivalent to about 6659 transgenic male mosquitoes Ac(DSM)2 and 8725 non-transgenic sibling mosquitoes called WT-Ac(DSM)2 sibling. All males were kept untouched and in the same cages throughout the whole process until being released.The marking process was performed inside the ACL2 insectary facility, and was carried out the day before field release to allow enough time for mosquito recovery, rest and feeding. The environmental conditions were similar to those used during mosquito production. The mosquitoes were marked directly in their cages by using a cloud dye dusting technique. Aside from being fast, this highly efficient marking procedure (100% of mosquitoes successfully marked) was developed to allow the dust-marking of males in their original emergence cages, thereby avoiding male handling and damage during the marking process. This marking technique consisted of injecting pressurised red fluorescent colour powder (Bioquip® Gladwick Rancho Dominguez, CA 90220, USA; Ref: 1162R) into the cages by using a 5 ml syringe and needle to create a cloud of powder. The cages were wrapped with aluminium foil on all sides to prevent the dust from escaping through the meshed walls. Forceful injection of small amounts of powder from different sides of the cages through the aluminium cover and side netting created a dense cloud of fluorescent powder inside the cages to mark all the mosquitoes. Following marking, sugar-water was available ad-libitum to all marked mosquitoes until field release.About 2 h before the release time, the marked mosquitoes within the mosquito cages were transferred from the IRSS insectary to the release site in Bana village. Before leaving the IRSS insectary, the mosquito cages were covered by a second layer of mosquito net for security purposes. The cages were then wrapped with damp towels and placed in lockable cool boxes dedicated to their transport into the field. After having been secured, the cool boxes containing marked mosquitoes were transported to the release site. The entire process complied carefully with all regulatory requirements related to the permissions received for maintenance, handling and the release of these genetically modified organisms in Burkina Faso.Release phaseAll marked mosquitoes were released on the same day at around 5 pm (about one hour before swarming) in the centre of Bana village by opening the travel cages and allowing free exodus. Mosquitoes that did not leave were counted and subtracted from the released total (n = 534, 3.5%). Taking into account mortality and based on the ratio of Ac(DSM)2 and their siblings previously established, a total of 14,850 male mosquitoes were effectively released, with estimated numbers of 6428 hemizygous transgenic male A. coluzzii mosquitoes Ac(DSM)2 and 8422 non-transgenic WT-Ac(DSM)2 siblings.Recapture phaseMosquito recapture activities started the same day of release (about 2 h after mosquito release) and took place daily for a period of 20 days after release. Two different recapture methods were used: swarm collections using sweep nets (SWN) and pesticides spray catches (PSC) inside houses.Swarm sampling started on the evening of the release day using a well-established sweep net collection method47,54. Previous surveys in the same village41 had allowed mapping of swarm location or natural markers where swarming repeatedly occurs. To ensure sampling across the whole study area, a stratified randomised sampling procedure was used to select and sample 15 mosquito swarms daily at dusk using the sweep net collection method. The area of Bana village and Bana Marché were divided in six and four zones, respectively. Zone 1 and 2 in Bana Village are areas of high swarm abundance and the design ensured that these were not over-represented in swarm collections. Each evening, the teams of capturers set-out to collect up to five swarms per zones depending on swarm availability (swarms are fewer and smaller in early July than later in the month). All mosquitoes captured in the swarms were transported in their sweep nets to the field laboratory and frozen until the next morning for processing. At this stage, a random sample of 15 swarms each day was picked for dust screening and genetic analyses.Pyrethroid spray catches started the morning following the release and continued for 19 days. A set of 20 compounds were sampled each day. The sampling design followed that established in baseline studies leading to the release and in previous MRR studies41. Ten of the compounds were selected completely randomly and the other ten are a fixed set of compounds distributed regularly across the whole village. For each compound selected, a single room (1 sleeping room) within one of the house of compound was chosen for sampling. Although some compounds were selected more than once during the recapture period days, a different room (from a different house inside the same compound when applicable) was selected and no room was sampled twice during the survey period.Pyrethroid spray catches started the morning following the release and continued for 19 days. A set of 20 compounds were selected randomly each day. For each compound selected, a single room (sleeping room) was chosen for sampling. Although some compounds were selected more than once during the seven days, a different room (from a different house inside the same compound when applicable) was selected and no room was sampled twice during the survey period.Captured mosquitoes were identified morphologically in the field using adult anopheline morphological identification keys developed by Holstein55 and a field stereomicroscope (Perfex Sciences® Zoom Pro, Reference: S0852Z5 Toulouse, France). All An. gambiae s.l. mosquitoes were counted, checked for fluorescent dust marking using a Biofinder portable ultraviolet illuminator (Vansky, Shenzhen, China) and preserved in 80% ethanol. The identification of each marked mosquito was confirmed independently by two well-trained members of the staff before conservation in individual 1.5 ml storage microtubes for further analysis. The non-dusted wild Anopheles mosquitoes were pooled (10 individuals per tube) and stored in similar conditions. The location of each collection was recorded and mapped using a GPS (Garmin GPS) device, series GPSMAP®62.2.3. For all recaptured mosquitoes, we calculated the straight line distance from the release point to the recapture location using a Euclidean dispersal distance56. In the present case, the space was assimilated to a two-dimensional orthogonal axis system where xl and yl represent the coordinates of the release point and xr and yr represent the coordinates of the recapture point56. Calculation of the estimated flight distance of the mosquitoes then used the following formula:$${EFD}=sqrt{{left({x}_{r}-{x}_{l}right)}^{2}+{left({y}_{r}-{y}_{l}right)}^{2}}$$
(1)
Ac(DSM)2 male identificationMolecular analysis of recaptured marked mosquitoes was performed by PCR, to identify the Ac(DSM)2 strain and distinguish them from their non-transgenic WT-Ac(DSM)2 siblings. This PCR analysis consisted of detecting the integration of the eGFP::I-PpoI of the DSM transgene which characterised the transgenic mosquito strain Ac(DSM)2. In addition, a molecular species-diagnostic was performed concomitantly using the PCR technique based on the detection of SINE 200× locus57 and this PCR served as a control for DNA integrity. Each mosquito was split into two parts (abdomen and thorax) using forceps. The abdomen was used for the PCR and processed for DNA extraction using ‘squish’ buffer (PCR reaction buffer). The thorax was stored in 80% ethanol at −20 °C. For each mosquito analysed, the same DNA extract was used for both eGFP::I-PpoI transgene detection (identification of Ac(DSM)2 transgenic mosquito) and SINE 200X locus detection (for specie identification and DNA quality control). The Ac(DSM)2 construct was detected using the primers: pBacR-fwd [ATCGGTCTGTATATCGAGGTTTATT] and pBacR-Rev [CTCTAATATTTTGCCAAATGAAGTGCC] targeting the piggyBacR region required for insertion of the transgene. PCR reactions used the Gotaq® PCR kit (GoTaq® G2 Flexi DNA Polymerase, reference: M829B, Promega Corporation, 2800 Woods Hollow Road·Madison, WI 53711-5399, USA).Monitoring of Ac(DSM)2 non-persistenceMonthly mosquito collections were carried out using PSC and swarm sampling to confirm the disappearance of the Ac(DSM)2 transgene from the release site. Monitoring collections were conducted monthly for seven months. This period of monitoring was justified by the regulatory requirement of describing the Ac(DSM)2 disappearance through failure to detect the Ac(DSM)2 transgene for a minimum period of three consecutive months and with high statistical power. During each month of survey, a randomised selection of 20 houses (one room per house) and 20 swarms was sampled. All collected mosquitoes were identified morphologically using identification keys and a field stereomicroscope. Mosquitoes from A. gambiae complex were counted and preserved in 80% (v/v) ethanol for subsequent molecular identification. Each month, a representative sample of collected mosquitoes (up to 300 when available, from both PSC and swarm sampling) was analysed using the Ac(DSM)2-specific and species-specific PCR diagnostics described above to detect whether any A. gambiae s.l. mosquitoes were carrying the DSM transgene.Bayesian inference of mosquito survival, movement and population sizeWe fitted the recapture data to a diffusion model to further investigate dispersal and survival of the marked Ac(DSM)2 and their sibling males, and also to estimate the number of mosquitoes in the background population. This model assumes that the released mosquitoes tend to move in a random manner, meaning they repeatedly take short randomly directed flights that are independent of one another and of the environment. As described below, however, our estimation procedure does also allow for small additional movements where mosquitoes are attracted into nearby swarms at swarming time (dusk), or nearby houses for resting behaviour.We write the diffusion equation as$${partial }_{t}u=D{partial }_{x}^{2}u,$$
(2)
where (u(x,t)) is the probability density of the location of a single marked mosquito at location (x) and time (t), conditional on the individual being alive, and (D) is the diffusion coefficient. Assuming a point release at time (t=0), the above equation has solution$$uleft(r,tright)=frac{{e}^{-frac{{r}^{2}}{4Dt}}}{4,pi D,t}$$
(3)
where (r) is the distance from the release point. We next assume that the released mosquitoes have a constant survival probability of (s) per day, so that the expected number of extant released mosquitoes on day (d) is (R{s}^{d}) where (R) is the number that were released. The expected number of released mosquitoes in a small area ({dA}) is then given by$$qleft(r,dright)=R{s}^{d}frac{{e}^{-frac{{r}^{2}}{4Dd}}}{4,pi D,d}{dA}.$$
(4)
We take three further steps to convert this equation for (q(r,t)) into a likelihood function for the spatio-temporal distribution of recaptures of either Ac(DSM)2 or their sibling males. First, we pool the recaptures on a given day, and made by a given method (either swarm sampling or PSC), by partitioning the study area into annuli centred on the release location. These annuli are the recapture regions, and the expected number of extant marked mosquitoes in a given annulus is the integral of (q(r,d)) over that annulus. This step, therefore, averages out the expected number of marked mosquitoes from the inner to the outer radius of each annulus, and the annulus widths set the scale at which small movements towards swarms or houses, where mosquitoes may be recaptured, are assumed to occur in addition to random movements that underpin the diffusion model. We set the width of each annulus to 50 m, based on our judgement that this distance balances the capacity to separate recaptures at different distances (this capacity reduces with width), with the confidence that movements towards swarms or houses will largely remain within annuli (this confidence increases with width).Second, we assume the observation probability of mosquitoes in a given sample (representing an annulus, capture method, and day), is the number of unmarked mosquitoes in the sample divided by the (unknown) unmarked population size in that annulus. The unmarked population is assumed to have a uniform density, that we will infer alongside the mobility and survival parameters. Finally, we assume the number of marked mosquitoes in a given sample is Poisson-distributed around the expected number.For the data from each recapture method, we used the likelihood function to sample a posterior distribution for the diffusion coefficients and survival rates of the two types of released male mosquitoes, and the density of the unmarked population. We assumed uniform priors with respect to all five parameters and used a Markov chain Monte Carlo algorithm based on Metropolis-Hastings sampling to sample the posterior distribution directly from the log-likelihood. For each analysis (swarm or PSC), we sampled for 100,000 iterations, of which we discarded the initial 20,000 as a transient and thinned the remainder by 100, giving 800 samples in total.Statistical analysisData were analysed using the software JMP 14 (SAS Institute, Inc.). All data were checked for deviations from normality and heterogeneity, and analyses were conducted using parametric and non-parametric methods as appropriate. General linear modelling with Poisson distribution was used to describe male recaptures as a function of genotype and time. Kruskall-Wallis and Mann-Whitney test was used to describe respectively male participation in swarm and Euclidian dispersal distances. General linear modelling with Poisson distribution was used to describe male recaptures as a function of genotype and time. Estimates of population size, survival, and mobility were calculated using a Bayesian approach as described above.Reporting summaryFurther information on research design is available in the Nature Research Reporting Summary linked to this article. More

Van Marle, M. J. E. et al. Historic global biomass burning emissions for CMIP6 (BB4CMIP) based on merging satellite observations with proxies and fire models (1750–2015). Geosci. Model Dev. 10, 3329–3357 (2017).
Google Scholar 
Erb, K. H. et al. Unexpectedly large impact of forest management and grazing on global vegetation biomass. Nature 553, 73–76 (2018).
Google Scholar 
Cook-Patton, S. C. et al. Mapping carbon accumulation potential from global natural forest regrowth. Nature 585, 545–550 (2020).
Google Scholar 
Bastin, J. F. et al. The global tree restoration potential. Science 365, 76–79 (2019).
Google Scholar 
Bowman, D. M. J. S. et al. Fire in the Earth system. Science 324, 481–485 (2009).
Google Scholar 
Archibald, S. et al. Biological and geophysical feedbacks with fire in the Earth system. Environ. Res. Lett. 13, 033003 (2018).
Google Scholar 
Mills, B. J. W., Belcher, C. M., Lenton, T. M. & Newton, R. J. A modeling case for high atmospheric oxygen concentrations during the Mesozoic and Cenozoic. Geology 44, 1023–1026 (2016).
Google Scholar 
Lenton, T. M. in Fire Phenomena and the Earth System: An Interdisciplinary Guide to Fire Science (ed. Belcher, C. M.) 298–308 (Wiley, 2013).Pechony, O. & Shindell, D. T. Driving forces of global wildfires over the past millennium and the forthcoming century. Proc. Natl Acad. Sci. USA 107, 19167–19170 (2010).
Google Scholar 
Marlon, J. R. et al. Reconstructions of biomass burning from sediment-charcoal records to improve data-model comparisons. Biogeosciences 13, 3225–3244 (2016).
Google Scholar 
Archibald, S., Staver, A. C. & Levin, S. A. Evolution of human-driven fire regimes in Africa.Proc. Natl Acad. Sci. USA 109, 847–852 (2012).
Google Scholar 
Santín, C. et al. Towards a global assessment of pyrogenic carbon from vegetation fires. Global Change Biol. 22, 76–91 (2016).
Google Scholar 
Jones, M. W., Santín, C., van der Werf, G. R. & Doerr, S. H. Global fire emissions buffered by the production of pyrogenic carbon. Nat. Geosci. 12, 742–747 (2019).
Google Scholar 
Bird, M. I., Wynn, J. G., Saiz, G., Wurster, C. M. & McBeath, A. The pyrogenic carbon cycle. Annu. Rev. Earth Planet. Sci. 43, 273–298 (2015).
Google Scholar 
Hammes, K. & Abiven, S. in Fire Phenomena and the Earth System: An Interdisciplinary Guide to Fire Science (ed. Belcher, C. M.) 157–176 (Wiley, 2013).Schmidt, M. W. I. et al. Persistence of soil organic matter as an ecosystem property. Nature 478, 49–56 (2011).
Google Scholar 
Lavallee, J. M. et al. Selective preservation of pyrogenic carbon across soil organic matter fractions and its influence on calculations of carbon mean residence times. Geoderma 354, 113866 (2019).
Google Scholar 
Coppola, A. I. et al. Global-scale evidence for the refractory nature of riverine black carbon. Nat. Geosci. 11, 584–588 (2018).
Google Scholar 
Kuzyakov, Y., Bogomolova, I. & Glaser, B. Biochar stability in soil: decomposition during eight years and transformation as assessed by compound-specific 14C analysis. Soil Biol. Biochem. 70, 229–236 (2014).
Google Scholar 
Singh, B. P., Cowie, A. L. & Smernik, R. J. Biochar carbon stability in a clayey soil as a function of feedstock and pyrolysis temperature. Environ. Sci. Technol. 46, 11770–11778 (2012).
Google Scholar 
Masiello, C. A. & Druffel, E. R. M. Black carbon in deep-sea sediments. Science 280, 1911–1913 (1998).
Google Scholar 
Santos, F., Torn, M. S. & Bird, J. A. Biological degradation of pyrogenic organic matter in temperate forest soils. Soil Biol. Biochem. https://doi.org/10.1016/j.soilbio.2012.04.005 (2012).Zimmermann, M. et al. Rapid degradation of pyrogenic carbon. Glob. Change Biol. 18, 3306–3316 (2012).
Google Scholar 
Jones, M. W. et al. Fires prime terrestrial organic carbon for riverine export to the global oceans. Nat. Commun. 11, 2791 (2020).
Google Scholar 
Qi, Y. et al. Dissolved black carbon is not likely a significant refractory organic carbon pool in rivers and oceans. Nat. Commun. 11, 5051 (2020).
Google Scholar 
Pausas, J. G. & Paula, S. Fuel shapes the fire-climate relationship: evidence from Mediterranean ecosystems. Glob. Ecol. Biogeogr. 21, 1074–1082 (2012).
Google Scholar 
Archibald, S., Lehmann, C. E. R., Gómez-Dans, J. L. & Bradstock, R. A. Defining pyromes and global syndromes of fire regimes. Proc. Natl Acad. Sci. USA 110, 6442–6447 (2013).
Google Scholar 
Abatzoglou, J. T., Williams, A. P., Boschetti, L., Zubkova, M. & Kolden, C. A. Global patterns of interannual climate-fire relationships. Glob. Change Biol. 24, 5164–5175 (2018).
Google Scholar 
Brando, P. M. et al. Prolonged tropical forest degradation due to compounding disturbances: implications for CO2 and H2O fluxes. Glob. Change Biol. 25, 2855–2868 (2019).
Google Scholar 
Silva, C. V. J. et al. Drought-induced Amazonian wildfires instigate a decadal-scale disruption of forest carbon dynamics. Phil. Trans. R. Soc. B 373, 20180043 (2018).
Google Scholar 
Withey, K. et al. Quantifying immediate carbon emissions from El Niño-mediated wildfires in humid tropical forests. Phil. Trans. R. Soc. B 373, 20170312 (2018).
Google Scholar 
Pellegrini, A. F. A. et al. Fire frequency drives decadal changes in soil carbon and nitrogen and ecosystem productivity. Nature 553, 194–198 (2018).
Google Scholar 
Reisser, M., Purves, R. S., Schmidt, M. W. I. & Abiven, S. Pyrogenic carbon in soils: a literature-based inventory and a global estimation of its content in soil organic carbon and stocks.Front. Earth Sci. 4, 80 (2016).
Google Scholar 
Wei, X., Hayes, D. J., Fraver, S. & Chen, G. Global pyrogenic carbon production during recent decades has created the potential for a large, long-term sink of atmospheric CO2. J. Geophys. Res. Biogeosci. 123, 3682–3696 (2018).
Google Scholar 
Guimberteau, M. et al. ORCHIDEE-MICT (v8.4.1), a land surface model for the high latitudes: model description and validation. Geosci. Model Dev. 11, 121–163 (2018).
Google Scholar 
Thonicke, K. et al. The influence of vegetation, fire spread and fire behaviour on biomass burning and trace gas emissions: results from a process-based model. Biogeosciences 7, 1991–2011 (2010).
Google Scholar 
Yue, C. et al. Modelling the role of fires in the terrestrial carbon balance by incorporating SPITFIRE into the global vegetation model ORCHIDEE—Part 1: simulating historical global burned area and fire regimes. Geosci. Model Dev. 7, 2747–2767 (2014).
Google Scholar 
Abiven, S. & Santín, C. Editorial: From fires to oceans: dynamics of fire-derived organic matter in terrestrial and aquatic ecosystems. Front. Earth Sci 7, 31 (2019).
Google Scholar 
Santín, C., Doerr, S. H., Preston, C. M. & González-Rodríguez, G. Pyrogenic organic matter production from wildfires: a missing sink in the global carbon cycle. Glob. Change Biol. 21, 1621–1633 (2015).
Google Scholar 
Santín, C. et al. Carbon sequestration potential and physicochemical properties differ between wildfire charcoals and slow-pyrolysis biochars. Sci. Rep. 7, 11233 (2017).
Google Scholar 
Andela, N. et al. A human-driven decline in global burned area. Science 356, 1356–1362 (2017).
Google Scholar 
Arora, V. K. & Melton, J. R. Reduction in global area burned and wildfire emissions since 1930s enhances carbon uptake by land. Nat. Commun. 9, 1326 (2018).
Google Scholar 
Mouillot, F. & Field, C. B. Fire history and the global carbon budget: a 1° × 1° fire history reconstruction for the 20th century. Global Change Biol. 11, 398–420 (2005).
Google Scholar 
Gibson, D. Grasses and Grassland Ecology. Annals of Botany (Oxford Univ. Press, 2009).Dixon, A. P., Faber-Langendoen, D., Josse, C., Morrison, J. & Loucks, C. J. Distribution mapping of world grassland types. J. Biogeogr. 41, 2003–2019 (2014).
Google Scholar 
Bond, W. J. Ancient grasslands at risk. Science 351, 120–122 (2016).
Google Scholar 
Retallack, G. J. Global cooling by grassland soils of the geological past and near future. Annu. Rev. Earth Planet. Sci. 41, 69–86 (2013).
Google Scholar 
Leys, B. A., Marlon, J. R., Umbanhowar, C. & Vannière, B. Global fire history of grassland biomes. Ecol. Evol. 8, 8831–8852 (2018).
Google Scholar 
Alvarado, S. T., Andela, N., Silva, T. S. F. & Archibald, S. Thresholds of fire response to moisture and fuel load differ between tropical savannas and grasslands across continents. Glob. Ecol. Biogeogr. 29, 331–344 (2020).
Google Scholar 
Buisson, E. et al. Resilience and restoration of tropical and subtropical grasslands, savannas and grassy woodlands. Biol. Rev. 94, 590–609 (2019).
Google Scholar 
Rodionov, A. et al. Black carbon in grassland ecosystems of the world. Glob. Biogeochem. Cycles 24, GB3013 (2010).
Google Scholar 
Haberl, H., Erb, K. H. & Krausmann, F. Human appropriation of net primary production: patterns, trends and planetary boundaries. Annu. Rev. Environ. Resources 39, 363–391 (2014).
Google Scholar 
Medan, D., Torretta, J. P., Hodara, K., de la Fuente, E. B. & Montaldo, N. H. Effects of agriculture expansion and intensification on the vertebrate and invertebrate diversity in the Pampas of Argentina. Biodivers. Conserv. 20, 3077–3100 (2011).
Google Scholar 
González-Roglich, M., Swenson, J. J., Villarreal, D., Jobbágy, E. G. & Jackson, R. B. Woody plant-cover dynamics in Argentine savannas from the 1880s to 2000s: the interplay of encroachment and agriculture conversion at varying scales. Ecosystems 18, 481–492 (2015).
Google Scholar 
Satir, O. & Erdogan, M. A. Monitoring the land use/cover changes and habitat quality using Landsat dataset and landscape metrics under the immigration effect in subalpine eastern Turkey. Environ. Earth Sci. 75, 1118 (2016).
Google Scholar 
Şekercioĝlu, Ç. H. et al. Turkey’s globally important biodiversity in crisis. Biol. Conserv. 144, 2752–2769 (2011).
Google Scholar 
Schierhorn, F. et al. Post-Soviet cropland abandonment and carbon sequestration in European Russia, Ukraine and Belarus. Glob. Biogeochem. Cycles 27, 1175–1185 (2013).
Google Scholar 
Jaglan, M. S. & Qureshi, M. H. Irrigation development and its environmental consequences in arid regions of India. Environ. Manage. 20, 323–336 (1996).
Google Scholar 
Joshi, A. A., Sankaran, M. & Ratnam, J. ‘Foresting’ the grassland: historical management legacies in forest-grassland mosaics in southern India, and lessons for the conservation of tropical grassy biomes. Biol. Conserv. 224, 144–152 (2018).
Google Scholar 
Huang, F., Wang, P. & Zhang, J. Grasslands changes in the Northern Songnen Plain, China during 1954–2000. Environ. Monit. Assess. 184, 2161–2175 (2012).
Google Scholar 
Zhou, Y., Hartemink, A. E., Shi, Z., Liang, Z. & Lu, Y. Land use and climate change effects on soil organic carbon in north and northeast China. Sci. Total Environ. 647, 1230–1238 (2019).
Google Scholar 
Williams, N. S. G. Environmental, landscape and social predictors of native grassland loss in western Victoria, Australia. Biol. Conserv. 137, 308–318 (2007).
Google Scholar 
Dowling, P. M. et al. Effect of continuous and time-control grazing on grassland components in south-eastern Australia. Aust. J. Exp. Agric. 45, 369–382 (2005).
Google Scholar 
DeLuca, T. H. & Zabinski, C. A. Prairie ecosystems and the carbon problem. Front. Ecol. Environ. 9, 407–413 (2011).
Google Scholar 
Ceballos, G. et al. Rapid decline of a grassland system and its ecological and conservation implications. PLoS ONE 5, e8562 (2010).
Google Scholar 
Haugo, R. et al. A new approach to evaluate forest structure restoration needs across Oregon and Washington, USA. For. Ecol. Manage. https://doi.org/10.1016/j.foreco.2014.09.014 (2015).DeLuca, T. H. & Aplet, G. H. Charcoal and carbon storage in forest soils of the Rocky Mountain West. Front. Ecol. Environ. 6, 18–24 (2008).
Google Scholar 
Walker, X. J. et al. Increasing wildfires threaten historic carbon sink of boreal forest soils. Nature 572, 520–523 (2019).
Google Scholar 
Bellè, S. L. et al. Key drivers of pyrogenic carbon redistribution during a simulated rainfall event. Biogeosciences 18, 1105–1126 (2021).
Google Scholar 
Abney, R. B., Jin, L. & Berhe, A. A. Soil properties and combustion temperature: controls on the decomposition rate of pyrogenic organic matter. Catena 182, 104127 (2019).
Google Scholar 
Bradstock, R. A., Hammill, K. A., Collins, L. & Price, O. Effects of weather, fuel and terrain on fire severity in topographically diverse landscapes of south-eastern Australia. Landsc. Ecol. 25, 607–619 (2010).
Google Scholar 
Rogers, B. M., Soja, A. J., Goulden, M. L. & Randerson, J. T. Influence of tree species on continental differences in boreal fires and climate feedbacks. Nat. Geosci. 8, 228–234 (2015).
Google Scholar 
Coppola, A. I. & Druffel, E. R. M. Cycling of black carbon in the ocean. Geophys. Res. Lett. 43, 4477–4482 (2016).
Google Scholar 
Stenzel, J. E. et al. Fixing a snag in carbon emissions estimates from wildfires. Glob. Change Biol. 25, 3985–3994 (2019).
Google Scholar 
Murphy, B. P., Prior, L. D., Cochrane, M. A., Williamson, G. J. & Bowman, D. M. J. S. Biomass consumption by surface fires across Earth’s most fire prone continent. Glob. Change Biol. 25, 254–268 (2019).
Google Scholar 
Brando, P. M. et al. Droughts, wildfires and forest carbon cycling: a pantropical synthesis. Annu. Rev. Earth Planet. Sci. 47, 555–581 (2019).
Google Scholar 
Appezzato-da-Glória, B., Cury, G., Soares, M. K. M., Rocha, R. & Hayashi, A. H. Underground systems of Asteraceae species from the Brazilian Cerrado. J. Torrey Bot. Soc. 135, 103–113 (2008).
Google Scholar 
Belcher, C. M. et al. The rise of angiosperms strengthened fire feedbacks and improved the regulation of atmospheric oxygen. Nat. Commun. 12, 503 (2021).
Google Scholar 
Barbero, R., Abatzoglou, J. T., Larkin, N. K., Kolden, C. A. & Stocks, B. Climate change presents increased potential for very large fires in the contiguous United States. Int. J. Wildl. Fire 24, 892–899 (2015).
Google Scholar 
Stephens, S. L. et al. Managing forests and fire in changing climates. Science 342, 41–42 (2013).
Google Scholar 
Trenberth, K. E. Changes in precipitation with climate change. Clim. Res. 47, 123–138 (2011).
Google Scholar 
Prein, A. F. et al. The future intensification of hourly precipitation extremes. Nat. Clim. Change 7, 48–52 (2017).
Google Scholar 
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).
Google Scholar 
Silveira, F. A. O. et al. Myth-busting tropical grassy biome restoration. Restor. Ecol. 28, 1067–1073 (2020).
Google Scholar 
Strassburg, B. B. N. et al. Global priority areas for ecosystem restoration. Nature 586, 724–729 (2020).
Google Scholar 
Schmidt, H. P. et al. Pyrogenic carbon capture and storage. GCB Bioenergy 11, 573–591 (2019).
Google Scholar 
Fu, Z. et al. Recovery time and state change of terrestrial carbon cycle after disturbance. Environ. Res. Lett. 12, 104004 (2017).
Google Scholar 
Zhu, D. et al. Improving the dynamics of Northern Hemisphere high-latitude vegetation in the ORCHIDEE ecosystem model. Geosci. Model Dev. 8, 2263–2283 (2015).
Google Scholar 
Zhu, D. et al. Simulating soil organic carbon in Yedoma deposits during the Last Glacial Maximum in a land surface model. Geophys. Res. Lett. 43, 5133–5142 (2016).
Google Scholar 
Krinner, G. et al. A dynamic global vegetation model for studies of the coupled atmosphere-biosphere system. Glob. Biogeochem. Cycles 19, GB1015 (2005).
Google Scholar 
Yue, C., Ciais, P., Cadule, P., Thonicke, K. & Van Leeuwen, T. T. Modelling the role of fires in the terrestrial carbon balance by incorporating SPITFIRE into the global vegetation model ORCHIDEE—Part 2: carbon emissions and the role of fires in the global carbon balance. Geosci. Model Dev. 8, 1321–1338 (2015).
Google Scholar 
Hantson, S. et al. The status and challenge of global fire modelling. Biogeosciences 13, 3359–3375 (2016).
Google Scholar 
Hantson, S. et al. Quantitative assessment of fire and vegetation properties in simulations with fire-enabled vegetation models from the Fire Model Intercomparison Project. Geosci. Model Dev. 13, 3299–3318 (2020).
Google Scholar 
Li, F. et al. Historical (1700–2012) global multi-model estimates of the fire emissions from the Fire Modeling Intercomparison Project (FireMIP). Atmos. Chem. Phys. 19, 12545–12567 (2019).
Google Scholar 
Forkel, M. et al. Emergent relationships with respect to burned area in global satellite observations and fire-enabled vegetation models. Biogeosciences 16, 57–76 (2019).
Google Scholar 
Parton, W. J., Stewart, J. W. B. & Cole, C. V. Dynamics of C, N, P and S in grassland soils: a model. Biogeochemistry 5, 109–131 (1988).
Google Scholar 
Singh, N. et al. Transformation and stabilization of pyrogenic organic matter in a temperate forest field experiment. Glob. Change Biol. 20, 1629–1642 (2014).
Google Scholar 
Viovy, N. CRUNCEP Version 7—Atmospheric Forcing Data for the Community Land Model (Research Data Archive at the National Center for Atmospheric Research, Computational and Information Systems Laboratory, 2018); https://doi.org/10.5065/PZ8F-F017Mckee, T. B. T. et al. The relationship of drought frequency and duration to time scales. In Proc. Eighth Conference on Applied Climatology 179–184 (American Meteorological Society, 1993).The NCAR Command Language, Version 6.6.2 (UCAR/NCAR/CISL/TDD, 2019).Freeborn, P. H., Wooster, M. J., Roy, D. P. & Cochrane, M. A. Quantification of MODIS fire radiative power (FRP) measurement uncertainty for use in satellite-based active fire characterization and biomass burning estimation. Geophys. Res. Lett. 41, 1988–1994 (2014).
Google Scholar 
Giglio, L. MODIS Collection 5 Active Fire Product User’s Guide Version 2.5 (Science Systems and Applications, 2013).Huang, N. et al. Spatial and temporal variations in global soil respiration and their relationships with climate and land cover. Sci. Adv. 6, eabb8508 (2020).
Google Scholar 
Warner, D. L., Bond-Lamberty, B., Jian, J., Stell, E. & Vargas, R. Spatial predictions and associated uncertainty of annual soil respiration at the global scale. Glob. Biogeochem. Cycles 33, 1733–1745 (2019).
Google Scholar  More