Ecology

Clare, E., Barber, B., Sweeney, B., Hebert, P. & Fenton, M. Eating local: Influences of habitat on the diet of little brown bats (Myotis lucifugus). Mol. Ecol. 20, 1772–1780 (2011).CAS 
PubMed 
Article 

Google Scholar 
Larter, N. C. & Gates, C. C. Diet and habitat selection of wood bison in relation to seasonal changes in forage quantity and quality. Can. J. Zool. 69, 2677–2685 (1991).Article 

Google Scholar 
Veloso, C. & Bozinovic, F. Dietary and digestive constraints on basal energy metabolism in a small herbivorous rodent. Ecology 74, 2003–2010 (1993).Article 

Google Scholar 
Hawke, T., Bates, H., Hand, S., Archer, M. & Broome, L. Dietary analysis of an uncharacteristic population of the Mountain Pygmy-possum (Burramys parvus) in the Kosciuszko National Park, New South Wales, Australia. PeerJ 7, e6307 (2019).PubMed 
PubMed Central 
Article 

Google Scholar 
Pearce-Higgins, J. W. Using diet to assess the sensitivity of northern and upland birds to climate change. Clim. Res. 45, 119–130 (2010).Article 

Google Scholar 
Eitzinger, B. et al. Assessing changes in arthropod predator-prey interactions through DNA-based gut content analysis—variable environment, stable diet. Mol. Ecol. 28, 266–280 (2019).CAS 
PubMed 
Article 

Google Scholar 
Edgar, G. J. Predator-prey interactions in seagrass beds. II. Distribution and diet of the blue manna crab Portunus pelagicus Linnaeus at Cliff Head, Western Australia. J. Exp. Mar. Biol. Ecol. 139, 23–32 (1990).Article 

Google Scholar 
Beck, J. L., Peek, J. M. & Strand, E. K. Estimates of elk summer range nutritional carrying capacity constrained by probabilities of habitat selection. J. Wildl. Manag. 70, 283–294 (2006).Article 

Google Scholar 
DeYoung, R. W., Hellgren, E. C., Fulbright, T. E., Robbins, W. F. Jr. & Humphreys, I. D. Modeling nutritional carrying capacity for translocated desert bighorn sheep in western Texas. Restor. Ecol. 8, 57–65 (2000).Article 

Google Scholar 
Hua, L. et al. Captive breeding of pangolins: current status, problems and future prospects. Zookeys 507, 99–114 (2015).Article 

Google Scholar 
Nielsen, J. M., Clare, E. L., Hayden, B., Brett, M. T. & Kratina, P. Diet tracing in ecology: Method comparison and selection. Methods Ecol. Evol. 9, 278–291 (2017).Article 

Google Scholar 
Galimberti, A. et al. DNA barcoding as a new tool for food traceability. Food Res. Int. 50, 55–63 (2013).CAS 
Article 

Google Scholar 
Soininen, E. M. et al. Shedding new light on the diet of Norwegian lemmings: DNA metabarcoding of stomach content. Polar Biol. 36, 1069–1076 (2013).Article 

Google Scholar 
Rees, G. N., Shackleton, M. E., Watson, G. O., Dwyer, G. K. & Stoffels, R. J. Metabarcoding demonstrates dietary niche partitioning in two coexisting blackfish species. Mar. Freshw. Res. 71(4), 512–517 (2019).Article 

Google Scholar 
Taberlet, P., Coissac, E., Pompanon, F., Brochmann, C. & Willerslev, E. Towards next-generation biodiversity assessment using DNA metabarcoding. Mol. Ecol. 21, 2045–2050 (2012).CAS 
PubMed 
Article 

Google Scholar 
Aylagas, E., Borja, Á., Irigoien, X. & Rodriguez-Ezpeleta, N. Benchmarking DNA metabarcoding for biodiversity-based monitoring and assessment. Front. Mar. Sci. 3, 96 (2016).
Google Scholar 
De Barba, M. et al. DNA metabarcoding multiplexing and validation of data accuracy for diet assessment: Application to omnivorous diet. Mol. Ecol. Resour. 14, 306–323 (2014).PubMed 
Article 

Google Scholar 
Kartzinel, T. R. et al. DNA metabarcoding illuminates dietary niche partitioning by African large herbivores. Proc. Natl. Acad. Sci. 112, 8019–8024 (2015).CAS 
PubMed 
PubMed Central 
ADS 
Article 

Google Scholar 
Lopes, C. et al. DNA metabarcoding diet analysis for species with parapatric vs sympatric distribution: A case study on subterranean rodents. Heredity 114, 525–536 (2015).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
Guillerault, N., Bouletreau, S., Iribar, A., Valentini, A. & Santoul, F. Application of DNA metabarcoding on faeces to identify European catfish Silurus glanis diet. J. Fish Biol. 90, 2214–2219 (2017).CAS 
PubMed 
Article 

Google Scholar 
Jakubavivciute, E., Bergström, U., Eklöf, J. S., Haenel, Q. & Bourlat, S. J. DNA metabarcoding reveals diverse diet of the three-spined stickleback in a coastal ecosystem. PLoS ONE 12, e0186929 (2017).Article 

Google Scholar 
Grant, T. & Fanning, D. Platypus 4th edn. (CSIRO Publishing, 2007).Book 

Google Scholar 
Hawke, T. et al. Long-term movements and activity patterns of platypus on regulated rivers. Sci. Rep. 11, 1–11 (2021).MathSciNet 
Article 

Google Scholar 
Gregory, J., Iggo, A., McIntyre, A. & Proske, U. Receptors in the bill of the platypus. J. Physiol. 400, 349 (1988).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
McLachlan-Troup, T., Dickman, C. & Grant, T. Diet and dietary selectivity of the platypus in relation to season, sex and macroinvertebrate assemblages. J. Zool. 280, 237–246 (2010).Article 

Google Scholar 
Harrop, C. & Hume, I. Digestive tract and digestive function in monotremes and nonmacropod marsupials. Compar. Physiol. Primitive Mamm. 4, 63–77 (1980).
Google Scholar 
Klamt, M., Davis, J. A., Thompson, R. M., Marchant, R. & Grant, T. R. Trophic relationships of the platypus: Insights from stable isotope and cheek pouch dietary analyses. Mar. Freshw. Res. 67, 1196–1204 (2016).CAS 
Article 

Google Scholar 
Faragher, R., Grant, T. & Carrick, F. Food of the platypus (Ornithorhynchus anatinus) with notes on the food of brown trout (Salmo trutta) in the Shoalhaven River, NSW. Austral. J. Ecol. 4, 171–179 (1979).Article 

Google Scholar 
Grant, T. R. Food of the platypus, Ornithorhynchus anatinus (Ornithorhynchidae: Monotremata) from various water bodies in New South Wales. Aust. Mammal. 5, 135–136 (1982).
Google Scholar 
Marchant, R. & Grant, T. The productivity of the macroinvertebrate prey of the platypus in the upper Shoalhaven River, New South Wales. Mar. Freshw. Res. 66, 1128–1137 (2015).Article 

Google Scholar 
Krueger, B., Hunter, S. & Serena, M. Husbandry, diet and behaviour of platypus Ornithorhynchus anatinus at Healesville Sanctuary. Int. Zoo Yearbook 31, 64–71 (1992).Article 

Google Scholar 
Thomas, J. L., Handasyde, K. A., Temple-Smith, P. & Parrott, M. L. Seasonal changes in food selection and nutrition of captive platypuses (Ornithorhynchus anatinus). Aust. J. Zool. 65, 319–327 (2018).Article 

Google Scholar 
Hawke, T., Bino, G. & Kingsford, R. T. Damming insights: Variable impacts and implications of river regulation on platypus populations. Aquat. Conserv. Mar. Freshwat. Ecosyst. 31, 504–519 (2021).Article 

Google Scholar 
Bino, G., Kingsford, R. T., Grant, T., Taylor, M. D. & Vogelnest, L. Use of implanted acoustic tags to assess platypus movement behaviour across spatial and temporal scales. Sci. Rep. 8, 5117 (2018).PubMed 
PubMed Central 
ADS 
Article 

Google Scholar 
Chinnadurai, S. K., Strahl-Heldreth, D., Fiorello, C. V. & Harms, C. A. Best-Practice guidelines for field-based surgery and anesthesia of free-ranging wildlife. I. Anesthesia and Analgesia. J. Wildl. Dis. 52(2 Suppl), S14–27. https://doi.org/10.7589/52.2S.S14 (2016).PubMed 
Article 

Google Scholar 
Fiorello, C. V., Harms, C. A., Chinnadurai, S. K. & Strahl-Heldreth, D. Best-Practice guidelines for field-based surgery and anesthesia on free-ranging wildlife. Ii. Surgery. J. Wildl. Dis. 52(2 Suppl), S28–39. https://doi.org/10.7589/52.2S.S28 (2016).PubMed 
Article 

Google Scholar 
Vogelnest, L. & Woods, R. Medicine of Australian mammals: CSIRO Publishing (2008).Geller, J., Meyer, C., Parker, M. & Hawk, H. Redesign of PCR primers for mitochondrial cytochrome c oxidase subunit I for marine invertebrates and application in all-taxa biotic surveys. Mol. Ecol. Resour. 13, 851–861 (2013).CAS 
PubMed 
Article 

Google Scholar 
Leray, M. et al. A new versatile primer set targeting a short fragment of the mitochondrial COI region for metabarcoding metazoan diversity: Application for characterizing coral reef fish gut contents. Front. Zool. 10, 1–14 (2013).Article 

Google Scholar 
Greenfield, P. Greenfield hybrid analysis pipeline (GHAP). v1 (CSIRO, 2017).
Google Scholar 
Shackleton, M. et al. How does molecular taxonomy for deriving river health indices correlate with traditional morphological taxonomy?. Ecol. Indic. 125, 107537 (2021).Article 

Google Scholar 
Hebert, P. D., Cywinska, A., Ball, S. L. & Dewaard, J. R. Biological identifications through DNA barcodes. Proc. R. Soc. Lond. Ser. B Biol. Sci. 270, 313–321 (2003).CAS 
Article 

Google Scholar 
Ostell, J. & Sayers, E. W. Dennis A. Benson, Mark Cavanaugh, Karen Clark, Ilene Karsch-Mizrachi, David J. Lipman.Wilcoxon, F. Individual comparisons by ranking methods. Biometrics Bulletin, 1. In Breakthroughs in Statistics, 196–202 (Springer, 1992).Oksanen, J. et al. Package “vegan”. Community ecology package, version. Vol 2, No. 9, 1–295. (2013).R Development Core Team. R: A Language and Environment for Statistical Computing (R Foundation for Statistical Computing, 2021).
Google Scholar  More

