More stories

  • in

    Sex-based differences in the use of post-fire habitats by invasive cane toads (Rhinella marina)

    Study speciesCane toads (Rhinella marina) are large (to  > 1 kg) bufonids (Fig. 1a). Although native to north-eastern South America, these toads have been translocated to many countries worldwide to control insect pests12. Adult cane toads forage at night for insect prey and retreat to moist shelter-sites per day13. Small body size (and thus, high desiccation rate) restricts young toads to the margins of natal ponds14, but adult toads can survive even in highly arid habitats if they have access to water13,15. Cane toads prefer open habitats for foraging12, and thus can thrive in post-fire landscapes16,17. Cane toads in post-fire landscapes tend to have lower parasite burdens, probably because free-living larvae of their lungworm parasites cannot survive either the fire or the more sun-exposed post-fire landscape18.Figure 1taken from study sites between Casino, Grafton, and surrounds, NSW, by S.W. Kaiser.The cane toad Rhinella marina (a), and unburned, (b) and burned (c) habitats in which toads were collected and radio-tracked. Photographs were Full size imageStudy areaEast of the Great Dividing Range, near-coastal Clarence Dry Sclerophyll Forests of north-eastern New South Wales (NSW) are dominated by Spotted gum (Corymbia variegata) and Pink bloodwood (Corymbia intermedia)19. Fires are common, but typically cover relatively small areas before they are extinguished. In the summer of 2019–2020, however, prolonged drought followed by an unusually hot summer resulted in massive fires across this region, burning almost 100,000 km2 of vegetation9. In the current study, the toads we measured and dissected came from several sites within 75 km of the city of Casino (for site locations, see Fig. 2, Table 1, and18). The impacts of fire on faunal abundance and attributes shift with time since fire; for example, the abundance of a particular species may be reduced by fire (due to mortality from flames) but then increase as individuals from surrounding areas migrate to the recently-burned site to exploit new ecological opportunities provided by that landscape8. We chose to study this system 1-year post-fire, to allow time for such longer-term effects to be manifested.Figure 2Sampling sites relative to fire history. Sample sites are burned (red circles), and unburned (green squares). See Table 1 for key to sites. The legend shows the extent of burn a year prior to our study. Map created in QGIS 3.22.3. Fire history available from https://datasets.seed.nsw.gov.au/dataset/fire-extent-and-severity-mapping-fesm CC BY 4.0.Full size imageTable 1 Sampling sites and sample sizes for dissected and radio-tracked cane toads (Rhinella marina) in New South Wales, Australia.Full size tableSurveys of toad abundanceTo quantify toad abundance in burned and unburned sites, one observer (MJG) walked 100-m transects along roads at night (N = 23 and 8 respectively), recording all toads and native frogs (both adult and juvenile). The smaller number of unburned sites reflects the massive spatial scale of the wildfires, which made it difficult to find unburned areas. The transect sites were not the same as those sampled by “toad-busters” (below). We sampled both burned and unburned sites on each night, to de-confound effects of weather conditions with fire treatment. We scored frogs as well as toads to provide an estimate of overall anuran abundance and activity, and so that we could examine toad abundance relative to frog abundance as well as absolute toad numbers.“Toad-buster” sampleBecause of their ecological impact on native fauna, cane toads are culled by community groups as well as by government authorities12,20. We asked “toad-buster” groups to record whether the sites at which they collected toads had been burned during the 2019–2020 fires, or had remained unburned (Table 1). The toads were humanely euthanized (cooled-then-pithed: see21). The euthanasia method is brief (a few hours in the refrigerator, followed by pithing) and thus should not have affected any of the traits that we measured. For all of these toads, we measured body length (snout-urostyle length = SUL) and mass, and determined sex based on external morphology (skin colour and rugosity, nuptial pads: see22). A subset of toads (chosen to provide relatively equal numbers of males and females, and with equal numbers from burned and unburned sites) was dissected to provide data on mass of internal organs (fat bodies, liver, ovaries), reproductive condition (state of ovarian follicle development) and diet (mass and identity of prey items). To select the subsample of toads for dissection, we took relatively equal numbers of male and female toads from each bag of toads that was provided to us by the “toad-busters”. For logistical reasons, we were unable to dissect all of the toads that had been collected. Overall, we obtained data on morphology, diets and other traits from 481 fully dissected and 1443 partially dissected cane toads.Radio-trackingTo explore habitat use and movement patterns, we radio-tracked 57 toads over the course of two fieldtrips (0900–1800 h from 20 Nov 2021 to 6 Dec 2021 and 25 Jan 2022 to 10 Feb 2022). We selected seven sites (4 burned, 3 unburned) within 28 km of Tabbimoble, NSW (see Table 1 for locations and sample sizes of tracked toads). We hand-captured toads found active at night. These were measured, and their sex determined by external morphology (see above) and behaviour (release calls, given only by males: see23). We then fitted the toads with radio-transmitters (PD-2; Holohil Systems, Ontario, Canada; weighing ≤ 3.8 g) on cotton waist-belts, and released them at the site of capture. Tracked toads were 88.2–160.9 mm SUL (mass 70.1–546.3 g); thus, transmitters weighed  20 mm thick) within the quadrat, and estimated exposure of the toad within its refuge (the percentage of the animal’s body exposed to the naked eye). We then selected a compass bearing at random and walked 20 m in that direction where we rescored all of the above habitat attributes, to quantify habitat features in the broader environment (i.e., not just in microhabitats used by toads). We used those “random” sites to quantify overall habitat attributes of burned and unburned sites. Temperature was recorded by directing a temperature gun (Digitech QM7221) on (or otherwise close-to) toads and at a random point on the ground for random replicates. In total, we gathered radio-tracking data on movements and habitat variables from 57 cane toads, each of which was tracked for 5 days. Recaptured toads were euthanized by cooling-then-pithing.Morphological traitsTo obtain an index of body condition of toads, we regressed ln mass against ln SUL, and used the residual scores from that general linear regression as our estimate of body condition. Negative residual scores show an individual that weighs less-than-expected based on its body length. Likewise, we regressed mass of the fat bodies, liver and stomach against body mass to obtain indices of energy stores and stomach-content volumes relative to body mass. We scored male secondary sexual characteristics using the system of Bowcock et al.22. In their system, three sexually dimorphic traits (nuptial pad size, skin roughness and skin colouration) are scored from 0 to 2, and the scores from those three traits are summed to create a final value (on a 6-point scale) for the degree of elaboration of male-specific secondary sexual characteristics. We scored reproductive condition in adult female toads based on whether or not egg masses were visible during dissection, based on dissected toads from both “toad-buster” and telemetry samples.Statistical methodsData were analysed in R version 4.2.025. We used Linear Mixed Models (LMMs), Generalised Linear Mixed Models (GLMMs) and logistic regressions for our analyses. The R packages ‘tidyverse’26, ‘lmerTest’27, and ‘performance’28 were used.Habitat dataWe compared habitat variables between burned and unburned sites, and attributes of toads in burned versus unburned sites, using GLMMs (with negative binomial distribution) for count data (models were checked for overdispersion29) and LMMs on distance data, using ln-transformations where required to achieve normality. LMMs were used on non-normal percentage data, which were ln- and then logit-transformed (using log[(P + e)/(1 − P + e)], where e is the lowest non-zero number, halved)30. We used toad id, site (sampling location) and sampling trip (2019 versus 2020) as random factors.Anuran transect dataCounts of toads in burned versus unburned areas were compared both directly via GLMMs with a negative binomial distribution and relative to the numbers of frogs sighted along the same transects (binding the columns in R as ‘number of toads, number of amphibians – number of toads’ and using a GLMM with a binomial distribution). We used site as a random factor.Telemetry dataFor telemetry data, we analysed response variables via LMMs, and ln-transformed data where relevant to achieve normality.Dissection dataWe used LMMs for SUL, body mass, body condition and organ mass residuals (e.g., fat body mass relative to body mass). For prey item data, we used a poisson distribution with row number as a random factor, as the negative binomial and beta distribution GLMMs were overdispersed (see31). We used LMM for number of prey items and number of prey groups, with site as a random factor. Where models failed to converge, we reduced or removed the error term(s). Analyses were restricted to toads ≥ 70 mm SUL, because animals below this size were difficult to sex. We also performed nominal logistic regression to explore variation in sex ratio and male secondary sexual traits.Reproductive conditionWe used LMM for male secondary sexual characteristic display, using site as a random factor. For ovary presence, we used a binomial GLMM with a logit link, using site as a random factor. We used a LMM of the residual values from ovary mass relative to body mass (ln-transformed), using site as a random factor.Ethics declarationsAll procedures were performed in accordance with the relevant guidelines and regulations approved by Macquarie University Animal Ethics Committee (ARA Number: 2019/040-2) and in accordance with ARRIVE guidelines. More

  • in

    Global hydro-environmental lake characteristics at high spatial resolution

    Shiklomanov, I. A. & Rodda, J. C. World water resources at the beginning of the twenty-first century. (Cambridge University Press, 2003).Biggs, J., von Fumetti, S. & Kelly-Quinn, M. The importance of small waterbodies for biodiversity and ecosystem services: implications for policy makers. Hydrobiologia 793, 3–39 (2017).Article 

    Google Scholar 
    Heino, J. et al. Lakes in the era of global change: moving beyond single-lake thinking in maintaining biodiversity and ecosystem services. Biol. Rev. 96, 89–106 (2021).PubMed 
    Article 

    Google Scholar 
    Janssen, A. B. G. et al. Shifting states, shifting services: linking regime shifts to changes in ecosystem services of shallow lakes. Freshw. Biol. 66, 1–12 (2021).Article 

    Google Scholar 
    Knoll, L. B. et al. Consequences of lake and river ice loss on cultural ecosystem services. Limnol. Oceanogr. Lett. 4, 119–131 (2019).Article 

    Google Scholar 
    Sterner, R. W. et al. Ecosystem services of Earth’s largest freshwater lakes. Ecosyst. Serv. 41, 101046 (2020).Article 

    Google Scholar 
    Reynaud, A. & Lanzanova, D. A global meta-analysis of the value of ecosystem services provided by lakes. Ecol. Econ. 137, 184–194 (2017).PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    Cooley, S. W., Ryan, J. C. & Smith, L. C. Human alteration of global surface water storage variability. Nature 591, 78–81 (2021).ADS 
    CAS 
    PubMed 
    Article 

    Google Scholar 
    Downing, J. A. Global limnology: up-scaling aquatic services and processes to planet Earth. SIL Proceedings, 1922–2010 30, 1149–1166 (2009).Article 

    Google Scholar 
    Tranvik, L. J., Cole, J. J. & Prairie, Y. T. The study of carbon in inland waters—from isolated ecosystems to players in the global carbon cycle. Limnol. Oceanogr. Lett. 3, 41–48 (2018).Article 

    Google Scholar 
    Balsamo, G. et al. On the contribution of lakes in predicting near-surface temperature in a global weather forecasting model. Tellus A Dyn. Meteorol. Oceanogr. 64, 15829 (2012).Article 

    Google Scholar 
    DelSontro, T., Beaulieu, J. J. & Downing, J. A. Greenhouse gas emissions from lakes and impoundments: upscaling in the face of global change. Limnol. Oceanogr. Lett. 3, 64–75 (2018).CAS 
    Article 

    Google Scholar 
    Beaulieu, J. J. et al. Methane and carbon dioxide emissions from reservoirs: controls and upscaling. J. Geophys. Res. Biogeosciences 125, e2019JG005474 (2020).ADS 
    CAS 
    Article 

    Google Scholar 
    Slater, J. A. et al. The SRTM data “finishing” process and products. Photogramm. Eng. Remote Sens. 72, 237–247 (2006).Article 

    Google Scholar 
    Pekel, J.-F., Cottam, A., Gorelick, N. & Belward, A. S. High-resolution mapping of global surface water and its long-term changes. Nature 540, 418–422 (2016).ADS 
    CAS 
    PubMed 
    Article 

    Google Scholar 
    Verpoorter, C., Kutser, T., Seekell, D. A. & Tranvik, L. J. A global inventory of lakes based on high-resolution satellite imagery. Geophys. Res. Lett. 41, 6396–6402 (2014).ADS 
    Article 

    Google Scholar 
    Pickens, A. H. et al. Mapping and sampling to characterize global inland water dynamics from 1999 to 2018 with full Landsat time-series. Remote Sens. Environ. 243, 111792 (2020).ADS 
    Article 

    Google Scholar 
    Messager, M. L., Lehner, B., Grill, G., Nedeva, I. & Schmitt, O. Estimating the volume and age of water stored in global lakes using a geo-statistical approach. Nat. Commun. 7, 13603 (2016).ADS 
    CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    Tickner, D. et al. Bending the curve of global freshwater biodiversity loss: an emergency recovery plan. Bioscience 70, 330–342 (2020).PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    Downing, J. A., Polasky, S., Olmstead, S. M. & Newbold, S. C. Protecting local water quality has global benefits. Nat. Commun. 12, 1–6 (2021).Article 
    CAS 

    Google Scholar 
    Hill, R. A., Weber, M. H., Debbout, R. M., Leibowitz, S. G. & Olsen, A. R. The Lake-Catchment (LakeCat) Dataset: characterizing landscape features for lake basins within the conterminous USA. Freshw. Sci. 37, 208–221 (2018).PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    Soranno, P. A. et al. LAGOS-NE: a multi-scaled geospatial and temporal database of lake ecological context and water quality for thousands of US lakes. Gigascience 6, 1–22 (2017).PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    Toptunova, O., Choulga, M. & Kurzeneva, E. Status and progress in global lake database developments. Adv. Sci. Res. 16, 57–61 (2019).Article 

    Google Scholar 
    Meyer, M. F., Labou, S. G., Cramer, A. N., Brousil, M. R. & Luff, B. T. The global lake area, climate, and population dataset. Sci. Data 7, 174 (2020).PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    Kling, G. W., Kipphut, G. W., Miller, M. M. & O’Brien, W. J. Integration of lakes and streams in a landscape perspective: the importance of material processing on spatial patterns and temporal coherence. Freshw. Biol. 43, 477–497 (2000).Article 

    Google Scholar 
    Fergus, C. E. et al. The freshwater landscape: lake, wetland, and stream abundance and connectivity at macroscales. Ecosphere 8, e01911 (2017).Article 

    Google Scholar 
    Lehner, B., Messager, ML., Korver, MC. & Linke, S. LakeATLAS Version 1.0, figshare, https://doi.org/10.6084/m9.figshare.19312001 (2022).Linke, S. et al. Global hydro-environmental sub-basin and river reach characteristics at high spatial resolution. Sci. data 6, 283 (2019).PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    Fergus, C. E. et al. National framework for ranking lakes by potential for anthropogenic hydro-alteration. Ecol. Indic. 122, 107241 (2021).Article 

    Google Scholar 
    Bracht-Flyr, B., Istanbulluoglu, E. & Fritz, S. A hydro-climatological lake classification model and its evaluation using global data. J. Hydrol. 486, 376–383 (2013).ADS 
    Article 

    Google Scholar 
    Soranno, P. A. et al. Using landscape limnology to classify freshwater ecosystems for multi-ecosystem management and conservation. Bioscience 60, 440–454 (2010).Article 

    Google Scholar 
    McCullough, I. M., Skaff, N. K., Soranno, P. A. & Cheruvelil, K. S. No lake left behind: how well do U.S. protected areas meet lake conservation targets? Limnol. Oceanogr. Lett. 4, 183–192 (2019).Article 

    Google Scholar 
    Stanley, E. H. et al. Biases in lake water quality sampling and implications for macroscale research. Limnol. Oceanogr. 64, 1572–1585 (2019).ADS 
    CAS 
    Article 

    Google Scholar 
    Hanson, P. C., Weathers, K. C. & Kratz, T. K. Networked lake science: how the Global Lake Ecological Observatory Network (GLEON) works to understand, predict, and communicate lake ecosystem response to global change. Inl. Waters 6, 543–554 (2016).Article 

    Google Scholar 
    Lottig, N. R. & Carpenter, S. R. Interpolating and forecasting lake characteristics using long-term monitoring data. Limnol. Oceanogr. 57, 1113–1125 (2012).ADS 
    Article 

    Google Scholar 
    Filazzola, A. et al. A database of chlorophyll and water chemistry in freshwater lakes. Sci. Data 2020 71 7, 1–10 (2020).
    Google Scholar 
    Lehner, B. & Messager, M. L. HydroLAKES – Technical Documentation Version 1.0. https://data.hydrosheds.org/file/technical-documentation/HydroLAKES_TechDoc_v10.pdf (2016).Natural Resources Canada. CanVec Hydrography: Waterbody Features. Version 12.0. https://ftp.maps.canada.ca/pub/nrcan_rncan/vector/canvec (2013).Lehner, B., Verdin, K. & Jarvis, A. New global hydrography derived from spaceborne elevation data. Eos, Trans. AGU 89, 93–94 (2008).ADS 
    Article 

    Google Scholar 
    Farr, T. G. & Kobrick, M. Shuttle radar topography mission produces a wealth of data. Eos, Trans. AGU 81, 583–585 (2000).ADS 
    Article 

    Google Scholar 
    Müller Schmied, H. et al. The global water resources and use model WaterGAP v2.2d: model description and evaluation. Geosci. Model Dev. 14, 1037–1079 (2021).ADS 
    Article 

    Google Scholar 
    Beck, H. E. et al. Global evaluation of runoff from 10 state-of-the-art hydrological models. Hydrol. Earth Syst. Sci. 21, 2881–2903 (2017).ADS 
    Article 

    Google Scholar 
    Alcamo, J. et al. Development and testing of the WaterGAP 2 global model of water use and availability. Hydrol. Sci. J. 48, 317–338 (2003).Article 

    Google Scholar 
    Döll, P., Kaspar, F. & Lehner, B. A global hydrological model for deriving water availability indicators: model tuning and validation. J. Hydrol. 270, 105–134 (2003).ADS 
    Article 

    Google Scholar 
    Lehner, B. & Grill, G. Global river hydrography and network routing: baseline data and new approaches to study the world’s large river systems. Hydrol. Process. 27, 2171–2186 (2013).ADS 
    Article 

    Google Scholar 
    Fick, S. E. & Hijmans, R. J. WorldClim 2: new 1-km spatial resolution climate surfaces for global land areas. Int. J. Climatol. 37, 4302–4315 (2017).Article 

    Google Scholar 
    Hengl, T. et al. SoilGrids250m: Global gridded soil information based on machine learning. PLoS One 12, e0169748 (2017).PubMed 
    PubMed Central 
    Article 
    CAS 

    Google Scholar 
    Zhang, X. et al. GLC_FCS30: Global land-cover product with fine classification system at 30 m using time-series Landsat imagery. Earth Syst. Sci. Data 13, 2753–2776 (2021).ADS 
    Article 

    Google Scholar 
    Buchhorn, M. et al. Copernicus Global Land Service: Land Cover 100m: Collection 3: epoch 2019: Globe, Zenodo, https://doi.org/10.5281/zenodo.3939050 (2020).ESRI. ArcGIS Desktop: Release 10.4.1 (Environmental Systems Research Institute, Redlands, CA, USA, 2016).Soranno, P. A., Cheruvelil, K. S., Wagner, T., Webster, K. E. & Bremigan, M. T. Effects of land use on lake nutrients: the importance of scale, hydrologic connectivity, and region. PLoS One 10, e0135454 (2015).PubMed 
    PubMed Central 
    Article 
    CAS 

    Google Scholar 
    Su, Z. H., Lin, C., Ma, R. H., Luo, J. H. & Liang, Q. O. Effect of land use change on lake water quality in different buffer zones. Appl. Ecol. Environ. Res. 13, 639–653 (2015).
    Google Scholar 
    Brakebill, J. W., Schwarz, G. E. & Wieczorek, M. E. An enhanced hydrologic stream network based on the NHDPlus medium resolution dataset. Scientific Investigations Report https://doi.org/10.3133/sir20195127 (2020).Carroll, M., Townshend, J., DiMiceli, C., Noojipady, P. & Sohlberg, R. Global raster water mask at 250 meter spatial resolution, Collection 5: MOD44W MODIS Water Mask. College Park, Maryland: University of Maryland (2009).Carroll, M. L., Townshend, J. R., DiMiceli, C. M., Noojipady, P. & Sohlberg, R. A. A new global raster water mask at 250 m resolution. Int. J. Digit. Earth 2, 291–308 (2009).ADS 
    Article 

    Google Scholar 
    European Environment Agency (EEA). European Catchments and Rivers Network System (ECRINS), https://www.eea.europa.eu/data-and-maps/data/european-catchments-and-rivers-network (2012).Ouellet Dallaire, C., Lehner, B., Sayre, R. & Thieme, M. A multidisciplinary framework to derive global river reach classifications at high spatial resolution. Environ. Res. Lett. 14, 024003 (2019).ADS 
    Article 

    Google Scholar 
    Global Runoff Data Centre (GRDC). River discharge data. Federal Institute of Hydrology, 56068 Koblenz, Germany, https://www.bafg.de/GRDC (2014).Openshaw, S. The modifiable areal unit problem. In Quantitative Geography: A British View (eds. Wrigley, N. & Bennett, R.) 60–69 (Routledge and Kegan Paul, Andover, 1981).United States Census Bureau. 2010 Census. ftp://ftp2.census.gov/geo/tiger (2010).Center for International Earth Science Information Network (CIESIN) & NASA Socioeconomic Data and Applications Center (SEDAC). Gridded Population of the World, Version 4 (GPWv4): Population Count and Density. https://doi.org/10.7927/H4JW8BX5 (2016).Grill, G. et al. Mapping the world’s free-flowing rivers. Nature 569, 215–221 (2019).ADS 
    CAS 
    PubMed 
    Article 

    Google Scholar 
    Allen, D. J. et al. The Diversity of Life in African Freshwaters: Under Water, Under Threat: an Analysis of the Status and Distribution of Freshwater Species Throughout Mainland Africa. (IUCN, 2011).Markovic, D. et al. Europe’s freshwater biodiversity under climate change: distribution shifts and conservation needs. Divers. Distrib. 20, 1097–1107 (2014).Article 

    Google Scholar 
    Fluet-Chouinard, E., Lehner, B., Rebelo, L.-M., Papa, F. & Hamilton, S. K. Development of a global inundation map at high spatial resolution from topographic downscaling of coarse-scale remote sensing data. Remote Sens. Environ. 158, 348–361 (2015).ADS 
    Article 

    Google Scholar 
    Lehner, B. et al. High‐resolution mapping of the world’s reservoirs and dams for sustainable river‐flow management. Front. Ecol. Environ. 9, 494–502 (2011).Article 

    Google Scholar 
    Fan, Y., Li, H. & Miguez-Macho, G. Global patterns of groundwater table depth. Science 339, 940–943 (2013).ADS 
    CAS 
    PubMed 
    Article 

    Google Scholar 
    Robinson, N., Regetz, J. & Guralnick, R. P. EarthEnv-DEM90: A nearly-global, void-free, multi-scale smoothed, 90m digital elevation model from fused ASTER and SRTM data. ISPRS J. Photogramm. Remote Sens. 87, 57–67 (2014).ADS 
    Article 

    Google Scholar 
    Metzger, M. J. et al. A high-resolution bioclimate map of the world: a unifying framework for global biodiversity research and monitoring. Glob. Ecol. Biogeogr. 22, 630–638 (2013).Article 

    Google Scholar 
    Hijmans, R. J., Cameron, S. E., Parra, J. L., Jones, P. G. & Jarvis, A. Very high resolution interpolated climate surfaces for global land areas. Int. J. Climatol. 25, 1965–1978 (2005).Article 

    Google Scholar 
    Zomer, R. J., Trabucco, A., Bossio, D. A. & Verchot, L. V. Climate change mitigation: a spatial analysis of global land suitability for clean development mechanism afforestation and reforestation. Agric. Ecosyst. Environ. 126, 67–80 (2008).Article 

    Google Scholar 
    Trabucco, A., Zomer, R. J., Bossio, D. A., van Straaten, O. & Verchot, L. V. Climate change mitigation through afforestation/reforestation: a global analysis of hydrologic impacts with four case studies. Agric. Ecosyst. Environ. 126, 81–97 (2008).Article 

    Google Scholar 
    Trabucco, A. & Zomer, R. J. Global soil water balance geospatial database. CGIAR Consortium for Spatial Information, https://cgiarcsi.community/data/global-high-resolution-soil-water-balance (2010).Hall, D. K., Riggs, G. A. & Salomonson, V. MODIS/Terra snow cover daily L3 global 500m grid, version 5, 2002–2015, https://doi.org/10.5067/MODIS/MOD10A1.006 (2016).Bartholomé, E. & Belward, A. S. GLC2000: a new approach to global land cover mapping from Earth observation data. Int. J. Remote Sens. 26, 1959–1977 (2005).Article 

    Google Scholar 
    Ramankutty, N. & Foley, J. A. Estimating historical changes in global land cover: Croplands from 1700 to 1992. Global Biogeochem. Cycles 13, 997–1027 (1999).ADS 
    CAS 
    Article 

    Google Scholar 
    Lehner, B. & Döll, P. Development and validation of a global database of lakes, reservoirs and wetlands. J. Hydrol. 296, 1–22 (2004).ADS 
    Article 

    Google Scholar 
    Ramankutty, N., Evan, A. T., Monfreda, C. & Foley, J. A. Farming the planet: 1. Geographic distribution of global agricultural lands in the year 2000. Global Biogeochem. Cycles 22, (2008).Siebert, S. et al. A global data set of the extent of irrigated land from 1900 to 2005. Hydrol. Earth Syst. Sci. 19, 1521–1545 (2015).ADS 
    Article 

    Google Scholar 
    GLIMS & NSIDC. Global land ice measurements from space (GLIMS) glacier database, v1. National Snow and Ice Data Center (NSIDC), https://doi.org/10.7265/N5V98602 (2012).Gruber, S. Derivation and analysis of a high-resolution estimate of global permafrost zonation. Cryosphere 6, 221–233 (2012).ADS 
    Article 

    Google Scholar 
    UNEP-WCMC & IUCN. The World Database on Protected Areas, http://www.protectedplanet.net (2014).Dinerstein, E. et al. An ecoregion-based approach to protecting half the terrestrial realm. Bioscience 67, 534–545 (2017).PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    Abell, R. et al. Freshwater ecoregions of the world: a new map of biogeographic units for freshwater biodiversity conservation. Bioscience 58, 403–414 (2008).Article 

    Google Scholar 
    Hengl, T. et al. SoilGrids1km—global soil information based on automated mapping. PLoS One 9, e105992 (2014).ADS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    Hartmann, J. & Moosdorf, N. The new global lithological map database GLiM: a representation of rock properties at the Earth surface. Geochem. Geophys. Geosyst. 13, Q12004 (2012).ADS 
    Article 

    Google Scholar 
    Williams, P. W. & Ford, D. C. Global distribution of carbonate rocks. Zeitschrift für Geomorphologie Suppl. 147, 1–2 (2006).
    Google Scholar 
    Borrelli, P. et al. An assessment of the global impact of 21st century land use change on soil erosion. Nat. Commun. 8, 1–13 (2017).ADS 
    CAS 
    Article 

    Google Scholar 
    Pesaresi, M. & Freire, S. GHS Settlement grid following the REGIO model 2014 in application to GHSL Landsat and CIESIN GPW v4-multitemporal (1975-1990-2000-2015). European Commission, Joint Research Centre (JRC), https://data.europa.eu/data/datasets/jrc-ghsl-ghs_smod_pop_globe_r2016a (2016).Doll, C. N. H. CIESIN thematic guide to night-time light remote sensing and its applications. CIESIN http://sedac.ciesin.columbia.edu/binaries/web/sedac/thematic-guides/ciesin_nl_tg.pdf (2008).Meijer, J. R., Huijbregts, M. A. J., Schotten, K. C. G. J. & Schipper, A. M. Global patterns of current and future road infrastructure. Environ. Res. Lett. 13, 64006 (2018).Article 

    Google Scholar 
    Venter, O. et al. Global terrestrial Human Footprint maps for 1993 and 2009. Sci. data 3, 160067 (2016).PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    University of Berkeley. Database of global administrative areas (GADM). University of Berkeley, Museum of Vertebrate Zoology and the International Rice Research Institute, http://www.gadm.org (2012).Kummu, M., Taka, M. & Guillaume, J. H. A. Gridded global datasets for gross domestic product and Human Development Index over 1990–2015. Sci. data 5, 180004 (2018).PubMed 
    PubMed Central 
    Article 

    Google Scholar  More

  • in

    Bateman gradients from first principles

    Model 1: Evolution of multiple mating and mate monopolisation under ancestral monogamyIn all models, I assume a large population with a 1:1 sex ratio. I begin with what is possibly the simplest model set-up for deriving Bateman functions in a scenario that is completely symmetrical aside from gamete number. Assume a monogamous, externally fertilising population where parents pair up and release their gametes into a nest. That is, every individual in the initial population participates in exactly one fertilisation event (the equivalent of a mating). Now consider a mutant individual that can attract multiple mates of the opposite type to release gametes into its nest, with no competition from other individuals of its own type. This simple set-up avoids asymmetries arising from internal fertilisation, and the complication of direct gamete competition for the multiply mating mutant individual (which is examined in Models 2–3), placing focus directly on the core of the problem: the asymmetry arising in fertilisation from imbalanced gamete numbers. All gametes are released in one burst by all individuals, but the focal individual may achieve ‘multiple matings’ simply by monopolising multiple mates at its nest. The reproductive success of the focal individual is then equivalent to the number of fertilisations that take place in that nest. Our aim is to understand how the reproductive success of an individual deviating from the monogamous population strategy and instead mating with (hat{m}) individuals of the opposite type is altered. A strong positive relationship between (hat{m}) and reproductive success then indicates a steep Bateman gradient. If Bateman’s assertion is correct, the resulting gradient should be steeper for the type that produces the larger number of gametes. Note that there is a game-theoretical25 flavour to this setting, where the focus is on the fitness of a rare mutant in a population with a fixed resident strategy.The two types are labelled with x and y, which could correspond to the two sexes, depending on what gamete numbers are assigned to them. The number of gametes produced by a single individual is labelled nx and ny, and the total number of gametes in a nest (or more generally, a fertilisation arena which could be internal or external) is labelled with Nx and Ny. To compute the number of fertilisations in a nest with a total of Nx and Ny gametes, I use a fertilisation function first derived by Togashi et al.24 purely from biophysical principles, treating the two gamete types symmetrically, with no pre-existing assumptions about differences between females and males or their gametes (for a broader context and comparison to other functions, see Table 1 and function F7 in19). Any sex-specific differences arise only retrospectively after different gamete numbers are assigned to x and y of which either one could be male or female. The fertilisation function is (fleft({N}_{x},{N}_{y}right)={N}_{x}{N}_{y}frac{{e}^{a{N}_{x}}-{e}^{a{N}_{y}}}{{{N}_{x}e}^{a{N}_{x}}-{N}_{y}{e}^{a{N}_{y}}}), where a is a parameter controlling fertilisation efficiency (for the special case Nx = Ny the function is defined as (fleft({N}_{x},{N}_{y}right)=frac{a{N}_{x}^{2}}{1+a{N}_{x}})19,24, which is also the limit of f when Ny → Nx).In a monogamous resident pair, we have simply Nx = nx and Ny = ny. But if a mutant individual of type x is able to attract (hat{m}) fertilisation partners of type y, then for that individual ({N}_{y}=hat{m}{n}_{y}), and the corresponding Bateman function is$${b}_{x}left(hat{m}right)=fleft({N}_{x},{N}_{y}right)=fleft({n}_{x},hat{m}{n}_{y}right)$$
    (1)
    where the fertilisation function f is as described above. Because of symmetry, the corresponding function for y is found simply by swapping x and y. This function can reproduce the characteristic Bateman gradient asymmetry as gamete numbers diverge (progressing from isogamy to anisogamy in Fig. 1), showing how Bateman’s assertion follows from biophysical effects that arise from unequal numbers of fusing particles: the fertilisation function f is derived solely from such biophysical effects, not from any sex-specific assumptions. Equation (1) makes no reference to sexes, and they only become specified when values are assigned to nx and ny. For example, if nx = 10 and ny = 10,000, the female Bateman function is ({b}_{x}left(hat{m}right)) and the male Bateman function ({b}_{y}left(hat{m}right)), where for the latter all xs in Eq. (1) are replaced with ys and vice versa. The labels x and y are truly just labels. While there are inevitably assumptions built into the equations, crucially we can be certain there are no sex-specific assumptions. Yet the typical shapes reminiscent of Bateman gradients arise from the model when different values are specified for nx and ny (Fig. 1).Fig. 1: The Bateman function of Eq. (1).This figure shows how the basic Bateman gradient asymmetry arises from simple biophysics and mathematics of fertilisation. The population is monogamous aside from a mutant individual, whose number of fertilisation partners (‘matings’) varies on the horizontal axes within panels. a–d show the effect of variation in sex-specific gamete numbers under efficient fertilisation, while e–h show the effect of variation in sex-specific gamete numbers under inefficient fertilisation. Parameter values used are shown in the figure. Females (gamete number nx) are indicated by blue crosses and connecting lines, while males (gamete number ny) are indicated by black dots and connecting lines. Under isogamy, females and males are undefined, and the two colours overlap. The typical sex-specific shapes of Bateman gradients arise from a single equation (which itself is not sex-specific) when a difference in gamete numbers is assigned to nx and ny, confirming Bateman’s intuition that the primary cause of the difference in selection is that females produce fewer gametes than males. The entire range of gamete number ratios presented in the figure is observed in nature, from equal gamete size in many unicellular organisms39 to vertebrates, where sperm count per ejaculate can commonly exceed 109 (see ref. 40 and Supplementary Information therein).Full size imageGamete limitation changes the results quantitatively so that under conditions of poor fertilisation efficiency a larger imbalance in gamete numbers is needed for Bateman gradients to diverge to a similar extent. However, even under inefficient fertilisation, the Bateman gradients do not reverse.Model 2: An external fertiliser model with population-level polygamy and gamete competitionModel 1 presented the simplest possible scenario, where all individuals except a rare mutant mate only once, and gamete competition (sperm competition26, but without assigning either gamete type to be sperm) was thus excluded for the focal mutant individual. Now I generalise from this to a situation that remains entirely symmetrical, but where the resident number of matings can take on any value, and then derive the Bateman function for a rare mutant that deviates from this population-level value. This set-up allows for gamete competition for the focal mutant individual, a crucial addition because of the empirical and theoretical importance of sperm competition26, as well as earlier theory suggesting that polyandry decreases the sex difference in Bateman gradients2.The biological set-up is such that there is a large population and a large number of patches (fertilisation arenas) where multiple individuals of both sexes can release their gametes for fertilisation. After all individuals have released their gametes, those in each patch mix freely and fertilisations take place randomly. Set up in this way, the model is again identical from the perspective of both sexes, and gamete number can be isolated as the sole possible causal factor in any subsequent differences that may arise, extending from the initially monogamous and gamete competition-free scenario of Model 1. All individuals of both sexes are assumed to initially have the same strategy: to divide their nx or ny gametes equally between m patches, and distribute themselves in such a way that gametes from m individuals of each type release gametes into each patch (the number of individuals of each sex per patch need not necessarily be strictly equal to m, but this is the simplest assumption to account for the fact that gamete competition tends to increase with multiple ‘matings’). Now, if a rare x mutant divides its gametes evenly into (hat{m}) randomly selected patches, its gamete number per patch and consequently competitiveness in each patch is altered. Therefore, gametes of a mutant of type x will gain, on average, a fraction ({c}_{x}=left({n}_{x}/hat{m}right)/{N}_{x}) of the fertilisations in that patch, where ({N}_{x}={n}_{x}/hat{m}+(m-1){n}_{x}/m). To compute the number of realised fertilisations in a patch, I use the same fertilisation function as in Model 1, where the mutant number of gametes in a patch is Nx as above and the number of gametes of the opposite type is ({N}_{y}=mfrac{{n}_{y}}{m}={n}_{y}). All the components are now in place to write down the Bateman function corresponding to this scenario, for a mutant of type x:$${b}_{x}left(hat{m},mright)=hat{m}{c}_{x}fleft({N}_{x},{N}_{y}right)$$
    (2)
    where cx, Nx and Ny are as defined above, and the fertilisation function f is as in Model 1. For completeness, define bx(0, m) = 0, which is necessarily true, but useful to define separately because division by 0 renders Eq. (2) formally undefined when (hat{m}=0).As in Model 1, Eq. (2) makes no reference to sexes, and they only become specified when values are assigned to nx and ny (Fig. 2).Fig. 2: The Bateman function of Eq. (2) for an externally fertilising population with potential for population-wide polygamy and gamete competition.Results are shown for two values of resident matings (m = 1 and m = 2). a–h show the effect of variation in sex-specific gamete numbers and in fertilisation efficiency with m = 1, while i–p show the same with m = 2. Parameter values used are shown in the figure. The value m = 2 is used here because it is comparable to the mean number of matings in Bateman’s1 work (see Fig. 3 for corresponding results with internal fertilisation, but note that the aim of the models is not to quantitatively reproduce Bateman’s results). Females (gamete number nx) are indicated by blue crosses and connecting lines, while males (gamete number ny) are indicated by black dots and connecting lines. Under isogamy, females and males are undefined, and the two colours overlap. Further variation in m is examined in Fig. 4.Full size imageModel 3: An internal fertiliser modelModels 1–2 were set up with the central aim of full symmetry and exclusion of any sex-specific assumptions. Internal fertilisation breaks this symmetry by introducing a sex-specific assumption other than gamete number. Bateman gradients are, however, most commonly applied to situations with internal fertilisation where females are gamete recipients and males are gamete donors27. I therefore construct a model accounting for internal fertilisation. Where Eqs. (1) and (2) allowed no sex differences aside from gamete number, here I additionally consider the fact that females receive gametes while males donate them.As in model 2, there is a very large population, and I assume that in the resident population, all females and males mate exactly m times. It is then considered how a rare mutant individual’s (of either sex) fitness depends on its number of matings (hat{m}).I use the same fertilisation function as in Models 1-2. Consider first the female perspective (labelled with x). A female produces nx gametes and retains them internally. Each female mates with m males, who also mate with m females, dividing their gametes evenly over these matings. Therefore a mutant female receives (hat{m}frac{{n}_{y}}{m}) male gametes, and her reproductive success is$${b}_{x}left(hat{m},mright)=fleft({n}_{x},hat{m}frac{{n}_{y}}{m}right)$$
    (3)
    A mutant male, on the other hand, mates with (hat{m}) females, each of which mates with m−1 additional males. Therefore, the mutant male’s mating partners will receive a total of ({{N}_{y}=n}_{y}/hat{m}+(m-1){n}_{y}/{m}) male gametes. Thus, the mutant male gains a fraction ({c}_{y}=left({n}_{y}/hat{m}right)/{N}_{y}) of the fertilisations with each female, while the total reproductive success per female is f(nx,Ny). The mutant male’s reproductive success is therefore$${b}_{y}left(hat{m},mright)=hat{m}{c}_{y}fleft({n}_{x},{N}_{y}right)$$
    (4)
    To avoid division by 0, we can again define by (0, m) = 0, analogous to Model 2. In contrast to Models 1–2, there are now separate equations for each sex because of the additional sex-specific assumption of internal fertilisation, but no further sex-specific assumptions are used in their derivation. Visually the Bateman functions (Fig. 3) are nevertheless very similar to Model 2, and again reproduce the sex-specific shapes first proposed by Bateman1 when fertilisation is efficient. However, an interesting exception arises when relatively weak asymmetry in gamete numbers is combined with inefficient fertilisation and gamete limitation. When these conditions are combined with internal fertilisation, Bateman gradients can theoretically be reversed.Fig. 3: The Bateman functions of Eqs. (3) and (4) for internal fertilisation.Where Figs. 1 and 2 show that the sex-specific shapes of Bateman functions are ultimately caused by differences in gamete number, Fig. 3 shows that internal fertilisation does not invalidate this outcome when fertilisation is efficient. As in Fig. 2, results are shown for two values of resident matings (1 and 2), and the value m = 2 is used because it is comparable to the mean number of matings in Bateman’s1 work. a–h show the effect of variation in sex-specific gamete numbers and in fertilisation efficiency with m = 1, while i–p show the same with m = 2. Parameter values used are shown in the figure. Inefficient fertilisation combined with relatively low asymmetry in gamete numbers and the added asymmetry of internal fertilisation can in principle reverse the Bateman gradients (second and fourth row). Females (gamete number nx) are indicated by blue crosses and connecting lines, while males (gamete number ny) are indicated by black dots and connecting lines.Full size image More

  • in

    Ecological networks of dissolved organic matter and microorganisms under global change

    Experimental designThe comparative field microcosm experiments were conducted on Laojun Mountain in China (26.6959 N; 99.7759 E) in September–October 2013, and on Balggesvarri Mountain in Norway (69.3809 N; 20.3483 E) in July 2013, designed to be broadly representative of subtropical and subarctic climatic zones, respectively, as first reported in Wang et al.29. In the Laojun Mountain region, mean annual temperatures ranged from 4.2 to 12.9 °C, with July mean temperatures of 17–25 °C. In the Balggesvarri Mountain region, mean annual temperatures ranged from −2.9 to 0.7 °C, with July mean temperatures of 8–16 °C. The experiments were characterised by an aquatic ecosystem with consistent initial DOM composition but different locally colonised microbial communities and newly produced endogenous DOM. While allowing us to minimise the complexity of natural ecosystems, the experiment provided a means for investigating DOM-microbe associations at large spatial scales by controlling the initial DOM supply. Briefly, we selected locations with five different elevations on each mountainside. The elevations were 3822, 3505, 2915, 2580 and 2286 m a.s.l. on Laojun Mountain in China, and 750, 550, 350, 170 and 20 m a.s.l. on Balggesvarri Mountain in Norway. At each elevation, we established 30 aquatic microcosms (1.5 L bottle) composed of 15 g of sterilised lake sediment and 1.2 L of sterilised artificial lake water at one of ten nutrient levels of 0, 0.45, 1.80, 4.05, 7.65, 11.25, 15.75, 21.60, 28.80 and 36.00 mg N L−1 of KNO3 in the overlying water. To compensate for nitrate additions shifting stoichiometric ratios, KH2PO4 was added to the bottles so that the N/P ratio of the initial overlying water was 14.93, which was similar to the annual average ratio in Taihu Lake during 2007 (that is, 14.49). Thus, we use “nutrient enrichment” to indicate a series of targeted nutrient levels of both nitrate and phosphate, the former of which was used to represent nutrient enrichment in the statistical analyses. Each nutrient level was replicated three times. The lake sediments were obtained from the centre of Taihu Lake, China, and were aseptically canned per bottle after autoclaving at 121 °C for 30 min. Nutrient levels for the experiments were selected based on conditions of the eutrophic Taihu Lake, and the highest nitrate concentration was based on the maximum total nitrogen in 2007 (20.79 mg L−1; Fig. S19). We chose the nutrient level of this year because a massive cyanobacteria bloom in Taihu Lake happened in May 2007 and initiated an odorous drinking water crisis in the nearby city of Wuxi.The microcosms were left in the field for one month allowing airborne bacteria to freely colonise the sediments and water. To keep the microbial dispersal events as natural as possible, we did not cover the experimental microcosms in case of rainfall. To avoid or minimize potential influence of extreme nature events, we (i) left the top 20% of each microcosm empty to prevent water from overflowing during heavy rains, and (ii) checked the experimental sites twice during each experimental period, and added sterilized water to obtain a final volume of approximately 1.2 L. The bottom of our microcosm was buried into the local soils by 10% of the bottle height, partly to reduce UV exposure to sediments. More considerations of the experimental design were detailed in the Supplementary Methods. To avoid the effects of daily temperature variation, we measured the water temperature and pH within 2 h before noon at all elevations in the day before the final sample collection. At the end of the experimental period, we aseptically sampled the water and sediments of the 300 bottles (that is, 2 mountains × 5 elevations × 10 nutrient levels × 3 replicates) for the following analyses of physiochemical variables, bacterial community and DOM composition.Physiochemical variables and bacterial communityWe measured environmental variables, namely, the total nitrogen (TN), total phosphorus (TP), dissolved nutrients (that is, NOx−, NO2−, NH4+ and PO43−), total organic carbon (TOC), dissolved organic carbon (DOC) and chlorophyll a (Chl a) in the sediments, and the NO3−, NO2−, NH4+, PO43− and pH in the overlying water (Table S2, Fig. S20), according to Wang et al.29.The sediment bacteria were examined using high-throughput sequencing of 16S rRNA genes. The sequences were processed in QIIME (v1.9)45 and OTUs were defined at 97% sequence similarity. The bacterial sequences were rarefied to 20,000 per sample. Further details on physicochemical and bacterial community analyses are available in Wang et al.29.ESI FT-ICR MS analysis of DOM samplesHighly accurate mass measurements of DOM within the sediment samples were conducted using a 15 Tesla solariX XR system, a ultrahigh-resolution Fourier transform ion cyclotron resonance mass spectrometer (FT-ICR MS, Bruker Daltonics, Billerica, MA) coupled with an electrospray ionization (ESI) interface, as demonstrated previously46 with some modifications. It should be noted that FT-ICR MS does not identify molecules, but only molecular formulae in terms of elemental composition and there can be many molecular structures sharing the same elemental compositions. DOM was solid-phase extracted (SPE) with Agilent VacElut resins before FT-ICR MS measurement47 with minor modifications. Briefly, an aliquot of 0.7 g freeze-dried sediment was sonicated with 30 ml ultrapure water for 2 h, and centrifuged at 5000 × g for 20 min. The extracted water was filtered through the 0.45 μm Millipore filter and further acidified to pH 2 using 1 M HCl. Cartridges were drained, rinsed with ultrapure water and methanol (ULC-MS grade), and conditioned with pH 2 ultrapure water. Calculated volumes of extracts were slowly passed through cartridges based on DOC concentration. Cartridges were rinsed with pH 2 ultrapure water and dried with N2 gas. Samples were finally eluted with methanol into precombusted amber glass vials, dried with N2 gas and stored at −20 °C until DOM analysis. The extracts were continuously injected into the standard ESI source with a flow rate of 2 μl min−1 and an ESI capillary voltage of 3.5 kV in negative ion mode. One hundred single scans with a transient size of 4 mega word (MW) data points, an ion accumulation time of 0.3 s, and within the mass range of m/z 150–1200, were co-added to a spectrum with absorption mode for phase correction, thereby resulting in a resolving power of 750,000 (FWHM at m/z 400). All FT-ICR mass spectra were internally calibrated using organic matter homologous series separated by 14 Da (-CH2 groups). The mass measurement accuracy was typically within 1 ppm for singly charged ions across a broad m/z range (150–1200 m/z).Data Analysis software (BrukerDaltonik v4.2) was used to convert raw spectra to a list of m/z values using FT-MS peak picker with a signal-to-noise ratio (S/N) threshold set to 7 and absolute intensity threshold to the default value of 100. Putative chemical formulae were assigned using the software Formularity (v1.0)48 following the Compound Identification Algorithm49. In total, 19,538 molecular formulas were putatively assigned for all samples (n = 300) based on the following criteria: S/N  > 7, and mass measurement error  0.80, P ≤ 0.001; Fig. S9). Similar conclusions were also obtained with either OTUs or genera when relating the pairwise distances of molecular traits with SparCC correlation coefficient ρ values among DOM molecules in Fig. 4c. To reduce type I errors in the correlation calculations created by low-occurrence genera or molecules, the majority rule was applied; that is, we retained genera or molecules that were observed in more than half of the total samples (≥75 samples) in China or Norway. The filtered table, including 1340 and 1246 DOM molecules, and 75 and 49 bacterial genera in China and Norway, respectively, was then used for pairwise correlation calculation of DOM and bacteria using SparCC with default parameters35.Finally, bipartite network analysis at a molecular level was performed to quantify the specialization of DOM-bacteria networks (Box 1). The specialization considers interaction abundance and is standardised to account for heterogeneity in the interaction strength and species richness, which describes the levels of “vulnerability” of DOM molecules and “generality” of bacterial taxa27. The threshold correlation for inclusion in bipartite networks was |ρ| = 0.30 to exclude weak interactions and we retained the adjacent matrix with only the interactions between DOM and bacteria. We then constructed two types of interaction networks (i.e., negative and positive networks) based on negative and positive correlation coefficients (SparCC ρ ≤ −0.30 and ρ ≥ 0.30, respectively). According to resource-consumer relationships, negative networks likely indicate the degradation of larger molecules into smaller structures, while positive networks may suggest the production of new molecules via degradation or biosynthetic processes. The SparCC ρ values were multiplied by 10,000 and rounded to integers, and the absolute values were taken for negative networks to enable the calculations of specialization indices. A separate negative and positive sub-network was obtained for each microcosm by selecting the DOM molecules and bacterial taxa in each sample based on its bacterial and DOM compositions. For the network level analysis, we calculated H2′, a measure of specialization27, for each network:$${H}_{2}=-mathop{sum }limits_{i{{mbox{=}}}1}^{i}mathop{sum }limits_{j{{mbox{=}}}1}^{j}{{mbox{(}}}{{{mbox{p}}}}_{{ij}}{{{{{{rm{ln}}}}}}}{{{mbox{p}}}}_{{ij}}{{mbox{)}}}$$
    (2)
    $${H}_{2}{prime} =frac{{H}_{2{max }}{-}{H}_{2}}{{H}_{2{max }}{-}{H}_{2{min }}}$$
    (3)
    where ({{{mbox{p}}}}_{{ij}}{{mbox{=}}}{{{mbox{a}}}}_{{ij}}{{mbox{/}}}m), represents the proportion of interactions in a i × j matrix. ({{{mbox{a}}}}_{{ij}}) is the number of interactions between DOM molecule i and bacterial genus j, which is also referred as “link weight”. m is the total number of interactions between all DOM molecules and bacterial genera. H2′ is the standardised H2 against the minimum (H2min) and maximum (H2max) possible for the same distribution of interaction totals. For the molecular level analysis, we calculated the specialization index Kullback–Leibler distance (d′) for DOM molecules (di′) and bacterial genera (dj′), which describes the levels of “vulnerability” of DOM molecules and “generality” of bacterial genera, respectively:$${d}_{i}=mathop{sum }limits_{j=1}^{j}left(frac{{{{mbox{a}}}}_{{ij}}}{{{{mbox{A}}}}_{i}}{{{mbox{ln}}}}frac{{{{mbox{a}}}}_{{ij}}m}{{{{mbox{A}}}}_{i}{{{mbox{A}}}}_{j}}right)$$
    (4)
    $${d}_{i}{prime} =frac{{d}_{i}-{d}_{{min }}}{{d}_{{max }}-{d}_{{min }}}$$
    (5)
    where ({A}_{i}) = (mathop{sum }limits_{j{{mbox{=}}}1}^{j}{{{mbox{a}}}}_{{ij}}) and ({A}_{j}) = (mathop{sum }limits_{i{{mbox{=}}}1}^{i}{{{mbox{a}}}}_{{ij}}), are the total number of interactions of DOM molecule i and bacterial genus j, respectively. di′ is the standardised di against the minimum (dmin) and maximum (dmax) possible for the same distribution of interaction totals. The equations of dj′ are analogous to di′, replacing j by i. Weighted means of d′ for DOM were calculated for each network as the sum of the product of d′ for each individual molecule i (di′) and relative intensity Ii divided by the sum of all intensities d′  = Ʃ(di′ × Ii)/Ʃ(Ii). Weighted means of d′ for bacteria were calculated as the sum of the d′ of each individual bacterial genus j (dj′) and relative abundance of bacterial genus Ij divided by the sum of all abundance. All calculations were performed using the R package FD V1.0.12. The observed H2′ and d′ values ranged from 0 (complete generalization) to 1 (complete specialization)28 (Fig. S21). Specifically, elevated H2′ or d′ values indicate a high degree of specialization, while lower values suggest increased generalization, that is, higher vulnerability of DOM and/or higher generality of microbes. To directly compare the network indices across the elevations or nutrient enrichment levels, we used a null modelling approach. We standardised the three observed specialization indices (Sobserved; that is, H2′, d′ of DOM, and d′ of bacteria) by calculating their z-scores63 using the equation:$${z}_{S}=({S}_{{{{{{rm{observed}}}}}}}-overline{{{S}}_{{{{{{rm{null}}}}}}}})/({sigma }_{S_{{{{{rm{null}}}}}}})$$
    (6)
    where (overline{{{S}}_{{{{{{rm{null}}}}}}}}) and ({sigma }_{S_{{{{{rm{null}}}}}}}) were, respectively, the mean and standard deviation of the null distribution of S (Snull). One hundred randomised null networks were generated for each bipartite network to derive Snull using the swap.web algorithm, which keeps species richness and the number of interactions per species constant along with network connectance. This null model analysis indicates that interactions between DOM and bacteria were non-random as the observed network specialization indices were generally significantly lower than expected by chance (P  0.05), which tests whether the model structure differs from the observed data, high comparative fit index (CFI  > 0.95) and low standardised root mean squared residual (SRMR  More

  • in

    Assessment of solar radiation resource from the NASA-POWER reanalysis products for tropical climates in Ghana towards clean energy application

    Geography and climatology of study areaThe area of study, Ghana, is on the coastal edge of tropical West African, bounded in latitude 4.5° N and 11.5° N and longitude 3.5° W and 1.5° E, and characterized by a tropical monsoon climate system23,24. Figure 1 shows map of the study area indicating the selected twenty two (22) sunshine measurement stations distributed across the four main climatological zones and Table 1 summarizes the geographical positions of selected stations.Figure 1Adapted from Asilevi27.Map of the study area showing all twenty two (22) synoptic stations distributed in four main climatological zones countrywide.Full size imageTable 1 Geographical position and elevation for study sites.Full size tableAtmospheric clarity over the area is closely connected to cloud amount distribution and rainfall activities, largely determined by the oscillatory migration of the Inter-Tropical Discontinuity (ITD), accounting for the West African Monsoon (WAM)25,26.Owing to the highly variable spatiotemporal distribution of cloud amount vis-à-vis rainfall activities, resulting in contrasting climatic conditions in different parts of the region, the country is partitioned by the Ghana Meteorological Agency (GMet) into four main agro-ecological zones namely, the Savannah, Transition, Forest and Coastal zones as shown in Fig. 123. As a result, the region experiences an estimated Global solar radiation (GSR) intensity peaks in April–May and then in October–November, with the highest monthly average of 22 MJm−2 day−1 over the savannah climatic zone and the lowest monthly average of 13 MJm−2 day−1 over the forest climatic zone27.Research datasetsGround-based measurement dataDaily sunshine duration measurement datasets (n) spanning 1983–2018 where derived for estimating Global solar radiation (GSR). The measurements were taken by the Campbell-Stokes sunshine recorder, mounted at the 22 stations shown in Fig. 1, under unshaded conditions to ensure optimum sunlight exposure. The device concentrates sunlight onto a thin strip of sunshine card, which causes a burnt line representing the total period in hours during which sunshine intensity exceeds 120.0 Wm−2 according to World Meteorological Organization (WMO) recommendations27. The as-received daily records were quality control checked by ensuring 0 ≤ n ≤ N, where N is the astronomical day length representing the possible maximum duration of sunshine in hours determined by Eq. 1 from the latitude (ϕ) of the site of interest and the solar declination (δ) computed by Eq. 227:$$ {text{N}} = frac{2}{15}cos^{ – 1} left[ { – tan phi tan {updelta }} right] $$
    (1)
    $$ {updelta } = 23.45sin left[ {360^{{text{o}}} times frac{{284 + {text{J}}}}{365}} right] $$
    (2)
    where J represents the number for the Julian day of the year (first January is 1 and second January is 2).NASA-POWER Global solar radiation (GSR) reanalysis dataThe satellite-based Global solar radiation (GSR) dataset for specific longitudes and latitudes of all 22 stations, assessed in the study, were retrieved from the National Aeronautics and Space Administration-Prediction of Worldwide Energy Resources (NASA-POWER) reanalysis repository based on the Modern Era Retrospective-Analysis for Research and Applications (MERRA-2) assimilation model products, developed from Surface Radiation Budget, and spanning equal study period (1983–2018). The datasets are accessible on a daily and monthly temporal resolution scales at 0.5° × 0.5° spatial coverage via a user friendly web-based mapping portal: https://power.larc.nasa.gov/data-access-viewer/17. The advantage of the NASA-POWER reanalysis GSR, is the wide spatial coverage, and thus can be used to develop a high spatial resolution of solar radiation across the study area.The POWER Project analyzes, synthesizes and makes available surface radiation related parameters on a global scale, primarily from the World Climate Research Programme (WCRP), Global Energy and Water cycle Experiment (GEWEX), Surface Radiation Budget (SRB) project (Version 2.9), the Clouds and the Earth’s Radiant Energy System (CERES), FLASHFlux (Fast Longwave and Shortwave Radiative Fluxes from CERES and MODIS), and the Global Modeling and Assimilation Office (GMAO)17. Table 2 shows the source satellites and the corresponding temporal coverage used in the development of NASA-POWER GSR products.Table 2 Satellites providing the NASA-POWER GSR datasets20.Full size tableThe monthly average NASA-POWER all-sky shortwave surface radiation reanalysis products are statistically validated, showing reasonable biases of − 6.6–13%, against a global network of surface radiation measurement metadata in an integrated database from the Baseline Surface Radiation Network (BSRN) of the World Radiation Monitoring Center (WRMC)20,22. The datasets are widely used in renewable energy application16,22, agricultural modelling of crop yields28, crop simulation exercises29, and plant disease modelling30.Furthermore, in order to assess the suitability of the NASA-POWER surface solar radiation products for the study area, a synthetic sunshine duration based Global solar radiation (GSR) is developed from the Angstrom-Prescott sunshine duration model by Eq. 3 for comparisons27.$$ {text{GSR}} = left[ {{text{a}} + {text{b}}frac{{text{n}}}{{text{N}}}} right]{text{H}}_{{text{o}}} $$
    (3)
    were Ho (kWhm−2 day−1) is the daily extraterrestrial solar radiation on an horizontal surface, n is the daily sunshine duration measurements obtained from the Ghana Meteorological Agency (GMet), and N is the maximum possible daily sunshine duration or the day length in hours determined by Eq. 1. Generalized regression constants a = 0.25 and b = 0.5 for the study area were determined by Asilevi27 from experimental radiometric data based on correlation regression analysis between atmospheric clarity index (GSR/Ho) and atmospheric cloudlessness index (n/N), for estimating solar radiation over the study area, and compared with other satellite data retrieved from the National Renewable Energy Laboratory (NREL) and the German Aerospace Centre (DLR)27. Ho was calculated from astronomical parameters by Eq. 4:$$ {text{H}}_{0} = frac{{24{ } cdot { }60}}{pi } cdot {text{G}}_{{{text{sc}}}} cdot {text{d}}_{{text{r}}} left[ {omega_{{text{s}}} sin varphi sin delta + cos varphi cos delta sin omega_{{text{s}}} } right] $$
    (4)
    where Gsc is the Solar constant in MJm−2 min−1, dr is the relative Earth–Sun distance in meters (m), (omega_{s}) is the sunset hour angle (angular distance between the meridian of the observer and the meridian whose plane contains the sun), (delta) is the angle of declination in degrees (°) and (varphi) is the local latitude. A detailed presentation of the calculation was published in a previous work27.Statistical assessment analysisFor the purpose of assessing the NASA-POWER derived monthly mean GSR (GSRn) datasets in comparison with the estimated Global Solar Radiation (GSRe) datasets used in this paper, the following deviation and correlation methods in Eqs. 5–11, each showing a complimentary result were used: Standard deviation (({upsigma })), residual error (RE), Root mean square error (RMSE), Mean bias error (MBE), Mean percentage error (MPE), Pearson’s correlation coefficient (r), and Willmott index of agreement (d) for n observations31,32,33,34,35. GSRe, GSRn, and RE represent the estimated GSR, NASA-POWER GSR, and the residual error between GSRe and GSRn respectively. A positive RE indicates that sunshine-based estimated GSR is larger than the NASA-POWER reanalysis dataset, while a negative RE indicates that sunshine-based estimated GSR is smaller than the NASA-POWER reanalysis dataset. The arithmetic mean of any dataset is µ.The standard deviation (({upsigma })) was used to check the upper and lower limits of distribution around the mean deviations between GSRe and GSRn in order to ascertain violations between both datasets33. The RMSE is a standard statistical metric to quantify error margins in meteorology and climate research studies, and by definition is always positive, representing zero in the ideal case, plus a smaller value signifying a good marginal deviation31. The MBE is a good indicator for under-or overestimation in observations, with MBE values closest to zero being desirable. The MPE further indicates the percentage deviation between the GSRe and GSRn individual datasets35.$$ {upsigma } = sqrt {frac{1}{{{text{n}} – 1}}mathop sum limits_{{{text{i}} = 1}}^{{text{n}}} left( {{text{GSR}} – {upmu }} right)^{2} } $$
    (5)
    $$ {text{RE}} = {text{GSR}}_{{text{e}}} – {text{GSR}}_{{text{n}}} $$
    (6)
    $$ {text{RMSE}} = sqrt {frac{1}{{text{n}}}mathop sum limits_{{{text{i}} = 1}}^{{text{n}}} left( {{text{RE}}} right)^{2} } $$
    (7)
    $$ {text{MBE}} = frac{1}{{text{n}}}mathop sum limits_{{{text{i}} = 1}}^{{text{n}}} left( {{text{RE}}} right) $$
    (8)
    $$ {text{MPE}} = frac{1}{{text{n}}}mathop sum limits_{{{text{i}} = 1}}^{{text{n}}} left( {frac{{{text{RE}}}}{{{text{GSR}}_{{text{e}}} }} times 100{text{% }}} right) $$
    (9)
    $$ {text{r}} = frac{{mathop sum nolimits_{{{text{i}} = 1}}^{{text{n}}} left( {{text{GSR}}_{{text{e}}} – {upsigma }_{{text{e}}} } right)left( {{text{GSR}}_{{text{n}}} – {upsigma }_{{text{n}}} } right)}}{{left( {{text{n}} – 1} right){upsigma }_{{text{e}}} {upsigma }_{{text{n}}} }} $$
    (10)
    $$ {text{d}} = 1 – left[ {frac{{mathop sum nolimits_{{{text{i}} = 1}}^{{text{n}}} left( {{text{GSR}}_{{text{e}}} – {text{GSR}}_{{text{n}}} } right)^{2} }}{{mathop sum nolimits_{{{text{i}} = 1}}^{{text{n}}} left( {left| {{text{GSR}}_{{text{e}}} – {text{GSR}}_{{{text{nave}}}} left| + right|{text{GSR}}_{{text{n}}} – {text{GSR}}_{{{text{nave}}}} } right|} right)^{2} }}} right] $$
    (11)
    Further, as with other statistical studies in meteorology36, the Pearson’s correlation coefficient (r) was used to quantify the strength of correlation between GSRe and GSRn. Finally, the Willmott index of agreement (d) commonly used in meteorological literature computed from Eq. 7 is used to assess the degree of GSRe/GSRn agreement34. More

  • in

    The crude oil biodegradation activity of Candida strains isolated from oil-reservoirs soils in Saudi Arabia

    Soil sample collectionSoil samples were collected from three different crude oil reservoirs et al. Faisaliyyah, Al Sina’iyah, and Ghubairah located in Riyadh, Saudi Arabia. Briefly, 400 g of soil samples were collected at 0–10 cm depth, under aseptic conditions. Samples were sieved by 2.5 mm pore size sieves, homogenized, and stored at 4ºC until use.Sources of different hydrocarbonsDifferent samples of crude oil, kerosene, diesel, and used oil were collected in sterile flasks from the tankers of Saudi Aramco Company (Dammam, Saudi Arabia). Additionally, another flask was prepared by mixing 1% of each oil in MSM liquid media to make up the mixed oil. The oil samples were sterilized by Millex® Syringe Filters (Merck Millipore co., Burlington, MA, United States) and stored at 4 °C for further usage.Isolation and identification of fungal speciesThe fungal species in the soil contaminated by crude oil were identified using the dilution method. Briefly, 10% of each soil sample was dissolved in distilled water and vortexed thoroughly. Then, 0.2 ml of each sample was cultured on a sterile PDA plate incubated at 28 °C for three days until the growth of different fungal colonies. Carefully, each colony was isolated, re-cultured on new PDA McCartney bottles of PDA slant, and incubated at 28 °C for three days. The fungi were identified microscopically using standard taxonomic keys based on typical mycelia growth and morphological characteristics provided in the mycological keys54. Besides, the taxonomy of the isolated yeast strains was confirmed by the API 20 C AUX kit (Biomerieux Corp., Marcy-l’Étoile, France) (data not shown). The morphology of pure cultures was tested and identified under a light microscope as described before55.The incidence of each strain was calculated as follows:$$ Incidence ;(% ) = frac{{{text{Number }};{text{of }};{text{samples }};{text{showed }};{text{microbial }};{text{growth}}}}{{{text{Total }};{text{samples}}}} times 100 $$Hydrocarbon tolerance testThe growth rate of isolated strains was tested in a liquid medium of MSM mixed with 1% of either crude oil, used oil, diesel, kerosene, or mixed oil. Furthermore, a control sample of MSM liquid medium without any of the oils tested and all culture media were autoclaved at 121 °C for 30 min. After cooling, 1 ml of each isolate was inoculated with one of the above mixtures and incubated at 25 °C on an orbital shaker. The growth rate was measured every three days for a month for each treatment versus the control. All experiments were performed in triplicates.Scanning electron microscopy (SEM)The morphology of different strains of the isolated fungi was tested by SEM, as previously described56, with some modifications. Briefly, 1 ml of each growing strain, in the liquid media, was centrifuged at the maximum speed (14,000 rpm) for 1 min, followed by fixation with 2.5% glutaraldehyde, and overnight incubation at 5 °C. Later, the sample was pelleted, washed with distilled water, then dehydrated with different ascending concentrations of ethanol (30, 50, 70, 90, 100 (v/v)) for 15 min at room temperature. Finally, samples were examined in the Prince Naif Research Centre (King Saud University, Riyadh, Saudi Arabia) by the JEOL JEM-2100 microscope (JEOL, Peabody, MA, United States), according to the manufacturer instructions.Crude oil degradation assayA modified version of the DCPIP assay57 was employed to assess the oil-degrading ability of the fungal isolates. For each strain, 100 ml of the autoclaved MSM was mixed with 1% (V/V) of one of the hydrocarbons (crude oil, used oil, diesel, kerosene, or mixed oil), 0.1% (v/v) of Tween 80, and 0.6 mg/mL of the redox indicator (DCPIP). Then, 1–2 ml of different fungi growing in liquid media (24–48 h) add to the Crude Oil Degradation media, prepared previously, and incubated for two weeks in a shaking incubator at 25 °C. All flasks were covered and protected from light, aeration, or temperature exchanges to reduce the effects of oil weathering (evaporation, photooxidation). The surfactant Tween 80 was used for bio-stimulation and acceleration of the biosurfactant production by increasing metabolism58. A non-inoculated Crude Oil Degradation media was used as the negative control. Afterward, the colorimetric analysis for the change in DCPIP color was estimated, spectrophotometrically, at 420 nm. All experiments were performed in triplicates.Preparation of cell-free supernatant (CFS)To prepare the Cell-Free Supernatant (CFS), all isolates were grown in MSM broth medium with 1% of either crude oil, used oil, diesel, kerosene, or mixed oil for 30 days in a shaking incubator at 25 °C. After incubation, the cells were removed by centrifugation at 10,000 rpm for 30 min at 4 °C. The supernatant (CFS) was collected and filter-sterilized with a 0.45 μm pore size sterile membrane. CFS was screened for the production of different biosurfactants. All the experiments were carried out in triplicates, and the average values were calculated.Drop-Collapse assayThe Drop-Collapse assay was performed as previously described9, with some modifications. 100 µl of crude oil was applied on glass slides, then 10 µl of each CFS was added to the center of the slide surface and incubated for a minute at room temperature. The slides were imaged by a light microscope using the 10X objective lenses. The spreading on the soil surface was scored by either « + » to indicate the level of positive spreading, biosurfactant production, or «—» for negative spreading. Biosurfactant production was considered positive at the drop diameter ≥ 0.5 mm, compared to the negative control (treated with distilled water).Oil spreading assayAn amount of 20 ml of water was added to the Petri plate (size of 100 mm) and mixed with 20 µl of crude oil or mixed oil, which created a thin layer on the water surface. Then, 10 µl of CFS was delivered onto the surface of the oil, and the clear zone surrounding the CFS drop was observed. The results were compared to the negative control (without CFS) and positive control of 1% SDS41. We have measured the clear zones diameter from images and calculate the actual values in regards to the diameter of the Petri dish (10 cm). The assay was performed in triplicates.Emulsification activity assayThe emulsification activity of each isolate was assessed by mixing equal volumes of MSM broth medium of each isolate with different oils in separate tubes. The samples were homogenized by vortex at high speed for two minutes at room temperature (25 °C) and allowed to settle for 24 h. The tests were performed in duplicate. Then, the emulsification index was calculated as follows59:$$ Emulsification; activity; left( % right) = frac{{{text{Height }};{text{of }};{text{emulsion }};{text{layer}}}}{{{text{Total }};{text{height}}}} times 100 $$Recovery of biosurfactantsThe recovery of biosurfactants from CFS was tested through different assays:Acid precipitation assay3 ml of each CFS was adjusted by 6 N HCl to pH 2 and incubated for 24 h at 4 °C. Later, equal volumes of chloroform/methanol mixture (2:1 v/v) were added to each tube, vortexed, and incubated overnight at room temperature. Afterward, the samples were centrifuged for 30 min at 10,000 rpm (4 °C), the precipitate (Light brown colored paste) was air-dried in a fume hood, and weighed53.Solvent extraction assayThe CFS containing biosurfactant was treated with a mixture of extraction solvents (equal volumes of methanol, chloroform, and acetone). Then, the new mixture was incubated in a shaking incubator at 200 rpm, 30 °C for 5 h. The precipitate was separated into two layers, in which the lower layer (White) was isolated, dried, weighed, and stored60.Ammonium sulfate precipitation assayThe CFS containing biosurfactant was precipitated with 40% (w/v) ammonium sulfate and incubated overnight at 4 °C. The samples were centrifuged at 10,000 rpm for 30 min (4 °C). The precipitate was collected and extracted with an amount of acetone equal to the volume of the supernatant. After centrifugation, the precipitate (Creamy-white) was isolated, air-dried in a fume hood, and weighed53.Zinc sulfate precipitation methodSimilarly, 40% (w/v) zinc sulfate was mixed with the CFS containing biosurfactant. Then, the mixture was incubated at 4 °C, overnight. The precipitate (Light Brown) was collected by centrifugation at 10,000 rpm for 30 min (4 °C), air-dried in a fume hood, and weighed53.Statistical analysisAll experiments were performed in triplicate, and the results were expressed as the mean values ± standard deviation (SD). One-way ANOVA and Dunnett’s tests were used to estimate the significance levels at P  More

  • in

    Social microbiota and social gland gene expression of worker honey bees by age and climate

    Evans, J. D. & Spivak, M. Socialized medicine: individual and communal disease barriers in honey bees. J. Invertebr. Pathol. 103, S62–S72 (2010).PubMed 
    Article 

    Google Scholar 
    Hughes, D. P., Pierce, N. E. & Boomsma, J. J. Social insect symbionts: evolution in homeostatic fortresses. Trends Ecol. Evol. 23, 672–677 (2008).PubMed 
    Article 

    Google Scholar 
    Simone, M., Evans, J. D. & Spivak, M. Resin collection and social immunity in honey bees. Evolution 63, 3016–3022 (2009).CAS 
    PubMed 
    Article 

    Google Scholar 
    Dalenberg, H., Maes, P., Mott, B., Anderson, K. E. & Spivak, M. Propolis envelope promotes beneficial bacteria in the honey bee (Apis mellifera) mouthpart microbiome. Insects 11, 1–12 (2020).Article 

    Google Scholar 
    Poulsen, M., Bot, A. N. M., Nielsen, M. G. & Boomsma, J. J. Experimental evidence for the costs and hygienic significance of the antibiotic metapleural gland secretion in leaf-cutting ants. Behav. Ecol. Sociobiol. 52, 151–157 (2002).Article 

    Google Scholar 
    Rosengaus, R. B., Traniello, J. F. A., Lefebvre, M. L. & Maxmen, A. B. Fungistatic activity of the sternal gland secretion of the dampwood termite Zootermopsis angusticollis. Insect. Soc. 51, 259–264 (2004).Article 

    Google Scholar 
    Kwong, W. K. & Moran, N. A. Gut microbial communities of social bees. Nat. Rev. Microbiol. 14, 374–384 (2016).CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    Maes, P. W., Floyd, A. S., Mott, B. M. & Anderson, K. E. Overwintering honey bee colonies: effect of worker age and climate on the hindgut microbiota. Insects 12, 1–16 (2021).Article 

    Google Scholar 
    Brown, B. P. & Wernegreen, J. J. Deep divergence and rapid evolutionary rates in gut-associated Acetobacteraceae of ants. BMC Microbiol. 16, 140 (2016).PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    Douglas, A. E. The microbial dimension in insect nutritional ecology. Funct. Ecol. 23, 38–47 (2009).Article 

    Google Scholar 
    Kešnerová, L. et al. Gut microbiota structure differs between honeybees in winter and summer. ISME J. 14, 801–814 (2020).PubMed 
    Article 

    Google Scholar 
    Raymann, K., Shaffer, Z. & Moran, N. A. Antibiotic exposure perturbs the gut microbiota and elevates mortality in honeybees. PLoS Biol. 15, 1–22 (2017).Article 
    CAS 

    Google Scholar 
    Anderson, K. E. & Ricigliano, V. A. Honey bee gut dysbiosis: a novel context of disease ecology. Curr. Opin. Insect Sci. 22, 125–132 (2017).PubMed 
    Article 

    Google Scholar 
    Maes, P. W., Rodrigues, P. A. P., Oliver, R., Mott, B. M. & Anderson, K. E. Diet-related gut bacterial dysbiosis correlates with impaired development, increased mortality and Nosema disease in the honeybee (Apis mellifera). Mol. Ecol. 25, 5439–5450 (2016).CAS 
    PubMed 
    Article 

    Google Scholar 
    Miller, D. L., Smith, E. A. & Newton, I. L. G. A bacterial symbiont protects honey bees from fungal disease. bioRxiv https://doi.org/10.1101/2020.01.21.914325 (2020).Article 
    PubMed 
    PubMed Central 

    Google Scholar 
    Motta, E. V. S., Raymann, K. & Moran, N. A. Glyphosate perturbs the gut microbiota of honey bees. Proc. Natl. Acad. Sci. USA 115, 10305–10310 (2018).CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    Grice, E. A. & Segre, J. A. The skin microbiome. Nat. Rev. Microbiol. 9, 244–253 (2011).CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    Corby-Harris, V. et al. Origin and effect of Alpha 2.2 Acetobacteraceae in honey bee larvae and description of Parasaccharibacter apium gen. nov., sp. nov.. Appl. Environ. Microbiol. 80, 7460–7472 (2014).ADS 
    PubMed 
    PubMed Central 
    Article 
    CAS 

    Google Scholar 
    Floyd, A. S. et al. Microbial ecology of european foul brood disease in the honey bee (Apis mellifera): towards a microbiome understanding of disease susceptibility. Insects 11, 1–16 (2020).MathSciNet 
    Article 

    Google Scholar 
    Babendreier, D., Joller, D., Romeis, J., Bigler, F. & Widmer, F. Bacterial community structures in honeybee intestines and their response to two insecticidal proteins. FEMS Microbiol. Ecol. 59, 600–610 (2007).CAS 
    PubMed 
    Article 

    Google Scholar 
    Sabree, Z. L., Hansen, A. K. & Moran, N. A. Independent studies using deep sequencing resolve the same set of core bacterial species dominating gut communities of honey bees. PLoS ONE 7, e41250 (2012).ADS 
    CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    Anderson, K. E. et al. Microbial ecology of the hive and pollination landscape: bacterial associates from floral nectar, the alimentary tract and stored food of honey bees (Apis mellifera). PLoS ONE 8, e83125 (2013).ADS 
    PubMed 
    PubMed Central 
    Article 
    CAS 

    Google Scholar 
    Rokop, Z. P., Horton, M. A. & Newton, I. L. G. Interactions between cooccurring lactic acid bacteria in honey bee hives. Appl. Environ. Microbiol. 81, 7261–7270 (2015).ADS 
    CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    Cox-foster, D. L. et al. A metagenomic survey of microbes in honey bee colony collapse disorder. Science 318, 283–287 (2007).ADS 
    CAS 
    PubMed 
    Article 

    Google Scholar 
    Anderson, K. E., Rodrigues, P. A. P., Mott, B. M., Maes, P. & Corby-Harris, V. Ecological succession in the honey bee gut: shift in lactobacillus strain dominance during early adult development. Microb. Ecol. 71, 1008–1019 (2016).CAS 
    PubMed 
    Article 

    Google Scholar 
    Powell, J. E., Martinson, V. G., Urban-Mead, K. & Moran, N. A. Routes of acquisition of the gut microbiota of the honey bee Apis mellifera. Appl. Environ. Microbiol. 80, 7378–7387 (2014).ADS 
    PubMed 
    PubMed Central 
    Article 
    CAS 

    Google Scholar 
    Zheng, H., Powell, J. E., Steele, M. I., Dietrich, C. & Moran, N. A. Honeybee gut microbiota promotes host weight gain via bacterial metabolism and hormonal signaling. Proc. Natl. Acad. Sci. USA 114, 4775–4780 (2017).CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    Anderson, K. E. et al. Hive-stored pollen of honey bees: many lines of evidence are consistent with pollen preservation, not nutrient conversion. Mol. Ecol. https://doi.org/10.1111/mec.12966 (2014).Article 
    PubMed 
    PubMed Central 

    Google Scholar 
    Ludvigsen, J. et al. Shifts in the midgut/pyloric microbiota composition within a honey bee apiary throughout a season. Microb. Environ. 30, 235–244 (2015).Article 

    Google Scholar 
    Corby-Harris, V., Maes, P. & Anderson, K. E. The bacterial communities associated with honey bee (Apis mellifera) foragers. PLoS ONE 9, e95056 (2014).ADS 
    PubMed 
    PubMed Central 
    Article 
    CAS 

    Google Scholar 
    Münch, D., Kreibich, C. D. & Amdam, G. V. Aging and its modulation in a long-lived worker caste of the honey bee. J. Exp. Biol. 216, 1638–1649 (2013).PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    Amdam, G. V. Social context, stress, and plasticity of aging. Aging Cell 10, 18–27 (2011).CAS 
    PubMed 
    Article 

    Google Scholar 
    Haddad, L. S., Kelbert, L. & Hulbert, A. J. Extended longevity of queen honey bees compared to workers is associated with peroxidation-resistant membranes. Exp. Gerontol. 42, 601–609 (2007).CAS 
    PubMed 
    Article 

    Google Scholar 
    Robinson, G. E. Hormonal and genetic control of honeybee division of labour. Behav. Physiol. Bees 14–27 (1991).Anderson, K. E. et al. The queen gut refines with age: longevity phenotypes in a social insect model. bioRxiv https://doi.org/10.1101/297507 (2018).Article 

    Google Scholar 
    Amdam, G. V., Norberg, K., Hagen, A. & Omholt, S. W. Social exploitation of vitellogenin. Proc. Natl. Acad. Sci. 100, 1799–1802 (2003).ADS 
    CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    Jones, B., Shipley, E. & Arnold, K. E. Social immunity in honeybees—density dependence, diet, and body mass trade-offs. Ecol. Evol. 8, 4852–4859 (2018).PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    Alaux, C., Ducloz, F., Crauser, D. & Le Conte, Y. Diet effects on honeybee immunocompetence. Biol. Lett. 6, 562–565 (2010).PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    Ohashi, K., Natori, S. & Kubo, T. Expression of amylase and glucose oxidase in the hypopharyngeal gland with an age-dependent role change of the worker honeybee (Apis mellifera L.). Eur. J. Biochem. 265, 127–133 (1999).CAS 
    PubMed 
    Article 

    Google Scholar 
    Vannette, R. L., Mohamed, A. & Johnson, B. R. Forager bees (Apis mellifera) highly express immune and detoxification genes in tissues associated with nectar processing. Sci. Rep. 5, (2015).Ohashi, K., Natori, S. & Kubo, T. Change in the mode of gene expression of the hypopharyngeal gland cells with an age-dependent role change of the worker honeybee Apis mellifera L.. Eur. J. Biochem. 249, 797–802 (1997).CAS 
    PubMed 
    Article 

    Google Scholar 
    Huang, Z. Y. & Robinson, G. E. Regulation of honey bee division of labor by colony age demography. Behav. Ecol. Sociobiol. 39, 147–158 (1996).Article 

    Google Scholar 
    Vojvodic, S. et al. The transcriptomic and evolutionary signature of social interactions regulating honey bee caste development. Ecol. Evol. 5, 4795–4807 (2015).PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    Ohashi, K. et al. Functional flexibility of the honey bee hypopharyngeal gland in a dequeened colony. Zool. Sci. 17, 1089–1094 (2000).CAS 
    Article 

    Google Scholar 
    Harwood, G., Salmela, H., Freitak, D. & Amdam, G. Social immunity in honey bees: royal jelly as a vehicle in transferring bacterial pathogen fragments between nestmates. J. Exp. Biol. 224 (2021).Santos, K. S. et al. Profiling the proteome complement of the secretion from hypopharyngeal gland of Africanized nurse-honeybees (Apis mellifera L.). Insect. Biochem. Mol. Biol. 35, 85–91 (2005).CAS 
    PubMed 
    Article 

    Google Scholar 
    Cremer, S., Armitage, S. A. O. & Schmid-Hempel, P. Social immunity. Curr. Biol. 17, 693–702 (2007).Article 
    CAS 

    Google Scholar 
    Mattila, H. R. & Otis, G. W. Dwindling pollen resources trigger the transition to broodless populations of long-lived honeybees each autumn. Ecol. Entomol. 32, 496–505 (2007).Article 

    Google Scholar 
    Crailsheim, K., Riessberger, U., Blaschon, B., Nowogrodzki, R. & Hrassnigg, N. Short-term effects of simulated bad weather conditions upon the behaviour of food-storer honeybees during day and night (Apis mellifera carnica Pollmann). Apidologie 30, 299–310 (1999).Article 

    Google Scholar 
    Ricigliano, V. A. et al. Honey bees overwintering in a southern climate: Longitudinal effects of nutrition and queen age on colony-level molecular physiology and performance. Sci. Rep. 8, 1–11 (2018).CAS 
    Article 

    Google Scholar 
    Ricigliano, V. A. et al. Honey bee colony performance and health are enhanced by apiary proximity to US Conservation Reserve Program (CRP) lands. Sci. Rep. 9, 1–11 (2019).CAS 
    Article 

    Google Scholar 
    Fukuda, H. S. K. Seasonal change of the honey bee worker longevity in Sapporo, North Japan with notes on some factors affecting life span. Ecol. Soc. Jpn. 16, 206–212 (1966).
    Google Scholar 
    Mattila, H. R., Harris, J. L. & Otis, G. W. Timing of production of winter bees in honey bee (Apis mellifera) colonies. Insect. Soc. 48, 88–93 (2001).Article 

    Google Scholar 
    Feliciano-Cardona, S. et al. Honey bees in the tropics show winter bee-like longevity in response to seasonal dearth and brood reduction. Front. Ecol. Evol. 8, 1–8 (2020).Article 

    Google Scholar 
    Döke, M. A., Frazier, M. & Grozinger, C. M. Overwintering honey bees: biology and management. Curr. Opin. Insect. Sci. 10, 185–193 (2015).PubMed 
    Article 

    Google Scholar 
    Liu, C. M. et al. BactQuant: an enhanced broad-coverage bacterial quantitative real-time PCR assay. BMC Microbiol. 12, 56 (2012).CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    Liu, C. M. et al. FungiQuant: a broad-coverage fungal quantitative real-time PCR assay. BMC Microbiol. 12, 1 (2012).CAS 
    Article 

    Google Scholar 
    Schloss, P. D. et al. Introducing mothur: open-source, platform-independent, community-supported software for describing and comparing microbial communities. Appl. Environ. Microbiol. 75, 7537–7541 (2009).ADS 
    CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    Evans, J. D. Beepath: an ordered quantitative-PCR array for exploring honey bee immunity and disease. J. Invertebr. Pathol. 93, 135–139 (2006).CAS 
    PubMed 
    Article 

    Google Scholar 
    Bourgeois, A. L., Rinderer, T. E., Beaman, L. D. & Danka, R. G. Genetic detection and quantification of Nosema apis and N. ceranae in the honey bee. J. Invertebr. Pathol. 103, 53–58 (2010).CAS 
    PubMed 
    Article 

    Google Scholar 
    Pearson, K. Mathematical contributions to the theory of evolution. On a form of spurious correlation which may arise when indices are used in the measurement of organs. Proc. R. Soc. Lond. 60, 489–498 (1986).Gloor, G. B. & Reid, G. Compositional analysis: a valid approach to analyze microbiome high throughput sequencing data. Can. J. Microbiol. 703, 0821 (2016).
    Google Scholar 
    Comas, M. CoDaPack 2.0: a stand-alone, multi-platform compositional software. Options 1–10 (2011).Větrovský, T. & Baldrian, P. The variability of the 16S rRNA gene in bacterial genomes and its consequences for bacterial community analyses. PLoS ONE 8, 1–10 (2013).Article 
    CAS 

    Google Scholar 
    Yek, S. H., Nash, D. R., Jensen, A. B. & Boomsma, J. J. Regulation and specificity of antifungal metapleural gland secretion in leaf-cutting ants. Proc. Biol. Sci. 279, 4215–4222 (2012).PubMed 
    PubMed Central 

    Google Scholar 
    Evans, J. D. et al. Immune pathways and defence mechanisms in honey bees Apis mellifera. Insect. Mol. Biol. 15, 645–656 (2006).CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    Steinmann, N., Corona, M., Neumann, P. & Dainat, B. Overwintering is associated with reduced expression of immune genes and higher susceptibility to virus infection in honey bees. PLoS ONE 10, 1–18 (2015).Article 
    CAS 

    Google Scholar 
    Seehuus, S.-C.C., Norberg, K., Gimsa, U., Krekling, T. & Amdam, G. V. Reproductive protein protects functionally sterile honey bee workers from oxidative stress. Proc. Natl. Acad. Sci. USA 103, 962–967 (2006).ADS 
    CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    Liu, J. R., Yang, Y. C., Shi, L. S. & Peng, C. C. Antioxidant properties of royal jelly associated with larval age and time of harvest. J. Agric. Food Chem. 56, 11447–11452 (2008).CAS 
    PubMed 
    Article 

    Google Scholar 
    Li-E, M., Jia, L., Yan, J., Xiao-Wen, L. & Xin, L. Isolation, purification and characterization of superoxide dismutase from royal jelly of the Italian worker bee, Apis mellifera. Acta Entomol. Sin. 47, 171–177 (2004).
    Google Scholar 
    Bottacini, F. et al. Bifidobacterium asteroides PRL2011 genome analysis reveals clues for colonization of the insect gut. 7, 1–14 (2012).Killer, J., Dubná, S., Sedláček, I. & Švec, P. Lactobacillus apis sp. nov., from the stomach of honeybees (Apis mellifera), having an in vitro inhibitory effect on the causative agents of American and European foulbrood. Int. J. Syst. Evol. Microbiol. 64, 152–157 (2014).Casteels, P. et al. Isolation and characterization of abaecin, a major antibacterial response peptide in the honeybee (Apis mellifera). Eur. J. Biochem. 187, 381–386 (1990).CAS 
    PubMed 
    Article 

    Google Scholar 
    Casteels, P., Ampe, C., Jacobs, F. & Tempst, P. Functional and chemical characterization of hymenoptaecin, an antibacterial polypeptide that is infection-inducible in the honeybee (Apis mellifera). J. Biol. Chem. 268, 7044–7054 (1993).CAS 
    PubMed 
    Article 

    Google Scholar 
    Barke, J. et al. A mixed community of actinomycetes produce multiple antibiotics for the fungus farming ant Acromyrmex octospinosus. BMC Biol. 8, 109 (2010).PubMed 
    PubMed Central 
    Article 
    CAS 

    Google Scholar 
    Lyapunov, Y. E., Kuzyaev, R. Z., Khismatullin, R. G. & Bezgodova, O. A. Intestinal enterobacteria of the hibernating Apis mellifera mellifera L. bees. Microbiology 77, 373–379 (2008).Paiva, C. N. & Bozza, M. T. Are reactive oxygen species always detrimental to pathogens?. Antioxid. Redox Signal. 20, 1000–1034 (2014).CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    Burritt, N. L. et al. Sepsis and hemocyte loss in honey bees (Apis mellifera) Infected with Serratia marcescens strain sicaria. PLoS ONE 11, 1–26 (2016).Article 
    CAS 

    Google Scholar 
    Bae, Y. S., Choi, M. K. & Lee, W. J. Dual oxidase in mucosal immunity and host-microbe homeostasis. Trends Immunol. 31, 278–287 (2010).CAS 
    PubMed 
    Article 

    Google Scholar 
    Ha, E. M., Oh, C. T., Bae, Y. S. & Lee, W. J. A direct role for dual oxidase in Drosophila gut immunity. Science 80(310), 847–850 (2005).ADS 
    Article 
    CAS 

    Google Scholar 
    Crailsheim, K., Hrassnigg, N., Gmeinbauer, R., Szolderits, M. J. & Schneider, L. H. W. Pollen utilization in non-breeding honeybees in Winter. J. Insect. Phys. 39, 369–373 (1993).Article 

    Google Scholar 
    Corona, M. & Robinson, G. E. Genes of the antioxidant system of the honey bee: annotation and phylogeny. 15, 687–701 (2006).Schwarz, R. S., Huang, Q. & Evans, J. D. Hologenome theory and the honey bee pathosphere. Curr. Opin. Insect. Sci. 10, 1–7 (2015).PubMed 
    Article 

    Google Scholar 
    Corona, M., Hughes, K. A., Weaver, D. B. & Robinson, G. E. Gene expression patterns associated with queen honey bee longevity. Mech. Age. Dev. 126, 1230–1238 (2005).CAS 
    Article 

    Google Scholar 
    Santos, D. E., Souza, A. D. O., Tibério, G. J., Alberici, L. C. & Hartfelder, K. Differential expression of antioxidant system genes in honey bee (Apis mellifera L.) caste development mitigates ROS-mediated oxidative damage in queen larvae. 20200173, (2020). More

  • in

    Evidence for a mixed-age group in a pterosaur footprint assemblage from the early Upper Cretaceous of Korea

    Wellnhofer, P. The Illustrated Encyclopedia of Pterosaurs (Crescent Books, 1991).Unwin, D. M. The pterosaurs from deep time (Pi Press, 2005).Witton, M. P. Pterosaurs: Natural History (Anatomy (Princeton University Press, 2013).Book 

    Google Scholar 
    Williams, C. J. et al. Helically arranged cross struts in azhdarchid pterosaur cervical vertebrae and their biomechanical implications. iScience 24, 102338 (2021).PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    Bestwick, J., Unwin, D. M., Butler, R. J. & Purnell, M. A. Dietary diversity and evolution of the earliest flying vertebrates revealed by dental microwear texture analysis. Nat. Commun. 11, 5293 (2020).CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    Ryang, W. H. Characteristics of strike-slip basin formation and sedimentary fills and the Cretaceous small basins of the Korean Peninsula. J. Geo. Soc. Korea 49, 31–45 (2013).CAS 

    Google Scholar 
    Kim, B. G. & Park, B. G. Geological report of the Dongbok sheet (1:50,000) (Geological Survey of Korea, Seoul, 1966).Lee, H., Sim, M. S. & Choi, T. Stratigraphic evolution of the northern part of the Cretaceous Neungju basin South Korea. Geosci. J. 23, 849–865 (2019).CAS 
    Article 

    Google Scholar 
    Paik, I. S., Huh, M., So, Y. H., Lee, J. E. & Kim, H. J. Traces of evaporites in Upper Cretaceous lacustrine deposits of Korea: Origin and paleoenvironmental implications. J. Asian Earth Sci. 30, 93–107 (2007).Article 

    Google Scholar 
    Cohen, K. M., Finney, S. M., Gibbard, P. L. & Fan, J.-X. The ICS international Chronostratigraphic chart. Episodes 36, 199–204 (2013).Article 

    Google Scholar 
    Calvo, J. O. & Lockley, M. G. The first pterosaur tracks from Gondwana. Cretac. Res. 22, 585–590 (2001).Article 

    Google Scholar 
    Kukihara, R. & Lockley, M. G. Fossil footprints from the dakota group (Cretaceous) john martin reservoir, bent county, Colorado: New insights into the paleoecology of the Dinosaur freeway. Cretac. Res. 33, 165–182 (2012).Article 

    Google Scholar 
    Lockley, M. & Schumacher, B. A new pterosaur swim tracks locality from the Cretaceous Dakota Group of eastern Colorado: implications for pterosaur swim track behavior. Fossil Footprints of Western North America. Bull. NM Mus. Nat. Hist. Sci, 365–371 (2014).Smith, R. E., Martill, D. M., Unwin, D. M. & Steel, L. Edentulous pterosaurs from the Cambridge Greensand (Cretaceous) of eastern England with a review of Ornithostoma Seeley, 1871. Proc. Geol. Assoc. (2020).Ibrahim, N., Unwin, D. M., Martill, D. M., Baidder, L. & Zouhri, S. A new pterosaur (Pterodactyloidea: Azhdarchidae) from the Upper Cretaceous of Morocco. PLoS ONE 5, e10875 (2010).PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    Martill, D. M. & Ibrahim, N. An unusual modification of the jaws in cf. Alanqa, a mid-Cretaceous azhdarchid pterosaur from the Kem Kem beds of Morocco. Cretac. Res. 53, 59–67 (2015).Article 

    Google Scholar 
    Jacobs, M. L., Martill, D. M., Ibrahim, N. & Longrich, N. A new species of Coloborhynchus (Pterosauria, Ornithocheiridae) from the mid-Cretaceous of North Africa. Cretac. Res. 95, 77–88 (2019).Article 

    Google Scholar 
    Jacobs, M. L. et al. New toothed pterosaurs (Pterosauria: Ornithocheiridae) from the middle Cretaceous Kem Kem beds of Morocco and implications for pterosaur palaeobiogeography and diversity. Cretac. Res. 110, 104413 (2020).Article 

    Google Scholar 
    McPhee, J. et al. A new ? Chaoyangopterid (Pterosauria: Pterodactyloidea) from the Cretaceous Kem Kem beds of southern Morocco. Cretac. Res. 110, 104410 (2020).Article 

    Google Scholar 
    Martill, D. M. et al. A new tapejarid (Pterosauria, Azhdarchoidea) from the mid-Cretaceous Kem Kem beds of Takmout, southern Morocco. Cretac. Res. 112, 104424 (2020).Article 

    Google Scholar 
    Martill, D. M., Unwin, D. M., Ibrahim, N. & Longrich, N. A new edentulous pterosaur from the Cretaceous Kem Kem beds of south eastern Morocco. Cretac. Res. 84, 1–12 (2018).Article 

    Google Scholar 
    Smith, R. E. et al. Small, immature pterosaurs from the Cretaceous of Africa: implications for taphonomic bias and palaeocommunity structure in flying reptiles. Cretac. Res. 130, 105061 (2022).Article 

    Google Scholar 
    Smith, R. E., Martill, D. M., Kao, A., Zouhri, S. & Longrich, N. A long-billed, possible probe-feeding pterosaur (Pterodactyloidea: ?Azhdarchoidea) from the mid-Cretaceous of Morocco North Africa. Cretac. Res. 118, 104643 (2021).Article 

    Google Scholar 
    Kellner, A. W. A. et al. First complete pterosaur from the Afro-Arabian continent: insight into pterodactyloid diversity. Sci. Rep. 9, 17875 (2019).PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    Elgin, R. A. & Frey, E. A new azhdarchoid pterosaur from the Cenomanian (Late Cretaceous) of Lebanon. Swiss J. Geosci. 104, 21–33 (2011).Article 

    Google Scholar 
    Averianov, A. O., Kurochkin, E. N., Pervushov, E. M. & Ivanov, A. V. Two bone fragments of ornithocheiroid pterosaurs from the Cenomanian of Volgograd Region, southern Russia. Acta Palaeontol. Pol. 50 (2005).Averianov, A. & Kurochkin, E. A new pterosaurian record from the Cenomanian of the Volga region. Paleontol. J. 44, 695–697 (2010).Article 

    Google Scholar 
    Nessov, L. Flying reptiles from the Jurassic and cretaceous of the USSR and significance of their remains for the reconstruction of paleogeographical conditions. Vestn. Leningr. Gos. Univ. Ser. 7, 28 (1990).
    Google Scholar 
    Bakhurina, N. N. & Unwin, D. M. A survey of pterosaurs from the Jurassic and Cretaceous of the former Soviet Union and Mongolia. (1995).Averianov, A. O. New records of azhdarchids (Pterosauria, Azhdarchidae) from the Late Cretaceous of Russia, Kazakhstan, and Central Asia. Paleontol. J. 41, 189–197 (2007).Article 

    Google Scholar 
    Averianov, A. Mid-Cretaceous ornithocheirids (Pterosauria, Ornithocheiridae) from Russia and Uzbekistan. Paleontol. J. 41, 79–86 (2007).Article 

    Google Scholar 
    Huh, M., Paik, I. S., Chung, C. H., Hwang, K. G. & Kim, B. S. Theropod tracks from Seoyuri in Hwasun, Jeollanamdo, Korea: occurrence and paleontological significance. J. Geo. Soc. Korea 39, 461–478 (2003).CAS 

    Google Scholar 
    Huh, M. et al. Well-preserved theropod tracks from the Upper Cretaceous of Hwasun County, southwestern South Korea, and their paleobiological implications. Cretac. Res. 27, 123–138 (2006).Article 

    Google Scholar 
    Lockley, M. G., Huh, M. & Kim, B. S. Ornithopodichnus and pes-only sauropod Trackways from the Hwasun tracksite Cretaceous of Korea. Ichnos 19, 93–100 (2012).Article 

    Google Scholar 
    Hwang, K. G., Huh, M. & Paik, I. S. A unique trackway of small theropod from Seoyu-ri, Hwasun-gun Jeollanam province. J. Geo. Soc. Korea 42, 69–78 (2006).CAS 

    Google Scholar 
    Kim, B. S. & Huh, M. Analysis of the acceleration phase of a theropod dinosaur based on a Cretaceous trackway from Korea. Palaeogeogr. Palaeoclimatol. Palaeoecol. 293, 1–8 (2010).Article 

    Google Scholar 
    Marchetti, L. et al. Defining the morphological quality of fossil footprints. Problems and principles of preservation in tetrapod ichnology with examples from the Palaeozoic to the present. Earth-Sci. Rev. 193, 109–145 (2019).Article 

    Google Scholar 
    Rodríguez-de La Rosa, R. A. Pterosaur tracks from the latest Campanian Cerro del Pueblo formation of southeastern Coahuila. Mexico. Geol. Soc. Spec. Publ. 271, 275–282 (2003).Article 

    Google Scholar 
    Lockley, M. G. & Meyer, C. Crocodylomorph trackways from the Jurassic to early cretaceous of North America and Europe: Implications for Ichnotaxonomy. Ichnos 11, 167–178 (2004).Article 

    Google Scholar 
    Ambroggi, R. & De Lapparent, A. Les empreintes de pas fossiles du Maestrichtien d’Agadir. Notes du Service Géologique du Maroc 10, 43–57 (1954).
    Google Scholar 
    Stokes, W. L. Pterodactyl tracks from the Morrison Formation. J. Paleontol. 31, 952–954 (1957).
    Google Scholar 
    Delair, J. Note on Purbeck fossil footprints, with descriptions of two hitherto unknown forms from Dorset. Proceedings of the Dorset Natural History and Archaeological Society. 92–100 (1963).Hwang, K.-G., Huh, M. I. N., Lockley, M. G., Unwin, D. M. & Wright, J. L. New pterosaur tracks (Pteraichnidae) from the Late Cretaceous Uhangri Formation, southwestern Korea. Geol. Mag. 139, 421–435 (2002).Article 

    Google Scholar 
    Mazin, J.-M. & Pouech, J. The first non-pterodactyloid pterosaurian trackways and the terrestrial ability of non-pterodactyloid pterosaurs. Geobios 58, 39–53 (2020).Article 

    Google Scholar 
    Masrour, M., de Ducla, M., Billon-Bruyat, J.-P. & Mazin, J.-M. Rediscovery of the Tagragra tracksite (Maastrichtian, Agadir, Morocco): Agadirichnus elegans Ambroggi and Lapparent 1954 is Pterosaurian Ichnotaxon. Ichnos 25, 285–294 (2018).Article 

    Google Scholar 
    Wright, J. L., Unwin, D. M., Lockley, M. G. & Rainforth, E. C. Pterosaur tracks from the Purbeck limestone formation of Dorset England. Proc. Geol. Assoc. 108, 39–48 (1997).Article 

    Google Scholar 
    Lockley, M. G. et al. The fossil trackway Pteraichnusis pterosaurian, not crocodilian: Implications for the global distribution of pterosaur tracks. Ichnos 4, 7–20 (1995).Article 

    Google Scholar 
    Billon-Bruyat, J.-P. & Mazin, J.-M. The systematic problem of tetrapod ichnotaxa: the case study of Pteraichnus Stokes, 1957 (Pterosauria, Pterodactyloidae). Geol. Soc. Spec. Publ. 217, 315–324 (2003).Article 

    Google Scholar 
    Pascual Arribas, C. & Sanz Pérez, E. Huellas de Pterosaurios en el grupo Oncala (Soria, España). Pteraichnus palaciei-saenzi, nov. icnosp. Estudios Geol. 56, 73–100 (2000).
    Google Scholar 
    Calvo, M. M., Vidarte, C. F., Fuentes, F. M. & Fuentes, M. M. Huellas de Pterosaurios en la Sierra de Oncala (Soria, España). Nuevas icnoespecies: pteraichnus vetustior, Pteraichnus parvus. Pteraichnus manueli. Celtiberia 54, 471–490 (2004).
    Google Scholar 
    Fuentes Vidarte, C., Meijide Calvo, M., Meijide Fuentes, F. & Meijide Fuentes, M. Pteraichnus longipodus nov. icnosp. en la Sierra de Oncala (Soria, España). Studia Geologica Salmanticensia, 103–114 (2004).Peng, B.-X., Du, Y.-S., Li, D.-Q. & Bai, Z.-C. The first discovery of the early Cretaceous Pterosaur track and its significance in Yanguoxia, Yongjing County, Gansu Province. Earth Sci.-J. China Univ. Geosci. 29, 21–24 (2004).
    Google Scholar 
    Lee, Y.-N., Lee, H.-J., Lü, J. & Kobayashi, Y. New pterosaur tracks from the Hasandong formation (Lower Cretaceous) of Hadong County South Korea. Cretac. Res. 29, 345–353 (2008).Article 

    Google Scholar 
    Lee, Y.-N., Azuma, Y., Lee, H.-J., Shibata, M. & Lü, J. The first pterosaur trackways from Japan. Cretac. Res. 31, 263–273 (2010).Article 

    Google Scholar 
    Chen, R. et al. Pterosaur tracks from the early late cretaceous of Dongyang City, Zhejiang Province China. Geol. Bull. China. 32, 693–698 (2013).CAS 

    Google Scholar 
    Li, Y., Wang, X. & Jiang, S. A new pterosaur tracksite from the Lower Cretaceous of Wuerho, Junggar Basin, China: inferring the first putative pterosaur trackmaker. PeerJ 9, e11361 (2021).PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    Ha, S. et al. Diminutive pterosaur tracks and trackways (Pteraichnus gracilis ichnosp. Nov.) from the lower Cretaceous Jinju formation, Gyeongsang basin. Korea. Cretac. Res. 131, 105080 (2021).Article 

    Google Scholar 
    Sánchez-Hernández, B., Przewieslik, A. G. & Benton, M. J. A reassessment of the Pteraichnus ichnospecies from the early Cretaceous of Soria Province Spain. J. Vertebr. Paleontol. 29, 487–497 (2009).Article 

    Google Scholar 
    Zhou, X. et al. A new darwinopteran pterosaur reveals arborealism and an opposed thumb. Curr. Biol. 31, 2429-2436.e2427 (2021).CAS 
    PubMed 
    Article 

    Google Scholar 
    Lü, J. et al. Dragons of the Skies (recent advances on the study of pterosaurs from China) (Zhejiang Science and Technology Press, 2013).
    Google Scholar 
    Beccari, V. et al. Osteology of an exceptionally well-preserved tapejarid skeleton from Brazil: Revealing the anatomy of a curious pterodactyloid clade. PLoS ONE 16, e0254789 (2021).CAS 
    PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    Lü, J. A new boreopterid pterodactyloid pterosaur from the Early Cretaceous Yixian Formation of Liaoning Province, northeastern China. Acta Geologica Sinica-English Edition 84, 241–246 (2010).Article 

    Google Scholar 
    Bennett, S. C. Terrestrial locomotion of pterosaurs: A reconstruction based on Pteraichnus trackways. J. Vertebr. Paleontol. 17, 104–113 (2010).Article 

    Google Scholar 
    Wang, X. & Lü, J. Discovery of a pterodactylid pterosaur from the Yixian Formation of western Liaoning China. Chin. Sci. Bull. 46, A3–A8 (2001).Article 

    Google Scholar 
    Frey, E. et al. A new specimen of nyctosaurid pterosaur, cf. Muzquizopteryx sp. from the Late Cretaceous of northeast Mexico. Revista mexicana de ciencias geológicas 29, 131–139 (2012).
    Google Scholar 
    Wu, W.-H., Zhou, C.-F. & Andres, B. The toothless pterosaur Jidapterus edentus (Pterodactyloidea: Azhdarchoidea) from the Early Cretaceous Jehol Biota and its paleoecological implications. PLoS ONE 12, e0185486 (2017).PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    Lü, J. et al. The toothless pterosaurs from China. Acta Geol. Sin. 90, 2513–2525 (2016).
    Google Scholar 
    Zhang, X., Jiang, S., Cheng, X. & Wang, X. New Material of Sinopterus (Pterosauria, Tapejaridae) from the Early Cretaceous Jehol Biota of China. An. Acad. Bras. Cienc. 91 (2019).Bestwick, J., Unwin, D. M., Butler, R. J., Henderson, D. M. & Purnell, M. A. Pterosaur dietary hypotheses: A review of ideas and approaches. Biol. Rev. 93, 2021–2048 (2018).PubMed 
    Article 

    Google Scholar 
    Chen, H. et al. New anatomical information on Dsungaripterus weii Young, 1964 with focus on the palatal region. PeerJ 8, e8741 (2020).PubMed 
    PubMed Central 
    Article 

    Google Scholar 
    Li, D. et al. A manus dominated pterosaur track assemblage from Gansu, China: Implications for behavior. Sci. Bull. 60, 264–272 (2015).Article 

    Google Scholar 
    Masrour, M., Pascual-Arribas, C., de Ducla, M., Hernández-Medrano, N. & Pérez-Lorente, F. Anza palaeoichnological site. Late Cretaceous. Morocco. Part I. The first African pterosaur trackway (manus only). J. African Earth Sci. 134, 766–775 (2017).Article 

    Google Scholar 
    Bramwell, C. D. & Whitfield, G. R. Biomechanics of Pteranodon. Phil. Trans. R. Soc. Lond. B. 267, 503–581 (1974).Article 

    Google Scholar 
    Bennett, S. C. Terrestrial locomotion of pterosaurs: a reconstruction based on Pteraichnus trackways. J. Vertebr. Paleontol. 17, 104–113 (1997).Article 

    Google Scholar 
    Mazin, J.-M., Billon-Bruyat, J.-P., Hantzpergue, P. & Lafaurie, G. Ichnological evidence for quadrupedal locomotion in pterodactyloid pterosaurs: Trackways from the Late Jurassic of Crayssac (southwestern France). Geol. Soc. Spec. Publ. 217, 283–296 (2003).Article 

    Google Scholar 
    Henderson, D. M. Pterosaur body mass estimates from three-dimensional mathematical slicing. J. Vertebr. Paleontol. 30, 768–785 (2010).Article 

    Google Scholar 
    Lockley, M. G. & Wright, J. L. Pterosaur swim tracks and other ichnological evidnce of behaviour and ecology. Geol. Soc. Spec. Publ. 217, 297–313 (2003).Article 

    Google Scholar 
    Lockley, M., Mitchell, L. & Odier, G. P. Small Theropod track assemblages from middle Jurassic Eolianites of eastern Utah: Paleoecological insights from dune Ichnofacies in a transgressive sequence. Ichnos 14, 131–142 (2007).Article 

    Google Scholar 
    Fiorillo, A. R., Hasiotis, S. T., Kobayashi, Y. & Tomsich, C. S. A pterosaur manus track from Denali National park, Alaska Range, Alaska United States. Palaios 24, 466–472 (2009).Article 

    Google Scholar 
    Bell, P. R., Fanti, F. & Sissons, R. A possible pterosaur manus track from the late Cretaceous of Alberta. Lethaia 46, 274–279 (2013).Article 

    Google Scholar 
    Stinnesbeck, W. et al. Theropod, avian, pterosaur, and arthropod tracks from the uppermost Cretaceous Las Encinas Formation, Coahuila, northeastern Mexico, and their significance for the end-Cretaceous mass extinction. Geol. Soc. Am. Bull. 129, 331–348 (2017).Article 

    Google Scholar 
    Xing, L. et al. Late Cretaceous ornithopod-dominated, theropod, and pterosaur track assemblages from the Nanxiong Basin, China: New discoveries, ichnotaxonomy, and paleoecology. Palaeogeogr. Palaeoclimatol. Palaeoecol. 466, 303–313 (2017).Article 

    Google Scholar 
    Lockley, M. G., Gierlinski, G. D., Adach, L., Schumacher, B. & Cart, K. Newly discovered tetrapod ichnotaxa from the Upper Blackhawk Formation Utah. Bull. N. M. M. Nat. Hist. Sci. 79, 469–480 (2018).
    Google Scholar 
    Lockley, M. G. & Gillette, D. Pterosaur and bird tracks from a new Late Cretaceous locality in Utah. Verteb. Paleontol. Utah 99, 355–359 (1999).
    Google Scholar 
    Bennett, S. C. The ontogeny of Pteranodon and other pterosaurs. Paleobiology 19, 92–106 (1993).Article 

    Google Scholar 
    Bennett, S. C. Year-classes of pterosaurs from the Solnhofen Limestone of Germany: taxonomic and systematic implications. J. Vertebr. Paleontol. 16, 432–444 (1996).Article 

    Google Scholar 
    Chiappe, L. M., Codorniú, L., Grellet-Tinner, G. & Rivarola, D. Argentinian unhatched pterosaur fossil. Nature 432, 571–572 (2004).CAS 
    PubMed 
    Article 

    Google Scholar 
    Codorniú, L., Chiappe, L. & Rivarola, D. Neonate morphology and development in pterosaurs: evidence from a Ctenochasmatid embryo from the Early Cretaceous of Argentina. Geol. Soc. Spec. Publ. 455, 83–94 (2018).Article 

    Google Scholar 
    Mickelson, D. L., Lockley, M. G., Bishop, J. & Kirkland, J. A New Pterosaur Tracksite from the Jurassic Summerville Formation, near Ferron Utah. Ichnos 11, 125–142 (2004).Article 

    Google Scholar  More