Epithemia isolation and cultureThe Epithemia cells were isolated from 0.5 L of seawater collected from depths of 25, 75, and 100 m in the North Pacific Subtropical Gyre (22°45′ N, 158°00′ W). Seawater was collected during the near-monthly Hawaii Ocean Time-series (HOT) expeditions to the long-term monitoring site Station ALOHA (water depth ca. 4800 m) in October 2014 (HOT cruise #266) and February–July 2019 (HOT cruises #310–313). Serial dilution (unialgal strains UHM3202, UHM3203, UHM3204) or micropipette isolation of single cells (clonal strains UHM3200, UHM3201, UHM3210, UHM3211) were used to establish the Epithemia cultures, which were grown in a seawater-based, low-nitrogen medium. Filtered (0.2 µm) and autoclaved, undiluted Station ALOHA seawater was amended with 2 μM EDTA, 50 nM ferric ammonium citrate, 7.5 μM phosphoric acid, trace metals (100 nM MnSO4, 10 nM ZnCl2, 10 nM Na2MoO4, 1 nM CoCl2, 1 nM NiCl2, 1 nM Na2SeO3), vitamins (50 μg/L inositol, 10 μg/L calcium pantothenate, 10 μg/L thiamin, 5 μg/L pyridoxine HCl, 5 μg/L nicotinic acid, 0.5 μg/L para-aminobenzoic acid, 0.1 μg/L folic acid, 0.05 μg/L biotin, 0.05 μg/L vitamin B12), and 106 μM Na2SiO3. Although not tested here, simpler formulations of diazotroph media such as PMP40 or RMP41 may also be suitable for growing Epithemia, when made with 100% seawater and adding Na2SiO3. The cultures were subsequently incubated at 24 °C on a 12:12 h light:dark cycle with 50–100 μmol quanta m−2 s−1 using cool white fluorescent bulbs. All E. pelagica and E. catenata symbioses were stable under these medium and incubation conditions. E. pelagica was successfully isolated from at least one of the three depths that were targeted during each sampling occasion.Morphological observationsEpithemia living and fixed cells were imaged by light and epifluorescence microscopy using a Nikon Eclipse 90i microscope at 40×–60× magnification. Diatom cell sizes were determined using >60 live, exponentially growing cells, imaged in either valve view (E. pelagica) or girdle view (E. catenata). Endosymbiont (spheroid body) cell sizes were averaged from DNA-stained cells for E. pelagica UHM3200 (n = 78) and E. catenata UHM3210 (n = 91), imaged by epifluorescence microscopy after preparing samples as follows: Epithemia cells were fixed in 4% glutaraldehyde for 30 min, pelleted at 1000 × g for 1 min, the supernatant was exchanged with 0.5% Triton X-100 (in autoclaved filtered seawater), samples were incubated for 10 min with gentle agitation, cells were then pelleted at 4000 × g for 1 min, supernatant was exchanged with autoclaved filtered seawater and fixed in 4% glutaraldehyde, and samples were stained with 1× final concentration of SYBR Gold nucleic acid stain (Invitrogen, cat. # S11494) for 2 h. For routine observations of endosymbionts (e.g., determining presence/absence and number per host cell), osmotic shock was used to disrupt the cell contents of diatom host cells and improve visualization of the endosymbionts. This was achieved by gently pelleting cells and exchanging the medium with either ultrapure water or 2–3 M NaCl solution, followed by immediate observation. While this is a simple technique for detecting and visualizing endosymbionts (Fig. 1c, f), it does not accurately represent the natural location of endosymbionts within the host cells, as seen when compared to fixed cell preparations for epifluorescence microscopy (Fig. 1n, o). To assess the presence of fluorescent photopigments in endosymbiont cells, live host cells were pelleted at 4000 × g for 5 min and crushed using a microcentrifuge tube pestle (SP Bel-Art, cat. # F19923-0000) to release the endosymbionts. The crushed pellet was resuspended in 75% glycerol containing live Synechococcus WH7803 cells (positive control for fluorescence), and samples were observed by epifluorescence microscopy using filter cubes appropriate for observing phycoerythrin (EX: 551/10, BS: 560, EM: 595/30) and chlorophyll (EX: 480/30, BS: 505, EM: 600LP).The loss of endosymbionts from Epithemia cultures (UHM3200 and UHM3210) was observed after propagating cells for four months in nitrogen-replete medium (K)18, where approximately 5–10% of the culture was transferred to fresh medium about every two weeks. Observations were only made at the end of the four-month period. Endosymbionts were not observed growing freely in these cultures, and the absence of endosymbionts within host cells was confirmed by the failure to observe spheroid bodies by light microscopy after osmotic shock of the diatoms, as well as a failure to amplify the endosymbiont SSU (16S rRNA) and nifH genes from cellular DNA extracts. PCR reactions were performed in parallel with DNA extracts from control cultures (grown in low-nitrogen medium), using the same template DNA amount (10 ng) and PCR conditions (see methods for Marker gene sequencing and phylogenetics).Ultrastructural observations by electron microscopy (EM) were conducted for E. pelagica UHM3200 and E. catenata UHM3210. EM preparations of diatoms typically involve the oxidative removal of organic matter to uncover the fine details of frustule ultrastructure. However, in the case of E. catenata, oxidatively cleaned cells lacked structural integrity, leading to collapsed frustules when dried and viewed by scanning EM (SEM). For this reason, both species were prepared for SEM with and without (Fig. 1a, d) the oxidative removal of organic matter, and cleaned E. catenata frustules were further analyzed by transmission EM (TEM). To remove organic matter, 100 mL of exponentially growing culture was pelleted by centrifugation at 1000 × g for 10 min and resuspended in 30% H2O2. Cells were boiled in H2O2 for 1–2 h, followed by rinsing cells six times in ultrapure water by sequential centrifugation at 1000 × g for 10 min and resuspension of cell pellets. Suspensions of the cleaned cells were dried on aluminum foil and mounted on aluminum stubs with double-sided copper tape. For some E. catenata SEM preparations, the cleaned frustules were dehydrated in an ethanol dilution series and exchanged into hexamethyldisilazane (HMDS) prior to drying on aluminum foil; this was to minimize the collapse of frustules resulting from drying. To prepare cells with organic matter intact, 25 mL of exponentially growing culture was mixed with an equal volume of fixative solution (5% glutaraldehyde, 0.2 M sodium cacodylate pH 7.2, 0.35 M sucrose, 10 mM CaCl2) and incubated overnight at 4 °C. Cells were gently filtered onto a 13 mm diameter 1.2 μm pore size polycarbonate membrane filter (Isopore, Millipore Sigma), washed with 0.1 M sodium cacodylate buffer (pH 7.4, 0.35 M sucrose), fixed with 1% osmium tetroxide in 0.1 M sodium cacodylate (pH 7.4), dehydrated in a graded ethanol series, and critical point dried. Filters were mounted on aluminum stubs with double-sided conductive carbon tape. All SEM stubs were sputter coated with Au/Pd, prior to observing on a Hitachi S-4800 field emission scanning electron microscope at the University of Hawai’i at Mānoa (UHM) Biological Electron Microscope Facility (BEMF). Cleaned E. catenata cells were prepared for TEM by drying a drop of sample on a formvar/carbon-coated grid and observing on a Hitachi HT7700 transmission electron microscope at UHM BEMF.Additional light microscopy of hydrogen-peroxide cleaned frustules was conducted for E. pelagica UHM3201 and E. catenata UHM3210. Samples were mounted in Naphrax (PhycoTech, Inc., cat. # P-Naphrax200) and observed at 100× using an Olympus BX41 Photomicroscope (Olympus America Inc., Center Valley, Pennsylvania) with differential interference contrast optics and an Olympus SC30 Digital Camera at California State University San Marcos.A key to the strains used in each micrograph is provided in Supplementary Table 2.Marker gene sequencing and phylogeneticsFor each Epithemia strain, 25–50 mL of culture was pelleted at 4000 × g for 10 min, and DNA was extracted from the pellet using the ZymoBIOMICS DNA Miniprep Kit (Zymo Research, cat. # D4300). Marker genes were amplified with the Expand High Fidelity PCR System (Roche, cat. # 4743733001), using conditions previously described for genes SSU encoding 18S rRNA (Euk328f/Euk329r)42, LSU encoding 28S rRNA (D1R/D2C)43, rbcL (rbcL66+/dp7−)44,45, psbC (psbC+/psbC−)44, and cob (Cob1f/Cob2r)21. For the endosymbionts, a partial sequence for the SSU (16S rRNA) gene was amplified using a primer set targeting unicellular cyanobacterial diazotrophs, CYA359F/Nitro821R46,47, and the nifH gene was amplified using new primers specific to the nifH of Cyanothece-like organisms, ESB-nifH-F (5′-TACGGAAAAGGCGGTATCGG-3′) and ESB-nifH-R (5′-CACCACCAAGRATACCGAAGTC-3′), with a 55 °C annealing temperature and 75 s extension time. All primers were synthesized by Integrated DNA Technologies (IDT). Amplified products were cloned and transformed into E. coli using the TOPO TA Cloning Kit for Sequencing (Invitrogen, cat. # K457501), and plated colonies were picked and grown in Circlegrow medium (MP Biomedicals, cat. # 113000132). Plasmids were extracted with the Zyppy Plasmid Miniprep kit (Zymo Research, cat. # D4019) and sequenced from the M13 vector primers using Sanger technology at GENEWIZ (South Plainfield, NJ). For the diatom SSU (18S rRNA) gene, sequencing reactions were also performed using the 502f and 1174r primers48.Phylogenetic trees (Fig. 2) were inferred using concatenated alignments for both diatom host genes (SSU encoding 18S rRNA, psbC, rbcL) and endosymbiont genes (SSU encoding 16S rRNA, nifH). For each gene, nucleotide sequences were aligned using MAFFT v7.45349 (L-INS-i method), and sites with gaps or missing data were removed. An appropriate nucleotide substitution model was selected for each gene alignment using jModelTest v2.1.1050. Bayesian majority consensus trees were inferred from the concatenated alignments using MrBayes v3.2.751 with two runs of 4–8 chains, until the average standard deviation of split frequencies dropped below 0.01. Maximum likelihood bootstrap values were generated for the Bayesian tree using RAxML v8.2.1252, implemented with 1000 iterations of rapid bootstrapping. To further analyze the phylogenetic position of the new Epithemia species in the broader context of Surirellales and Rhopalodiales diatoms, individual gene trees (SSU encoding 18S rRNA, LSU, rbcL, psbC, and cob; Supplementary Figs. 13–19) were constructed from sequences aligned using MAFFT (automatic detection method) and trimmed using trimAl v1.253 (gappyout method). rRNA gene phylogenies were also inferred using sequences aligned according to the global SILVA alignment for SSU and LSU genes using SINA54, which were either left untrimmed in the case of the LSU gene or trimmed to remove highly variable positions (SINA’s “012345” positional variability filter) and gappy positions (trimAL v1.2, gappyout method) in the case of the SSU gene. These trimming strategies were selected based on their ability to maximize the monophyly of the previously described Rhopalodiales clade and minimize the separation of known conspecific strains, such as the strains of E. pelagica described here. All gene phylogenies were inferred using the Bayesian methods described above. To investigate the level of support for constrained tree topologies placing E. catenata within or outside of the genus Epithemia and family Rhopalodiaceae, SH55 and AU56 statistical tests were performed in IQ-TREE 257 (implementing ModelFinder58) using all alignments from the individual gene trees (Supplementary Table 3).Given E. catenata’s unusual morphology, test trees were inferred with the inclusion of diatom sequences from orders Bacillariales (Nitzschia, Pseudo-nitzschia), Cymbellales (Didymosphenia), Naviculales (Amphiprora, Navicula, Pinnularia), and Thalassiophysales (Amphora, Halamphora, Thalassiophysa); however, E. catenata was consistently placed within Rhopalodiales, and these trees were not pursued further.An additional nifH phylogeny was constructed using all environmental sequences from NCBI’s non-redundant nucleotide (nt) database >300 bp and sharing >95% nucleotide sequence identity with EpSB and EcSB nifH sequences (Supplementary Fig. 23), including 51 environmental sequences from prior studies investigating marine diazotrophs34,59,60,61,62,63,64,65,66. Environmental nifH sequences were aligned to the previously generated nifH sequence alignment using MAFFT (automatic method detection and addfragments options), and the best-scoring maximum likelihood phylogeny was inferred using RAxML with 1000 iterations of rapid bootstrapping. NCBI accession numbers for all tree sequences are in the Source Data file.Analysis of Epithemia endosymbiont nifH sequences in environmental datasetsNucleotide sequences for EpSB and EcSB nifH were queried against NCBI’s non-redundant nucleotide (nt) database using webBLAST67 (megablast; https://blast.ncbi.nlm.nih.gov/) and SRA databases for nifH amplicon sequencing projects from the marine environment using the SRA Toolkit68 (dc-megablast, with database validation using vdb-validate; https://github.com/ncbi/sra-tools). Database hits with 98–100% nucleotide identity over an alignment of the entire subject sequence (BLAST alignment length = subject sequence length) were identified, and the associated sample’s latitude and longitude coordinates (where available) were mapped. Coordinates were also mapped for metagenome and metatranscriptome samples containing matches to unigene MATOU-v1_93255274 from the Marine Atlas of Tara Oceans Unigenes69, a unigene that shares 100% identity over the entire length of the EpSB UHM3202 nifH sequence and >99.4% identity with all other EpSB nifH sequences.The presence of EpSB and EcSB nifH sequences was examined in metagenomes prepared from sinking particles collected at 4000 m depth at Station ALOHA27,28. The sinking particles were collected during intervals of 12, 10, and 8 days during 2014, 2015, and 2016, respectively, using a McLane sediment trap equipped with a 21-sample bottle carousel. The presence of EpSB and EcSB nifH sequences in the metagenomes was assessed by blastn70, after first removing low quality bases from metagenomic reads using Trimmomatic v0.3971 (parameters: LEADING:20 TRAILING:20 MINLEN:100). For each sediment trap metagenome, the total number of reads matching EpSB or EcSB nifH nucleotide sequences with 100% identity were tallied and normalized to the total number of reads in the database. Only EpSB-matching reads were detected in this analysis.Quantitative PCRSpecific PCR primers were designed targeting a 102 bp region of E. pelagica’s LSU gene (Epel-LSU-F, 5′-GAAACCAGTGCAAGCCAAC-3′; Epel-LSU-R, 5′-AGGCCATTATCATCCCTTGTC-3′) and an 85 bp region EpSB’s nifH gene (EpSB-nifH-F, 5′-CACACTAAAGCACAAACTACC-3′; EpSB-nifH-R, 5′-CAAGTAGTACTTCGTCTAGCTC-3′) and were synthesized by IDT. Gene copy concentrations were quantified for Station ALOHA water samples (~2 L) collected by Niskin bottles at 5, 25, 45, 75, 100, 125, 150, and 175 m on January 16 and July 1 (except 5 m), 2014, during HOT cruises #259 and #264. Samples were filtered onto 25 mm diameter, 0.02 μm pore size aluminum oxide filters (Anotop; Whatman, cat. # WHA68092102) and stored at −80 °C until extracting DNA using the MasterPure Complete DNA and RNA Purification Kit (Epicentre, cat. # MC85200) according to Mueller et al.72. Briefly, a 3-mL syringe filled with 1 mL of tissue and cell lysis solution (MasterPure) containing 100 μg mL−1 proteinase K was attached to the outlet of the filter, and the filter inlet was sealed with a second 3-mL syringe. The lysis solution was pulled halfway through to saturate the filter membrane, and the entire assembly was incubated at 65 °C for 15 min while attached to a rotisserie in a hybridization oven rotating at ca. 16 rpm. The lysis buffer was then drawn fully into the inlet syringe, transferred to a microcentrifuge tube, and placed on ice. The remaining steps for protein precipitation and removal and nucleic acid precipitation were carried out following the manufacturer’s instructions. For each sample, DNA was resuspended in a final volume of 100 μL. Quantitative PCR (qPCR) was performed using the PowerTrack SYBR Green Master Mix system (Applied Biosystems, cat. # A46109) and run on an Eppendorf Mastercycler epgradient S realplex2 real-time PCR machine. Reactions (20 µL total volume) were prepared according to the manufacturer’s protocol, containing 500 nM of each primer. Sample reactions (four replicates) contained 2 μL of environmental DNA extract (24–76 ng DNA), while standards (three replicates) contained 2 μL of gBlocks Gene Fragments (IDT) that were prepared at 1, 2, 3, 4, 5, and 6 log gene copies/μL. The gBlocks Gene Fragments were 500 bp in length and encompassed the entire E. pelagica UHM3201 LSU sequence and positions 1–500 of the EpSB UHM3201 nifH sequence, respectively. The main cycling conditions consisted of an initial denaturation and enzyme activation step of 95 °C for 2 min, followed by 40 cycles of 95 °C for 5 s and 57 °C or 55 °C for 30 s for the LSU and nifH genes, respectively. Melting curves were analyzed to verify the specificity of the amplifications, and reactions containing Epithemia catenata DNA extract were included as negative controls. Reaction efficiencies were 104.23% and 95.15% for the LSU and nifH genes, respectively. The limit of detection for these assays was not empirically determined. gBlocks sequences, qPCR threshold cycle values, and conversion equations are provided in the Source Data file.Physiology experimentsThe daily patterns of N2 fixation were quantified for E. pelagica UHM3200 and E. catenata UHM3210 using two techniques: acetylene (C2H2) reduction to ethylene (C2H4) and argon induced dihydrogen (H2) production (AIHP). Both analyses were conducted using a gaseous flow-through system that quantified the relevant trace gas on the sample outlet line with a temporal resolution of 10 min73. To conduct the measurements, a 10-mL subsample of each Epithemia culture was placed in a 20-mL borosilicate vial and closed using gas-tight rubber stoppers and crimp seals. Separate bottles were used for H2 production and C2H2 reduction. During the experimental period, the temperature was maintained at 25 ± 0.2 °C using a benchtop incubator (Incu-Shaker; Benchmark Scientific) and light exposure was 200 μmol photons m−2 s−1 at wavelengths of 380–780 nm with a 12:12 h square light:dark cycle (Prime HD+; Aqua Illumination). To conduct the AIHP method, the sample vial containing the culture was flushed with a high purity gas mixture consisting of argon (makeup gas; 80%), oxygen (20%), and carbon dioxide (0.04%). In the absence of N2, all of the electrons that would have been used to reduce N2 to NH3 are diverted to H2 production, thereby providing a measure of Total Nitrogenase Activity (TNA). The C2H2 reduction assay also represents a measure of TNA. Our analytical set-up introduced C2H2 at a 1% addition (vol/vol) to the high purity air with a total flow rate (13 mL min−1) identical to the AIHP method. The gas emissions were analyzed using separate reductive trace gas analyzers that were optimized for the quantification of H2 and C2H4. To verify the observed daily patterns in N2 fixation, 15N2 assimilation measurements were conducted on triplicate samples of Epithemia cultures at targeted time points. Five milliliters of 15N-enriched seawater was added to the subsamples, which were subsequently crimp sealed and incubated for a 2 h period with the same light and temperature conditions as the daily gas measurements. At the end of the incubation, the contents of each vial were filtered onto a pre-combusted glass fiber filter. The concentration and isotopic composition (δ15N) of particulate nitrogen for incubated and non-incubated (i.e., natural abundance) samples was measured using an elemental analyzer/isotope ratio mass spectrometer (Carlo-Erba EA NC2500 coupled with a ThermoFinnigan Delta Plus XP). For each of the described analyses, cell-specific rates were calculated based on the average of triplicate cell concentration measurements, obtained from cell samples preserved at 4 °C with Lugol’s iodine solution and quantified within a week using a Sedgwick-Rafter counting chamber (Electron Microscopy Sciences, cat. # 68050-52). All rate measurement data is provided in the Source Data file.Reporting summaryFurther information on research design is available in the Nature Research Reporting Summary linked to this article. 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

Chemicals and reagentsEOs of basil (Ocimum basilicum L.), bergamot (Citrus bergamia Risso & Poit), camphor [Cinnamomum camphora (L.) J. Presl.], cinnamon (Cinnamomum zeylanicum Blume), citronella [Cymbopogon nardus (L.) Rendle], clove (Eugenia caryophyllus Wight), eucalyptus (Eucalyptus globulus Labill.), jasmine (Jasminum officinale L.), lavender (Lavandula angustifolia Mill.), lemon grass [Cymbopogan citratus (DC.) Stapf], mentha (Mentha piperita L.), rosemary (Rosmarinus officinalis L.), patchouli (Pogostemon patchouli Benth), and wild turmeric (Curcuma aromatica Salisb.) were procured from Talent Technologies (Talent Technologies, Kanpur, India). Acetylcholinesterase (AChE) activity assay kit, Anti-OBP2A antibody, ELISA kits, 1,1-diphenyl-2-picrylhydrazyl (DPPH), radioimmunoprecipitation (RIPA) buffer and phosphate buffer saline (PBS) were purchased from Sigma Aldrich (Sigma Aldrich Chemical Co., St. Luis, USA). TRPV1 antibody was purchased from Santa Cruz (Santa Cruz, California, USA). 1-chloro-2,4-dinitrobenzene (CDNB) was purchased from Cayman (Cayman Chemical Company, Michigan, USA). Human normal lung cell line (L-132) was obtained from the National Centre for Cell Sciences (NCCS), Pune, India. High performance liquid chromatography (HPLC) grade acetone was purchased from Merck (Merck Pvt. Ltd., Mumbai, India). All other chemicals used were of the highest analytical grade available.Test insects5–7 days old adult female Ae. albopictus mosquitoes were housed at the laboratory insectary, Division of Pharmaceutical Technology, Defence Research Laboratory, Tezpur, Assam, India. Mosquitoes were reared by maintaining temperature at 27 ± 2 °C, relative humidity: 75 ± 5% RH and 14L:10D h of light–dark alternative cycles in standard-sized wooden cages (75 cm × 60 cm × 60 cm) with a sleeve opening on one side as described previously63. 10% sucrose solution ad libitum were provided for nourishment. Before testing, the mosquitoes were starved for 24 h.Screening of EOsDose response study was performed to evaluate the best oils among the fourteen EOs. This study was approved (approval number: 032/2021TMCH, 28/08/2018) by the Institutional Human Ethical Committee (IHEC), of the Tezpur Medical College & Hospital (TMCH), Tezpur, Assam, India, and all experiments were performed in accordance with relevant guidelines and regulations. Five volunteers are chosen, not allergic to mosquito bite and all volunteers provided written informed consent. A volunteer’s thigh was marked according to the door opening hole of the K&D module as described by Klun and Debboun64. It is made of Plexiglas and the base of the rectangular cage (26 cm × 5 cm × 5 cm) has six holes, each with rectangular 3 × 4 cm holes that are opened and closed by a sliding door (Supplementary Fig. S8: Provide the photograph of K&D module). The flexor region of the forearms of a human volunteer was outlined with four rectangular (3 cm × 4 cm) test areas. A volume of 25 µL of each concentration of the EOs in soybean oil (40, 4 and 0.4 µg/cm2) and 25 µL of the soybean oil (diluent) as control was applied to the marked areas. After air drying for 5 min, a K&D module with matching cut outs in its floor was placed over the treated areas, containing five nulliparous 5–7 days old female mosquitoes in each hole. The doors of the cells were opened and the number of mosquitoes biting in each cell was recorded within a 2 min exposure, after which the doors were closed. After completion of each observation, mosquitoes were freed by opening cells of the K&D module in a sleeved screened cage. For each test, fresh sets of mosquitoes are used. Five replications for each test were carried out. The efficacy of EOs were determined by the percentage repellency against mosquitoes, using the formula or Eq. (2) described by WHO46.$$% ;{text{repellency}} = frac{C – T}{C} times 100$$
(2)
where, C is the number of mosquitoes landing, or biting at the control area; T is the number of mosquitoes landing or biting at the treated area.Fourier transform-infra red spectroscopy (FT-IR)Study of chemical compatibility for each formulation ingredients are necessary. All formulation ingredients possess specific value of vibrational frequency and have varied functional groups in their chemical structures. For compatibility study, each EOs, excipients to be used in cream formulation, and their physical mixture was placed one by one over the sample plate of the FT-IR instrument (Bruker, ALPHA, Billerica, MA, USA). The covering probe was placed over the sample and IR spectra was obtained over a wavelength of 2.5–25 μm at room temperature. Functional groups possessed by each individual ingredient should be identical in their physical mixture which confirms their compatibility37.Thermogravimetric analysis (TGA)The thermal behaviour of citronella oil, clove oil, lemon grass oil, their mixture and EO-MRC were evaluated using a thermal analyser (TG 209 F1 Libra®, NETZSCH-Gerätebau GmbH, 95100 Selb, Germany). Approximately about 10 mg sample weight was placed in the crucible each time. Nitrogen was used as a shielding gas. Heating program was fixed as 30–600 °C at a rate of 10 °C/min.Formulation development and optimizationFor optimization, a 17-run, 3-factor, 3-level Box-Behnken design (BBD) was utilized. A second order polynomial model was constructed by quadratic response surface methodology (RSM) using Design-Expert software (Version 6.0.8, Stat-Ease Inc., USA). Total seventeen formulations were obtained using EO concentrations as dependent variables against complete protection time (CPT) as independent variable or response variable. Analysis of variance (ANOVA) was performed using the same software to obtain the most effective formulation.Preparation of creamPhase inversion temperature method was applied for the preparation of EO-based mosquito repellent cream (EO-MRC). About 50 g cream sample was prepared in order to get enough for performing the various qualitative and quantitative assay. The oil phase (phase B) was prepared by dissolving the oil soluble excipients, except phase A (mosquito repellent active ingredients) under mild heating at 200 rpm in a hot magnetic plate stirrer (Magnetic Stirrer IKA RCT basic) and heated to 65 °C. The aqueous phase was prepared by mixing various aqueous soluble ingredients (phase C) under gentle heating and stirring. Temperature of the aqueous phase was raised to 65 °C. Phase A was gently added to the oil phase at a stirring speed of 200 rpm and 55 ± 2 °C. The mixture was then emulsified by adding phase C slowly and kept for 1 h at a stirring rate of 800 rpm and 60 ± 2 °C. The formulated EO-MRC was then kept for natural cooling.Efficacy assessmentCPT of the developed cream (EO-MRC) formulation was carried out by arm in cage bioassay. 1 mL EO-MRC was applied to ≈ 600 cm2 area of the forearm skin between the wrist and elbow and 1 mL of the 12% N, N-di ethyl benzamide (DEBA) based marketed cream (DBMC) was compared on the other arm. Two mosquito cages (size: 40 × 40 × 40 cm) each containing 200–250 non-blood-fed female Ae. Albopictus were used. One cage is designated for testing the EO-MRC and the other for the positive control (DBMC). During testing, hands were protected by surgical gloves for which the mosquitoes cannot bite while the volunteer avoids movement of the arm. EO-MRC and DBMC treated arms were exposed for 3 min at 30 min intervals to determine landing and/or probing activity. A single landing or probing of mosquito within a 3 min test interval concludes the test. CPT was calculated as the time (min) required for the first mosquito landing or probing after repellent application to the treated area. The median CPT and confidence intervals were estimated from the Kaplan–Meier Survival Function46.Efficacy was correlated with DEBA based marketed cream (DBMC). The inclusion of the specific commercial product DBMC is for comparison and does not constitute any recommendations.CharacterizationGas chromatography-mass spectroscopy (GC–MS)Qualitative studyDifferent chemical components in fourteen EOs and the selected blend were identified by a GC–MS system of Agilent Technologies (5301 Stevens Creek Blvd. Santa Clara, CA 95051, United States). Test sample concentration of 500 μg/mL was prepared in GC grade acetone. A sample volume of 1 μL was introduced into the injector held at 250 °C. Oven temperature of 40–300 °C was programmed at 20 °C/min. Helium was used as carrier gas at flow rate 1 mL/min. The injector and detector temperature were set at 250 °C and 230 °C (quad) and 150 °C (core) respectively37. Standard C7–C30 saturated alkanes were purchased from Sigma Aldrich Chemicals Co., St. Louis, USA. Retention indices (RI) of the identified components were determined for identification of the detected components.% Assay by GC–MS studyCalibration samples of eugenol and citronellol were prepared by dissolving an appropriate amount in GC grade acetone to get concentrations of 62.5 μg/mL, 125 μg/mL, 250 μg/mL and 500 μg/mL. Test samples of EO-MRC, clove oil and citronella oil were prepared by dissolving a required amount in acetone to quantify the EO components in the final formulation. A sample volume of 1 μL was introduced into the injector as described in ‘Qualitative study’ section.Physicochemical parametersPhysical parameters of the EO-MRC and placebo formulations were determined in order to establish aesthetic compliance and consumer acceptability. To determine the viscosity, a programmable viscometer was used (Model: DV2T, Ametek Brookfield, Middleboro, MA, USA); combined with software Rheo3000, version 1.2.2019.1 [R]. Sample volume was fixed at 30 g and viscosities were determined at 10 rpm for 40 s at room temperature using a T-Bar spindle (B-92) (Helipath spindle set, Brookfield Engineering Labs. Inc). Density was determined by using a pycnometer. pH of EO-MRC was checked by using digital pH meter (Labman Scientific instruments, Tamil Nadu, India).Spread ability of EO-MRC was determined as per the method reported earlier by Sabale65. In brief, 1 g of EO-MRC was placed on 1 cm2 pre-marked circular area on the glass slide (7.5 cm × 2.5 cm). EO-MRC was compressed using another glass slide placed from edge to centre of primary slide. 200 g of commercial weight was placed on the set up and allowed the gel to spread for the period of 1 min. The spread diameter was calculated with the aid of graph paper and spread ability was evaluated using formula expressed as Eq. (3):$$mathrm{Spread, ability}=mathrm{m}times frac{mathrm{l}}{mathrm{t}}$$
(3)
where, m is the commercial weight placed on the setup; l is the length of cream spread; and t is the time.Safety assessmentCytotoxicity by MTT assayThe reduction of tetrazolium salts is now widely accepted as a reliable way to examine cell proliferation. The yellow tetrazolium MTT (3-(4,5-dimethylthiazolyl-2)-2,5-diphenyltetrazolium bromide) is reduced by metabolically active cells, in part by the action of dehydrogenase enzymes, to generate reducing equivalents such as NADH and NADPH. With the help of spectrophotometric means, the resulting intracellular purple formazan can be quantified. The assay measures the cell proliferation rate and conversely, when metabolic events cause apoptosis or necrosis, the reduction in cell viability66.Cells cultured in T-25 flasks were trypsinized and aspirated into a 5 mL centrifuge tube. Cell pellet was obtained by centrifugation at 3000 rpm. The cell count was adjusted, using DMEM HG medium, such that 200 μL of suspension contained approximately 10,000 cells. To each well of the 96 well microtiter plate, 200 μL of the cell suspension was added and the plate was incubated at 37 ℃ and 5% CO2 atmosphere for 24 h. After 24 h, the spent medium was aspirated. 200 μL of different test concentrations viz. 62 µg/mL, 125 µg/mL, 250 µg/mL, 500 µg/mL, and 1000 µg/mL, of EO-MRC were added to the respective wells. The plate was then incubated at 37 °C and 5% CO2 atmosphere for 24 h. The plate was removed from the incubator and the drug containing media was aspirated. 200 μL of medium containing10% MTT reagent was then added to each well to get a final concentration of 0.5 mg/mL and the plate was incubated at 37 ℃ and 5% CO2 atmosphere for 3 h. Without disturbing the crystals formed in the wells, culture medium was completely removed. 100 μL of solubilisation solution (DMSO) was added to each well and the plate was then gently shake in a rocking shaker (ROCKYMAX™, Tarsons, Kolkata, India) to solubilize the formed formazan. The absorbance was measured at a wavelength of 570 nm and also at 630 nm using a microplate reader. The percentage growth inhibition was calculated and concentration of EO-MRC needed to inhibit cell growth by 50% (IC50) was generated from the dose–response curve for the cell line.Animals and ethics statementAll experimenting protocols using animal were performed according to the “Principles of Laboratory Animal care” (NIH publication 85–23, revised 1985) and approved by the Institutional Animal Ethical Committee (IAEC) of Defence Research Laboratory (DRL), Tezpur, Assam, India (approval no. CPCSEA/DRL/Protocol no. 3, 20/06/2018). All studies involving animals are reported in accordance with the ARRIVE guidelines for reporting experiments involving animals67. All efforts were made during the study period to minimize the suffering of animals and to reduce the number of animals used.5–8 weeks old, about 210–250 g of male healthy adult Wistar rats (Rattus norvegicus) and young and healthy New Zealand albino rabbits (Oryctolagus cuniculus) were obtained from the institutional animal housing facility and allowed to acclimatize for 7 days prior to the study. Standard food and purified water ad libitum were provided in clean and hygienic condition at 22–25 ℃, 40–70% RH with 12 h light–dark cycles.Acute dermal irritation studyAcute dermal irritation study was conducted on healthy New Zealand albino rabbits following the OECD test guidelines 40468. Approximately 24 h before the test, fur was removed from the dorsal area of the trunk. 0.5 g EO-MRC, was directly applied to the skin and after 4 h exposure period, residual EO-MRC was removed by using water without disturbing the integrity of the epidermis and examined for signs of erythema and oedema, at 60 min, and then at 24 h, 48 h and 72 h after EO-MRC removal. Dermal reactions are graded and recorded according to the grades in the Table 8. As per the method described by Banerjee et al.69; primary irritation index (PII) was calculated. Further, we have followed the Draize method of classification for PII scoring as non-irritant (if PII  More

1.Ceballos, G., Ehrlich, P. R. & Raven, P. H. Vertebrates on the brink as indicators of biological annihilation and the sixth mass extinction. Proc. Natl Acad. Sci. USA 117, 13596–13602 (2020).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
2.Committee, T. I. R. L. The IUCN Red List of Threatened Species – Strategic Plan 2017-2020. Tech. Rep., IUCN (2017).3.Witmer, G. W. Wildlife population monitoring: some practical considerations. Wild. Res. 32, 259–263 (2005).
Google Scholar 
4.McEvoy, J. F., Hall, G. P. & McDonald, P. G. Evaluation of unmanned aerial vehicle shape, flight path and camera type for waterfowl surveys: disturbance effects and species recognition. PeerJ 4, e1831 (2016).PubMed 
PubMed Central 

Google Scholar 
5.Burghardt, G. M. et al. Perspectives–minimizing observer bias in behavioral studies: a review and recommendations. Ethology 118, 511–517 (2012).
Google Scholar 
6.Giese, M. Effects of human activity on Adelie penguin Pygoscelis adeliae breeding success. Biol. Conserv. 75, 157–164 (1996).
Google Scholar 
7.Köndgen, S. et al. Pandemic human viruses cause decline of endangered great apes. Curr. Biol. 18, 260–264 (2008).PubMed 

Google Scholar 
8.Weissensteiner, M. H., Poelstra, J. W. & Wolf, J. B. W. Low-budget ready-to-fly unmanned aerial vehicles: an effective tool for evaluating the nesting status of canopy-breeding bird species. J. Avian Biol. 46, 425–430 (2015).
Google Scholar 
9.Sasse, D. B. Job-related mortality of wildlife workers in the united states, 1937–2000. Wildl. Soc. Bull. 31, 1015–1020 (2003).10.Kays, R., Crofoot, M. C., Jetz, W. & Wikelski, M. Terrestrial animal tracking as an eye on life and planet. Science 348, aaa2478 (2015).11.Altmann, J. Observational study of behavior: sampling methods. Behaviour 49, 227–266 (1974).CAS 
PubMed 

Google Scholar 
12.Hodgson, J. C. et al. Drones count wildlife more accurately and precisely than humans. Methods Ecol. Evolution 9, 1160–1167 (2018).
Google Scholar 
13.Betke, M. et al. Thermal imaging reveals significantly smaller Brazilian free-tailed bat colonies than previously estimated. J. Mammal. 89, 18–24 (2008).
Google Scholar 
14.Rollinson, C. R. et al. Working across space and time: nonstationarity in ecological research and application. Front. Ecol. Environ. 19, 66–72 (2021).
Google Scholar 
15.Junker, J. et al. A severe lack of evidence limits effective conservation of the world’s primates. BioScience 70, 794–803 (2020).PubMed 
PubMed Central 

Google Scholar 
16.Sherman, J., Ancrenaz, M. & Meijaard, E. Shifting apes: Conservation and welfare outcomes of Bornean orangutan rescue and release in Kalimantan, Indonesia. J. Nat. Conserv. 55, 125807 (2020).
Google Scholar 
17.O’Donoghue, P. & Rutz, C. Real-time anti-poaching tags could help prevent imminent species extinctions. J. Appl. Ecol. 53, 5–10 (2016).PubMed 

Google Scholar 
18.Lahoz-Monfort, J. J. & Magrath, M. J. L. A comprehensive overview of technologies for species and habitat monitoring and conservation. BioScience biab073. https://academic.oup.com/bioscience/advance-article/doi/10.1093/biosci/biab073/6322306 (2021).19.Gottschalk, T., Huettmann, F. & Ehlers, M. Thirty years of analysing and modelling avian habitat relationships using satellite imagery data: a review. Int. J. Remote Sens. 26, 2631–2656 (2005).
Google Scholar 
20.Steenweg, R. et al. Scaling-up camera traps: monitoring the planet’s biodiversity with networks of remote sensors. Front. Ecol. Environ. 15, 26–34 (2017).
Google Scholar 
21.Hausmann, A. et al. Social media data can be used to understand tourists’ preferences for nature-based experiences in protected areas. Conserv. Lett. 11, e12343 (2018).
Google Scholar 
22.Sugai, L. S. M., Silva, T. S. F., Ribeiro, J. W. & Llusia, D. Terrestrial passive acoustic monitoring: review and perspectives. BioScience 69, 15–25 (2018).
Google Scholar 
23.Wikelski, M. et al. Going wild: what a global small-animal tracking system could do for experimental biologists. J. Exp. Biol. 210, 181–186 (2007).PubMed 

Google Scholar 
24.Belyaev, M. Y. et al. Development of technology for monitoring animal migration on Earth using scientific equipment on the ISS RS. in 2020 27th Saint Petersburg International Conference on Integrated Navigation Systems (ICINS), 1–7 (IEEE, 2020).25.Harel, R., Loftus, J. C. & Crofoot, M. C. Locomotor compromises maintain group cohesion in baboon troops on the move. Proc. R. Soc. B 288, 20210839 (2021).PubMed 
PubMed Central 

Google Scholar 
26.Farley, S. S., Dawson, A., Goring, S. J. & Williams, J. W. Situating ecology as a big-data science: current advances, challenges, and solutions. BioScience 68, 563–576 (2018).
Google Scholar 
27.Lasky, M. et al. Candid critters: Challenges and solutions in a large-scale citizen science camera trap project. Citizen Science: Theory and Practice 6, https://doi.org/10.5334/cstp.343 (2021).28.Hastie, T., Tibshirani, R. & Friedman, J. The Elements of Statistical Learning: Data Mining, Inference, and Prediction (Springer, 2001).29.Christin, S., Hervet, É. & Lecomte, N. Applications for deep learning in ecology. Methods Ecol. Evolution 10, 1632–1644 (2019).
Google Scholar 
30.Kwok, R. Ai empowers conservation biology. Nature 567, 133–135 (2019).ADS 
CAS 
PubMed 

Google Scholar 
31.Kwok, R. Deep learning powers a motion-tracking revolution. Nature 574, 137–139 (2019).ADS 
CAS 
PubMed 

Google Scholar 
32.LeCun, Y., Bengio, Y. & Hinton, G. Deep learning. Nature 521, 436–444 (2015).ADS 
CAS 
PubMed 

Google Scholar 
33.Pichler, M., Boreux, V., Klein, A.-M., Schleuning, M. & Hartig, F. Machine learning algorithms to infer trait-matching and predict species interactions in ecological networks. Methods Ecol. Evolution 11, 281–293 (2020).
Google Scholar 
34.Knudby, A., LeDrew, E. & Brenning, A. Predictive mapping of reef fish species richness, diversity and biomass in Zanzibar using IKONOS imagery and machine-learning techniques. Remote Sens. Environ. 114, 1230–1241 (2010).ADS 

Google Scholar 
35.Rey, N., Volpi, M., Joost, S. & Tuia, D. Detecting animals in African savanna with UAVs and the crowds. Remote Sens. Environ. 200, 341–351 (2017).ADS 

Google Scholar 
36.Beery, S., Morris, D. & Yang, S. Efficient pipeline for camera trap image review. in Proceedings of the Workshop Data Mining and AI for Conservation, Conference for Knowledge Discovery and Data Mining (2019).37.Kellenberger, B., Marcos, D. & Tuia, D. When a few clicks make all the difference: improving weakly-supervised wildlife detection in UAV images. in Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition Workshops (2019).38.Schofield, D. et al. Chimpanzee face recognition from videos in the wild using deep learning. Sci. Adv. 5, eaaw0736 (2019).ADS 
PubMed 
PubMed Central 

Google Scholar 
39.Ditria, E. M. et al. Automating the analysis of fish abundance using object detection: optimizing animal ecology with deep learning. Front. Mar. Sci. 7, 429 (2020).
Google Scholar 
40.Kellenberger, B., Veen, T., Folmer, E. & Tuia, D. 21 000 birds in 4.5 h: efficient large-scale seabird detection with machine learning. Remote Sens. Ecol. Conserv. 7, 445–460 (2021).
Google Scholar 
41.Ahumada, J. A. et al. Wildlife insights: a platform to maximize the potential of camera trap and other passive sensor wildlife data for the planet. Environ. Conserv. 47, 1–6 (2020).MathSciNet 

Google Scholar 
42.Eikelboom, J. A. J. et al. Improving the precision and accuracy of animal population estimates with aerial image object detection. Methods Ecol. Evolution 10, 1875–1887 (2019).
Google Scholar 
43.Weinstein, B. G. A computer vision for animal ecology. J. Anim. Ecol. 87, 533–545 (2018).PubMed 

Google Scholar 
44.Valletta, J. J., Torney, C., Kings, M., Thornton, A. & Madden, J. Applications of machine learning in animal behaviour studies. Anim. Behav. 124, 203–220 (2017).
Google Scholar 
45.Peters, D. P. C. et al. Harnessing the power of big data: infusing the scientific method with machine learning to transform ecology. Ecosphere 5, art67 (2014).
Google Scholar 
46.Yu, Q. et al. Study becomes insight: ecological learning from machine learning. Methods Ecol. Evol. 12, 2117–2128 (2021).47.Lucas, T. C. D. A translucent box: interpretable machine learning in ecology. Ecol. Monogr. 90, https://doi.org/10.1002/ecm.1422 (2020).48.Reichstein, M. et al. Deep learning and process understanding for data-driven Earth system science. Nature 566, 195–204 (2019).ADS 
CAS 
PubMed 

Google Scholar 
49.Camps-Valls, G., Tuia, D., Zhu, X. X. & Reichstein, M. Deep Learning for the Earth Sciences: A Comprehensive Approach to Remote Sensing, Climate Science and Geosciences (Wiley & Sons, 2021).50.Karpatne, A. et al. Theory-guided data science: A new paradigm for scientific discovery from data. IEEE Trans. Knowl. Data Eng. 29, 2318–2331 (2017).
Google Scholar 
51.Oliver, R. Y., Meyer, C., Ranipeta, A., Winner, K. & Jetz, W. Global and national trends, gaps, and opportunities in documenting and monitoring species distributions. PLoS Biol 19, e3001336 https://doi.org/10.1371/journal.pbio.3001336 (2021).52.Beery, S., Wu, G., Rathod, V., Votel, R. & Huang, J. Context R-CNN: long term temporal context for per-camera object detection. in 2020 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 13075–13085 (2020).53.Norouzzadeh, M. S. et al. Automatically identifying, counting, and describing wild animals in camera-trap images with deep learning. Proc. Natl Acad. Sci. USA 115, E5716–E5725 (2018).CAS 
PubMed 
PubMed Central 

Google Scholar 
54.Schneider, S., Taylor, G. W., Linquist, S. & Kremer, S. C. Past, present and future approaches using computer vision for animal re-identification from camera trap data. Methods Ecol. Evolution 10, 461–470 (2019).
Google Scholar 
55.Beery, S., Van Horn, G. & Perona, P. Recognition in terra incognita. in 2018 European Conference on Computer Vision (ECCV), 456–473 (2018).56.Sugai, L. S. M., Silva, T. S. F., Ribeiro Jr, J. W. & Llusia, D. Terrestrial passive acoustic monitoring: review and perspectives. BioScience 69, 15–25 (2019).
Google Scholar 
57.Wrege, P. H., Rowland, E. D., Keen, S. & Shiu, Y. Acoustic monitoring for conservation in tropical forests: examples from forest elephants. Methods Ecol. Evolution 8, 1292–1301 (2017).
Google Scholar 
58.Desjonquères, C., Gifford, T. & Linke, S. Passive acoustic monitoring as a potential tool to survey animal and ecosystem processes in freshwater environments. Freshw. Biol. 65, 7–19 (2020).
Google Scholar 
59.Davis, G. E. et al. Long-term passive acoustic recordings track the changing distribution of North Atlantic right whales (eubalaena glacialis) from 2004 to 2014. Sci. Rep. 7, 1–12 (2017).
Google Scholar 
60.Wood, C. M. et al. Detecting small changes in populations at landscape scales: a bioacoustic site-occupancy framework. Ecol. Indic. 98, 492–507 (2019).
Google Scholar 
61.Kahl, S., Wood, C. M., Eibl, M. & Klinck, H. Birdnet: a deep learning solution for avian diversity monitoring. Ecol. Inform. 61, 101236 (2021).
Google Scholar 
62.Stowell, D., Wood, M. D., Pamuła, H., Stylianou, Y. & Glotin, H. Automatic acoustic detection of birds through deep learning: the first bird audio detection challenge. Methods Ecol. Evolution 10, 368–380 (2019).
Google Scholar 
63.Ford, J. K. B. in Encyclopedia of Marine Mammals 253–254 (Elsevier, 2018).64.Hughey, L. F., Hein, A. M., Strandburg-Peshkin, A. & Jensen, F. H. Challenges and solutions for studying collective animal behaviour in the wild. Philos. Trans. R. Soc. B: Biol. Sci. 373, 20170005 (2018).
Google Scholar 
65.Williams, H. J. et al. Optimizing the use of biologgers for movement ecology research. J. Anim. Ecol. 89, 186–206 (2020).PubMed 

Google Scholar 
66.Korpela, J. et al. Machine learning enables improved runtime and precision for bio-loggers on seabirds. Commun. Biol. 3, 1–9 (2020).
Google Scholar 
67.Yu, H. An evaluation of machine learning classifiers for next-generation, continuous-ethogram smart trackers. Mov. Ecol. 9, 14 (2021).
Google Scholar 
68.Browning, E. et al. Predicting animal behaviour using deep learning: GPS data alone accurately predict diving in seabirds. Methods Ecol. Evolution 9, 681–692 (2018).
Google Scholar 
69.Liu, Z. Y.-C. et al. Deep learning accurately predicts white shark locomotor activity from depth data. Anim. Biotelemetry 7, 1–13 (2019).
Google Scholar 
70.Wang, G. Machine learning for inferring animal behavior from location and movement data. Ecol. Inform. 49, 69–76 (2019).
Google Scholar 
71.Wijeyakulasuriya, D. A., Eisenhauer, E. W., Shaby, B. A. & Hanks, E. M. Machine learning for modeling animal movement. PLoS ONE 30, e0235750 (2020).72.Linchant, J., Lisein, J., Semeki, J., Lejeune, P. & Vermeulen, C. Are unmanned aircraft systems (UASs) the future of wildlife monitoring? A review of accomplishments and challenges. Mammal. Rev. 45, 239–252 (2015).
Google Scholar 
73.Hodgson, J. C., Baylis, S. M., Mott, R., Herrod, A. & Clarke, R. H. Precision wildlife monitoring using unmanned aerial vehicles. Sci. Rep. 6, 1–7 (2016).
Google Scholar 
74.Mathis, A. et al. DeepLabCut: markerless pose estimation of user-defined body parts with deep learning. Nat. Neurosci. 21, 1281–1289 (2018).CAS 
PubMed 

Google Scholar 
75.Graving, J. M. et al. DeepPoseKit, a software toolkit for fast and robust animal pose estimation using deep learning. Elife 8, e47994 (2019).CAS 
PubMed 
PubMed Central 

Google Scholar 
76.Mathis, A., Schneider, S., Lauer, J. & Mathis, M. W. A primer on motion capture with deep learning: principles, pitfalls, and perspectives. Neuron 108, 44–65 (2020).CAS 
PubMed 

Google Scholar 
77.Kellenberger, B., Marcos, D. & Tuia, D. Detecting mammals in UAV images: best practices to address a substantially imbalanced dataset with deep learning. Remote Sens. Environ. 216, 139–153 (2018).ADS 

Google Scholar 
78.Kellenberger, B., Veen., T., Folmer, E. & Tuia, D. 21,000 birds in 4.5 hours: efficient large-scale seabird detection with machine learning. Remote Sens. Ecol. Conserv. https://doi.org/10.1002/rse2.200 (2021).79.Andrew, W., Greatwood, C. & Burghardt, T. Aerial animal biometrics: individual Friesian cattle recovery and visual identification via an autonomous UAV with onboard deep inference. in International Conference on Intelligent Robots and Systems (IROS) (2019).80.Schroeder, N. M., Panebianco, A., Gonzalez Musso, R. & Carmanchahi, P. An experimental approach to evaluate the potential of drones in terrestrial mammal research: a gregarious ungulate as a study model. R. Soc. open Sci. 7, 191482 (2020).ADS 
PubMed 
PubMed Central 

Google Scholar 
81.Bennitt, E., Bartlam-Brooks, H. L. A., Hubel, T. Y. & Wilson, A. M. Terrestrial mammalian wildlife responses to Unmanned Aerial Systems approaches. Sci. Rep. 9, 1–10 (2019).CAS 

Google Scholar 
82.Deneu, B., Servajean, M., Botella, C. & Joly, A. Evaluation of deep species distribution models using environment and co-occurrences. in International Conference of the Cross-Language Evaluation Forum for European Languages, 213–225 (Springer, 2019).83.Zhu, X. et al. Deep learning in remote sensing: A comprehensive review and list of resources. IEEE Geosci. Remote Sens. Mag. 5, 8–36 (2017).
Google Scholar 
84.Guirado, E., Tabik, S., Rivas, M. L., Alcaraz-Segura, D. & Herrera, F. Whale counting in satellite and aerial images with deep learning. Sci. Rep. 9, 1–12 (2019).CAS 

Google Scholar 
85.Duporge, I., Isupova, O., Reece, S., Macdonald, D. W. & Wang, T. Using very-high-resolution satellite imagery and deep learning to detect and count African elephants in heterogeneous landscapes. Remote Sens. Ecol. Conserv. https://doi.org/10.1002/rse2.195 (2020).86.Fretwell, P. T. & Trathan, P. N. Discovery of new colonies by Sentinel2 reveals good and bad news for emperor penguins. Remote Sens. Ecol. Conserv. https://doi.org/10.1002/rse2.176 (2020).87.Brodrick, P. G., Davies, A. B. & Asner, G. P. Uncovering ecological patterns with convolutional neural networks. Trends Ecol. Evolution 34, 734–745 (2019).
Google Scholar 
88.Audebert, N., Le Saux, B. & Lefèvre, S. Deep learning for classification of hyperspectral data: a comparative review. IEEE Geosci. Remote Sens. Mag. 7, 159–173 (2019).
Google Scholar 
89.McKinley, D. C. et al. Citizen science can improve conservation science, natural resource management, and environmental protection. Biol. Conserv. 208, 15–28 (2017).
Google Scholar 
90.Wäldchen, J. & Mäder, P. Machine learning for image based species identification. Methods Ecol. Evolution 9, 2216–2225 (2018).MATH 

Google Scholar 
91.Torney, C. J. et al. A comparison of deep learning and citizen science techniques for counting wildlife in aerial survey images. Methods Ecol. Evolution 10, 779–787 (2019).
Google Scholar 
92.Parham, J., Crall, J., Stewart, C., Berger-Wolf, T. & Rubenstein, D. I. Animal population censusing at scale with citizen science and photographic identification. in AAAI Spring Symposium-Technical Report (2017).93.Kühl, H. S. & Burghardt, T. Animal biometrics: quantifying and detecting phenotypic appearance. Trends Ecol. Evolution 28, 432–441 (2013).
Google Scholar 
94.Yu, X. et al. Automated identification of animal species in camera trap images. EURASIP J. Image Video Process. 2013, 1–10 (2013).ADS 

Google Scholar 
95.Mac Aodha, O. et al. Bat detective–deep learning tools for bat acoustic signal detection. PLoS Computat. Biol. 14, e1005995 (2018).
Google Scholar 
96.Schindler, F. & Steinhage, V. Identification of animals and recognition of their actions in wildlife videos using deep learning techniques. Ecol. Inform. 61, 101215 (2021).97.Avise, J. C. Molecular Markers, Natural History and Evolution (Springer Science & Business Media, 2012).98.Vidal, M., Wolf, N., Rosenberg, B., Harris, B. P. & Mathis, A. Perspectives on Individual Animal Identification from Biology and Computer Vision. Integr. Comp. Biol. 61, 900–916 https://doi.org/10.1093/icb/icab107 (2021).Article 
PubMed 
PubMed Central 

Google Scholar 
99.Berger-Wolf, T. Y. et al. Wildbook: crowdsourcing, computer vision, and data science for conservation. Preprint at https://arxiv.org/abs/1710.08880 (2017).100.Parham, J. et al. An animal detection pipeline for identification. in IEEE Winter Conference on Applications of Computer Vision (WACV), 1075–1083 (IEEE, 2018).101.Weideman, H. et al. Extracting identifying contours for African elephants and humpback whales using a learned appearance model. in Proceedings of the IEEE/CVF Winter Conference on Applications of Computer Vision (WACV) (2020).102.Brust, C.-A. et al. Towards automated visual monitoring of individual gorillas in the wild. in 2017 IEEE International Conference on Computer Vision Workshops (ICCVW), 2820–2830 (2017).103.Li, S., Li, J., Tang, H., Qian, R. & Lin, W. ATRW: a benchmark for Amur tiger re-identification in the wild. in 2020 ACM International Conference on Multimedia, 2590–2598 (2020).104.Bendale, A. & Boult, T. E. Towards open set deep networks. in 2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 1563–1572 (2016).105.Mathis, M. W. & Mathis, A. Deep learning tools for the measurement of animal behavior in neuroscience. Curr. Opin. Neurobiol. 60, 1–11 (2020).CAS 
PubMed 

Google Scholar 
106.Sanakoyeu, A., Khalidov, V., McCarthy, M. S., Vedaldi, A. & Neverova, N. Transferring dense pose to proximal animal classes. in 2020 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 5233–5242 (2020).107.Zuffi, S., Kanazawa, A., Jacobs, D. W. & Black, M. J. 3D menagerie: modeling the 3D shape and pose of animals. in 2017 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 6365–6373 (2017).108.Biggs, B., Roddick, T., Fitzgibbon, A. & Cipolla, R. Creatures great and smal: recovering the shape and motion of animals from video. in 2018 Asian Conference on Computer Vision (ACCV), 3–19 (Springer, 2018).109.Biggs, B., Boyne, O., Charles, J., Fitzgibbon, A. & Cipolla, R. Who left the dogs out? 3D animal reconstruction with expectation maximization in the loop. in 2020 European Conference on Computer Vision (ECCV), 195–211 (Springer, 2020).110.Zuffi, S., Kanazawa, A., Berger-Wolf, T. & Black, M. J. Three-D safari: learning to estimate zebra pose, shape, and texture from images” in the wild”. in 2019 IEEE International Conference on Computer Vision (ICCV), 5359–5368 (2019).111.Wang, Y., Kolotouros, N., Daniilidis, K. & Badger, M. Birds of a feather: capturing avian shape models from images. in 2021 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 14739–14749 (2021).112.Haalck, L., Mangan, M., Webb, B. & Risse, B. Towards image-based animal tracking in natural environments using a freely moving camera. J. Neurosci. methods 330, 108455 (2020).PubMed 

Google Scholar 
113.Pettorelli, N. et al. Satellite remote sensing for applied ecologists: opportunities and challenges. J. Appl. Ecol. 51, 839–848 (2014).
Google Scholar 
114.Davies, A. B., Tambling, C. J., Kerley, G. I. H. & Asner, G. P. Effects of vegetation structure on the location of lion kill sites in African thicket. PLoS ONE 11, e0149098 (2016).PubMed 
PubMed Central 

Google Scholar 
115.Froidevaux, J. S. P., Zellweger, F., Bollmann, K., Jones, G. & Obrist, M. K. From field surveys to LiDAR: shining a light on how bats respond to forest structure. Remote Sens. Environ. 175, 242–250 (2016).ADS 

Google Scholar 
116.Risse, B., Mangan, M., Stürzl, W. & Webb, B. Software to convert terrestrial LiDAR scans of natural environments into photorealistic meshes. Environ. Model. Softw. 99, 88–100 (2018).
Google Scholar 
117.Haalck, L. & Risse, B. Embedded dense camera trajectories in multi-video image mosaics by geodesic interpolation-based reintegration. in 2021 IEEE Winter Conference on Applications of Computer Vision (WACV), 1849–1858 (2021).118.Schonberger, J. L. & Frahm, J.-M. Structure-from-motion revisited. in 2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 4104–4113 (2016).119.Mur-Artal, R. & Tardós, J. D. ORB-SLAM2: an open-source SLAM system for monocular, stereo, and RGB-D cameras. IEEE Trans. Robot. 33, 1255–1262 (2017).
Google Scholar 
120.Kuppala, K., Banda, S. & Barige, T. R. An overview of deep learning methods for image registration with focus on feature-based approaches. Int. J. Image Data Fusion 11, 113–135 (2020).ADS 

Google Scholar 
121.Lisein, J., Linchant, J., Lejeune, P., Bouché, P. & Vermeulen, C. Aerial surveys using an unmanned aerial system (UAS): comparison of different methods for estimating the surface area of sampling strips. Tropical Conserv. Sci. 6, 506–520 (2013).
Google Scholar 
122.Wu, C. Critical configurations for radial distortion self-calibration. in 2014 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 25–32 (2014).123.Ferrer, J., Elibol, A., Delaunoy, O., Gracias, N. & Garcia, R. Large-area photo-mosaics using global alignment and navigation data. in Mts/IEEE Oceans Conference, 1–9 (2007).124.Guisan, A. & Zimmermann, N. E. Predictive habitat distribution models in ecology. Ecol. Model. 135, 147–186 (2000).
Google Scholar 
125.Lehmann, A., Overton, J. M. & Austin, M. P. Regression models for spatial prediction: their role for biodiversity and conservation. Biodivers. Conserv. 11, 2085–2092 (2002).126.Breiman, L. Random forests. Mach. Learn. 45, 5–32 (2001).MATH 

Google Scholar 
127.Parravicini, V. et al. Global patterns and predictors of tropical reef fish species richness. Ecography 36, 1254–1262 (2013).
Google Scholar 
128.Smoliński, S. & Radtke, K. Spatial prediction of demersal fish diversity in the baltic sea: comparison of machine learning and regression-based techniques. ICES J. Mar. Sci. 74, 102–111 (2017).
Google Scholar 
129.Čandek, K., Čandek, U. P. & Kuntner, M. Machine learning approaches identify male body size as the most accurate predictor of species richness. BMC Biol. 18, 1–16 (2020).
Google Scholar 
130.Baltensperger, A. P. & Huettmann, F. Predictive spatial niche and biodiversity hotspot models for small mammal communities in Alaska: applying machine-learning to conservation planning. Landscape Ecol. 30, 681–697 (2015).131.Faisal, A., Dondelinger, F., Husmeier, D. & Beale, C. M. Inferring species interaction networks from species abundance data: a comparative evaluation of various statistical and machine learning methods. Ecol. Inform. 5, 451–464 (2010).
Google Scholar 
132.Van Horn, G. et al. The inaturalist species classification and detection dataset. in 2018 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 8769–8778 (2018).133.Copas, K. et al. Training machines to improve species identification using GBIF-mediated datasets. in AGU Fall Meeting Abstracts, Vol. 2019, IN53C–0758 (2019).134.Lennox, R. J. et al. A novel framework to protect animal data in a world of ecosurveillance. BioScience 70, 468–476 (2020).
Google Scholar 
135.Strubell, E., Ganesh, A. & McCallum, A. Energy and policy considerations for deep learning in NLP. in Proceedings of the 57th Annual Meeting of the Association for Computational Linguistics, 3645–3650 (2019).136.Samek, W., Montavon, G., Vedaldi, A., Hansen, L. K. & Müller, K.-R. Explainable AI: Interpreting, Explaining and Visualizing Deep Learning, Vol. 11700 (Springer Nature, 2019).137.Swanson, A. et al. Snapshot Serengeti, high-frequency annotated camera trap images of 40 mammalian species in an African savanna. Sci. data 2, 1–14 (2015).
Google Scholar 
138.de Lutio, R. et al. Digital taxonomist: identifying plant species in community scientists’ photographs. ISPRS J. Photogramm. Remote Sens. 182, 112–121 (2021).139.Mac Aodha, O., Cole, E. & Perona, P. Presence-only geographical priors for fine-grained image classification. in Proceedings of the IEEE/CVF International Conference on Computer Vision, 9596–9606 (2019).140.Gurumurthy, S. et al. Exploiting Data and Human Knowledge for Predicting Wildlife Poaching. in Proceedings of the 1st ACM SIGCAS Conference on Computing and Sustainable Societies, 1–8, https://doi.org/10.1145/3209811.3209879 (ACM, 2018).141.Datta, S., Anderson, D., Branson, K., Perona, P. & Leifer, A. Computational neuroethology: a call to action. Neuron 104, 11–24 (2019).CAS 
PubMed 
PubMed Central 

Google Scholar 
142.Joska, D. et al. AcinoSet: a 3D pose estimation dataset and baseline models for Cheetahs in the wild. 2021 IEEE International Conference on Robotics and Automation (ICRA) Preprint at https://arxiv.org/abs/2103.13282 (IEEE, Xi’an, China, 2021).143.Chen, Q. & Koltun, V. Photographic image synthesis with cascaded refinement networks. in 2019 IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 1511–1520 (2017).144.Lee, J., Hwangbo, J., Wellhausen, L., Koltun, V. & Hutter, M. Learning quadrupedal locomotion over challenging terrain. Sci. Robot. 5, eabc5986 (2020).145.Botella, C., Joly, A., Bonnet, P., Munoz, F. & Monestiez, P. Jointly estimating spatial sampling effort and habitat suitability for multiple species from opportunistic presence-only data. Methods Ecol. Evolution 12, 933–945 (2021).
Google Scholar 
146.Beery, S., Cole, E., Parker, J., Perona, P. & Winner, K. Species distribution modeling for machine learning practitioners: a review. in Proceedings of the 4th ACM SIGCAS Conference on Computing and Sustainable Societies (2021).147.Arzoumanian, Z., Holmberg, J. & Norman, B. An astronomical pattern-matching algorithm for computer-aided identification of whale sharks Rhincodon typus. J. Appl. Ecol. 42, 999–1011 (2005).
Google Scholar 
148.de Knegt, H. J., Eikelboom, J. A. J., van Langevelde, F., Spruyt, W. F. & Prins, H. H. T. Timely poacher detection and localization using sentinel animal movement. Sci. Rep. 11, 1–11 (2021).
Google Scholar 
149.Walter, T. & Couzin, I. D. TRex, a fast multi-animal tracking system with markerless identification, and 2D estimation of posture and visual fields. eLife 10, e64000 (2021).PubMed 
PubMed Central 

Google Scholar 
150.Kellenberger, B., Tuia, D. & Morris, D. AIDE: accelerating image-based ecological surveys with interactive machine learning. Methods Ecol. Evolution 11, 1716–1727 (2020).
Google Scholar 
151.Settles, B. Active learning. Synth. lectures Artif. Intell. Mach. Learn. 6, 1–114 (2012).MathSciNet 
MATH 

Google Scholar 
152.Ofli, F. et al. Combining human computing and machine learning to make sense of big (aerial) data for disaster response. Big Data 4, 47–59 (2016).PubMed 

Google Scholar 
153.Simpson, R., Page, K. R. & De Roure, D. Zooniverse: observing the world’s largest citizen science platform. in Proceedings of the 23rd International Conference on World Wide Web 1049–1054 (2014).154.Pocock, M. J. O., Roy, H. E., Preston, C. D. & Roy, D. B. The biological records centre: a pioneer of citizen science. Biol. J. Linn. Soc. 115, 475–493 (2015).
Google Scholar  More

The improved performance of the mangrove LUE model considering coastal environments in this study was mainly attributed to the determination of environmental scalars. Parameters determining environmental stressors (e.g., Topt, Tmin, Tmax, VPDmin, and VPDmax) were set based on the general characteristics of mangroves worldwide. It may not be as accurate for the mangroves in our study sites, but it generally reflects the response of mangroves to environmental changes. Furthermore, as can be seen in Fig. S1, it is applicable to our study sites. Despite the specific characteristics of each mangrove ecosystem at different sites being preferred, this study first offers the possibility to estimate mangrove productivity at a larger scale to track GPP, thus emphasizing the role of mangrove ecosystems nationally or worldwide.The validation results showed that the LUE values of the mangrove model agreed well with those estimated by EC method (Fig. 3) and indicated improved performance (slope = 0.8218–1.0108, intercept = -0.0006–0.0052, R2 = 0.54–0.64, RMSE = 0.0051–0.0068, Pearson’s r = 0.73–1), compared to the MOD17 model (slope = 0.4993–0.5566, intercept = 0.0311–0.0313, R2 = 0.24–0.45, RMSE = 0.0217–0.0220, Pearson’s r = 0.45–0.49). Firstly, the RS-based LUE model for terrestrial ecosystems (MOD17) considers only the environmental stressors of Tair and VPD. The photosynthesis in mangrove forests is influenced by other unique environmental factors caused by tidal inundation. According to Fig. S3, PAR caused the most significant effect on LUE, which is consistent with previous studies14,30,32. The impact of SST has not been quantitatively assessed, however, SST is a critical control that determines the upper latitudinal range of mangrove ecosystems12,33. In our study, the effects of SST and salinity on the mangrove LUE were quantified and helped improve LUE modeling.Secondly, LUEmax was typically defined for different land covers, however, there were no specific values for mangrove forests. In this study, the LUEmax of mangroves was first determined. It is worth noting that daytime NEE responses to PAR vary depending on the Tair23,30,34 so that LUEmax was determined separately at high, optimal, and low temperatures. The results showed that LUEmax reached a maximum when Tair was within the optimal range for mangroves, which represents the high productivity of mangrove ecosystems. Furthermore, the estimated LUEmax of mangrove forests (0.057) was larger than most terrestrial forests35,36,37, which could contribute to the high production and carbon sequestration in mangrove forests.Lastly, the relatively low stomatal conductance of mangroves leads to low LSP compared with terrestrial forests, which could result in the high-irradiance stress for photosynthesis38,39. Mangrove LSP ranges from about 0.2–1.2 mmol/m2/s, depending on the species and environments40,41,42. LUE was relatively low in April and May when seasonal PAR was high, as photosynthesis is more likely to reach saturation. Therefore, we assumed the LUE of mangroves decreased with increasing PAR. In addition, we found that the downscaling effect of PAR on LUE was not constant, but varied with increasing PAR. As follows, different PAR scalars were set for mangroves according to different PAR values. This is a first attempt at refining PARscalar considering different solar radiation, which represents a significant departure from the assumption of a constant downscaling effect of PAR in RS-driven models14,43. The accuracy of the LUE model was improved by refining the PARscalar with different downscaling slopes, especially in periods of high PAR values.Compared with the results obtained from flux-tower measurements, the modeled GPP was basically within the confidence interval of the measured results. The annual averages of GPP in Zhangjiang were 1729 g C/m2/year and 1924 g C/m2/year, in 2012 and 2016, and the annual mean value of GPP in Zhanjiang was 1434 g C/m2/year in 2015. The previous study showed that the GPP in Zhangjiang ranged from 1763 to 1919 g C/m2/year with a mean value of 1871 g C/m2/year32,44,45, which is in good agreement with the estimated values obtained in this study. Liu and Lai46 reported that the GPP of the Mai Po mangrove reserve was 2827 g C/m2/year. Rodda, et al.20 found a GPP value of 1271 g C/m2/year for Sunderbans mangroves in India. Gnanamoorthy, et al.47 estimated a GPP of 2305 g C/m2/year for Pichavaram mangroves. Variations in these estimates across sites were possibly caused by different climate-hydrological conditions, mangrove species, and ages. Differences in the same location may be due to different time scales and different methods of data gap filling and flux partitioning.In a similar way to the GPP model for terrestrial ecosystems48, the effect of the mangrove GPP model on the accuracy of GPP estimates can vary considerably under different environmental conditions. However, in comparison with the accuracy of models built for other vegetation types, the GPP model in this study performed substantially in two sites with RMSE of 2.54–3.41 g C/m2/day. Wang et al.49 adopted different models to estimate GPP for global vegetation and validation results showed the RMSE ranged from 1.79 to 2.33 g C/m2/day. Xiao, et al.50 demonstrated that the deviation between observed and predicted GPP was about 35–282 g C/m2 in an evergreen needleleaf forest. Also, the absolute GPP errors were 7.94–20.92% and 9.97–13.70% for maize cropland and degraded grassland36. Despite the discrepancy, our results were generally consistent with previous studies and were verified by field observations near the flux towers.The comparison of MODIS GPP and EC-estimated GPP showed that the MODIS GPP had a large fluctuation and weakly reflected productivity, being overestimated in 2012 and underestimated in 2015. Different meteorological inputs, different environmental scalars and fraction of absorbed photosynthetic active radiation (fAPAR) products in MODIS GPP and our mangrove GPP model can explain their different results. However, the improvements in our GPP model may help to obtain more accurate GPP estimates. The response of mangrove productivity to Tair has not been well-calibrated in the MODIS GPP product, which may partly account for the poor correlation between the MODIS GPP and EC estimates. Besides, MODIS GPP product was developed based on the International Geosphere-Biosphere Programme (IGBP) land cover map, which doesn’t include mangroves as a specific land cover37. Therefore, LUEmax and environmental parameters were not defined for mangroves, which varied with different environments. This may lead to uncertainty in MODIS GPP product for mangrove forests14. However, the GPP model generated in our study showed similar trends to the field measurements, capturing seasonal variations. The increase in the difference between MODIS GPP and EC estimates may be due to the assumption that the increase in GPP is linear with respect to PAR. In our model, the response of GPP to PAR was suppressed, resulting in seasonal changes in GPP that better match the observations. In addition, the GPP derived from this study was in higher agreement with measured values compared with GPP estimated from the vegetation photosynthesis model (VPM), as shown in Fig. S4. The improvement of this model was more obvious in winter (December to February), which may be due to the environmental stress of SST and PAR. The VPM without considering SSTscalar and PARscalar overestimated GPP in winter. It is indicated that the performance of the mangrove GPP model in this study varied with season. It is recommended to improve the estimation of GPP in the future by considering the seasonal variation of mangrove forests when determining environmental variables.Most studies provide EC-based estimates of GPP that are measurements from a limited footprint. It is possible to extrapolate results across similar vegetation types and geographic settings, but not to areas of heterogeneous vegetation. The RS-based GPP model offers spatial-scale estimates that can be directly incorporated into ecosystem-type models. PAR, SST, and salinity are the key environmental parameters of this RS-based mangrove GPP model. SST and salinity data were derived from the satellite images, while PAR was generated from the reconstructed PAR data, since it is more accurate than the existing RS data and has historical year data. However, PAR products from Hamawari-8, MERIS, and SeaWiFS are available now, which provide an opportunity to obtain large-scale PAR data using RS in the future. In addition to this, GPP of two mangrove forests was assessed and validated with three-year measurements. Validation at different sites and years showed similar results, which indicated the model has similar performance across mangrove forests. Nonetheless, these estimates need to be corroborated with EC databases, which are relatively accurate and provide many additional variables that are currently beyond the scope of higher spatial-resolution RS estimates. The proposed GPP model considering coastal environments was well suited to extend the study area by incorporating RS information and meteorological data. Currently, there are still few mangrove carbon flux towers worldwide. The LUE and GPP models proposed in this study are difficult to validate with measurements from flux towers in other countries. However, local measurements are available in many countries with large mangrove forests, such as Thailand, Vietnam, India, and Bangladesh. Therefore, it is expected that comparisons with measurements from previous studies can be conducted to show the consistency and applicability.The LUE model considering the effects of SST, salinity, and PAR performed well, however, the GPP estimated from the LUE, fAPAR, and PAR showed discrepancies and were generally lower than the measured values. Although the results are better than MODIS GPP products, limitations exist still.Firstly, the effects of salinity and SST on mangrove productivity were directly related to tidal activities. The soil pore water and surface water salinity could affect the osmotic pressure of mangroves especially for the submerged parts which would control the stomatal conductance. In the same way, SST could influence the temperature of mangrove root systems and soil sediments which has impacts on mangrove roots’ respiration and transpiration. Although, theoretically, salinity and SST should be considered as environmental variables affecting mangrove LUE, our results (Fig. S3) indicated that salinity and SST have little influence on mangrove productivity51. To date, the quantitative impact of SST has not been comprehensively unfolded, but it is a global control that determines the upper limit of the latitudinal range of mangroves12,33. The weak relationships between salinity, SST, and mangrove GPP could be due to the uncertainty caused by tidal inundation. Tide duration, tide height, and tide cycle would determine the effect of salinity and SST on the mangrove LUE and GPP. However, quantifying the influence from the tidal cycle remains a challenging task, which could result in the relatively poor performance of Salinityscalar and SSTscalar as shown in Fig. S3. Quantifying the soil temperature and surface water salinity considering the tidal cycle will contribute to model the LUE and GPP of mangrove forests.Secondly, mangroves of different species and ages exhibit diverse structural and physical conditions, resulting in different LUEmax, and optimal growing conditions such as Topt and VPDmin. The environmental settings would also vary from region to region. Liu and Lai46 found that LUE increased slightly with the increasing salinity below 15 ppt (R2 = 0.16). However, it was noted that photosynthetic activity of mangroves would be inhibited when the surface water salinity was high30,51,52,53. Probably, the mutual relationship between LUE and salinity depends on the salinity level and mangrove species. However, we have not specified the variables for different mangrove species, ages and locations which could be improved in the future. Besides, there are multicollinearities between different environmental variables. For example, Tair may have effects on SST and VPD, but as shown in Fig. S5, they are all important for mangrove photosynthesis. However, the correlations between them are not clear and need to be quantified in the future.Thirdly, the relatively low spatial and temporal resolution of the environmental data from RS would influence the accuracy of the model. The datasets have a relatively coarse resolution (usually 500 m–1 km and daily) and are thereby less suitable for smaller nature reserves, especially in the narrow patches of mangrove areas that are rapidly being exploited in coastal China. Moreover, the variability in LUE decreases with increasing temporal scale54. In our study, we determined the PARscalar based on the response of LUE to hourly-scale PAR and found the different down-regulation effects with increasing PAR. However, this phenomenon is not obvious in previous studies. Most RS-based LUE models were developed at a daily or 8-day temporal scale6,50,55,56,57. In terrestrial forests, the light saturated effect caused by increasing PAR was neglectable with coarse temporal scale because the average PAR was usually lower than the LSP. However, as the time scale increases, the effect of light saturation on LUE becomes more pronounced32,58,59. More importantly, this effect is more obvious in mangroves due to their lower LSP18,38, which makes it important in mangrove LUE modeling. The results in Fig. 3 show similar performances of LUE model on hourly and daily scale. Thus, we suggested that our model can be adopted in hourly and daily temporal resolution. However, the PARscalar developed in this study was based on the mangrove forests in one study site which may be influenced by the mangrove species with different LSP and light conditions. What’s more, VPD was on a monthly scale, which cannot reflect environmental dynamics. However, the hourly and daily VPD data are currently not available for coastal areas in China. Therefore, we used monthly averages to represent daily VPD, which may lead to uncertainty in the derived GPP estimates (Figs. 6 and 7). Besides, porewater salinity is controlled by sea surface salinity, precipitation, and river discharge. However, currently, pore water salinity was expressed in terms of sea surface salinity, which may lead to an underestimation of Salinityscalar. More systematic study is necessary to make it more applicable and accurate on a large scale, of which modeling the LUE for different mangrove species and locations is inevitable. However, serving as a fundamental and preliminary step, our study aims to provide a framework for RS-based mangrove GPP modeling. Recently, with the advancement of satellite imagery, hourly-scale RS data for PAR, temperature and SST are available. It can be expected that our current work could be further improved by investigating the light saturation effects in different mangrove forests and adopt higher temporal resolution RS products such as Himawari-8 and GCOM-C in the future.Lastly, the overall underestimation of GPP was mainly caused by the underestimation of fAPAR. Even though the fAPAR computed from Sentinel-2 had higher resolution and accuracy than MODIS fAPAR products, future improvements are still needed. Sentinel-2 fAPAR products (fAPAR-S2) was calculated as the instantaneous fAPAR obtained at 10:00 local solar time which only roughly represented the daily average but was not accurate. Besides, RS-derived fAPAR only considers the absorptions by living green vegetation elements, whereas the ground measured fAPAR refers to the contributions from all absorbing components60. The lower fAPAR-S2 values in mangrove forests may be due to the exposed-to-air root systems which absorb the radiation. Moreover, the spatial distribution of PAR was determined by Co-Kriging interpolation. The elevation was taken as the covariate to estimate spatial PAR. There are many other variables affecting the incoming PAR (e.g., slope and clearness)61. A more comprehensive set of variables needs to be included in the Co-kriging interpolation to improve the PAR estimation.The spatial and seasonal variations of the mangrove GPP were related to environmental changes along the shoreline. The low summer GPP was explained by the lower fAPAR in summer compared with other seasons, which was principally due to the underestimation of fAPAR in summer. Furthermore, PARscalar took a mean value of LSP as 1 mmol/m2/s, however, LSP varied with different species and environmental conditions. In summer, mangroves are more likely to obtain light saturation, and thus PARscalar may lead to an underestimation of LUE and thus GPP. On the contrary, PAR values in winter were relatively low but increased slightly with decreasing latitude. Thus, the inhibitory effect of PAR on LUE was not significant, and GPP increased with decreasing latitude. Salinity and VPD were more stable across years and locations and had no noticeable effect on the mangrove LUE and GPP. The seasonal latitudinal patterns and effects on mangrove productivity were similar for Tair and SST. Tair and SST were lower in winter, especially at high latitudes where mangroves were more sensitive to cold weather. Therefore, the GPP of mangroves at high latitudes in winter was the lowest throughout the year. However, hot weather in summer also limited the photosynthesis in mangroves, especially at low latitudes, where Tair and SST were higher. Nevertheless, there were some correlations among these environmental constants. For example, the Tair affects the vapor pressure and SST. There was a positive correlation between PAR and Tair. The multicollinearity among these variables and the various conditions of mangroves may affect the performance of the model and show variations along the coastline, which would be improved in future studies.Additionally, the GPP of mangroves increased from 2007 to 2018, which was mainly due to the expansion of mangrove forests in the coastal areas. As mangroves grow, canopy size and tree density increase, which may lead to higher LUE and less underestimation of fAPAR, thus contributing to high productivity. However, Zhejiang province (27° 02′ N–31° 11′ N) experienced extremely cold weather in January 2016 caused by the East Asia cold wave62,63, and large areas of mangrove forests died or became sick, leading to a decline in the mangrove GPP at high latitudes in 2018. More