Ecology

AbstractDetecting C4 plants consumption is central to investigating animal ecology, agriculture, dietary transitions, and socio-environmental adaptations, and can be done using carbon isotope analysis. The conventional δ¹³C threshold used to identify C4 plant intake does not consider substantial ecological variability across Europe. By analyzing over 4,000 δ13C values from archaeological C3 and C4 grains, we present a European-wide C3 grain δ13C baseline and establish adjusted δ13C threshold estimations for C4 consumption from the site to the ecozone scale using multicomponent environmental models and ecozone cluster analysis. We show that a fixed threshold lead to under- or overestimation of C4 plant consumption, particularly in northern/humid and southern/arid regions, where the threshold needs to be revised downwards or upwards by up to 2‰. This refined framework offers a more accurate baseline for interpreting human and animal diet and enhances our understanding of the spread, adoption and consumption of C4 crops across Europe.

Similar content being viewed by others

Contracting eastern African C4 grasslands during the extinction of Paranthropus boisei

Article
Open access
30 March 2021

Isotopic evidence of high reliance on plant food among Later Stone Age hunter-gatherers at Taforalt, Morocco

Article
Open access
29 April 2024

The Pleistocene high-elevation environments between 2.02 and 0.6 Ma at Melka Kunture (Upper Awash Valley, Ethiopia) based upon stable isotope analysis

Article
Open access
19 March 2024

IntroductionEstimating the proportion of C3 versus C4 plants in human and animal diet is a key part of bioarchaeological, palaeontological and ecological research. Scholars worldwide have been investigating the spread of millet across Eurasia because this highly nutritious and drought-resistant C4 plant can address various questions about past societies1,2. This includes complex social structures, the adoption of new subsistence strategies, the adaptation to challenging climatic and environmental settings, as well as mobility and health status3,4,5,6,7. In ecology and palaeontology, identifying C3 and C4 diets provides insight into (past) habitats, niche partitioning, and animal behaviours8,9,10,11,12.Stable carbon (C) isotope analyses—paired with nitrogen (N) isotope analyses when investigating collagen—represent the most efficient and preferred method to identify C4 plant consumption from bioarchaeological and palaeontological skeletal remains. Due to different photosynthetic pathways between C3 and C4 plants, the ratio of 13C/12C isotopes (expressed as δ13C in ‰) is significantly more negative in C3 plants (−35 to −23‰) compared to C4 plants (−14 to −10‰)13,14,15. Plant δ13C is transferred to the consumer’s body tissues along the food chain, following well-described and quantified fractionation processes16,17, which leads to enriched δ13C values in the tissues of C4 plant consumers compared to C3 plant consumers. In general, bone or dentine collagen δ13C values above −18‰ and enamel or bone apatite δ13C values above −10‰ indicate a mixed diet of C3 and C4 plants13,18,19. In contrast, δ13C values above −12‰ and −4‰, respectively, represent a pure C4 diet13,18,19.However, specific environmental and climatic settings influence the plant’s δ13C value20. For instance, C3 plant δ13C values are more depleted under oceanic or Mediterranean climates, forest soil, dense canopy, elevated humidity or increased CO2 concentration20,21. Conversely, continental climate, aridity, salinity (including sea-spray effect), elevated temperature or high altitude tend to enrich plant δ¹³C values, which are the basis of the terrestrial food chain20,21. The geographical location of the investigated site is also determinant. On average, skeletal tissues can be 1 to 2‰ lower in high latitudes compared to low latitudes in Europe21,22,23. Although C4 plants react differently to environmental and climatic factors14,24, their δ13C values can vary across latitude as well25,26. This implies that the threshold value for C4 plant consumption needs to be adapted depending on the geographical location and environmental settings to avoid misinterpretations of body tissue isotopic composition. This paper fills this research gap to avoid an over- or underestimation of C4 plants dietary intakes of premodern animals and human communities across Europe.The main research questions addressed here are: (i) In which regions of Europe is it required to apply an environmentally adjusted δ13C threshold value for identifying C3 versus C4 plant consumption? (ii) What is the magnitude of this adjustment? (iii) Can we identify specific biogeographic parameters related to this isotope variability in Europe? This study draws on over 400027 published δ13C values from charred archaeological C3 and C4 grains derived from Isotope-Ratio Mass Spectrometry (IRMS)28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,101. We present an innovative and, to the best of our knowledge, unprecedented ecozone-based model framework that integrates multivariate environmental data to facilitate the identification of C3 versus C4 plant diets in bioarchaeological, palaeontological and ecological research.ResultsEcozone cluster modelUsing multicomponent environmental datasets from topographical and climatic variables, we applied k-means cluster analysis to determine zones of similar environmental conditions across Europe. The model provides 20 spatial clusters including one cluster with unclassified values where not all conditions were equally met (Figs. 1 and  S1, Table 1). Due to numeric quantization of k observations, we labelled the clusters based on the percentiles of the respective data ranges using expressional combinations of temperature, climatic moisture index (CMI), and topography. Some of them are geographically restricted to specific regions, such as the cold and humid ecozones 4 and 5 in North-Eastern Europe, the mild and very humid ecozone 8 at the Atlantic coast, or the hot and arid ecozones in the south of Europe, northern Africa and the Near East (e.g., clusters 10, 16 and 18). Other ecozones are spatially more scattered, for example European high mountain ranges (e.g., clusters 11 and 15). The results can be compared to the biome-based ecoregions from the literature102 despite the reduced number of modelled ecozones.Fig. 1: Site distribution over the European Ecozone clusters.20 clusters based on temperature (TMP), moisture availability (CWB), and topography (DEM) were defined using k-means cluster analysis (with k = 20), including unclassified NA values (i.e., inland water). See the methods and material section for a description of the open source TMP, CWB and DEM data and their provenience. The sites are distributed over 15 clusters. See Table 1 for the ecozone descriptions and numbering and Fig. S1 for the ecozones displayed without sites. Figure by Michael Kempf, created using the open source R and QGIS software.Full size imageTable 1 K-means cluster ecozone cluster summary table including descriptionFull size tableThe sites from which the investigated grains originate cover 15 out of 20 modelled ecozones (Fig. 1, Supplementary data 1), representing most of the European geographical and ecological diversity. However, their distribution is biased by past and modern human activity, archaeological sampling and analytical strategies. Not all ecozones are equally represented in the sample and the ecozones 3, 12, 14, 15 and 18 were excluded from the analyses due to small sample sizes.Isotope diversity among grain speciesDespite the isotopic diversity between the main C3 grain species included in this study (Fig. S2A, B, Supplementary data 2), the two dominant crops of this sample, i.e. the Triticum (wheat, n = 1923) and Hordeum (barley, n = 1843) species, are largely overlapping in the northern and southern parts of Europe (Fig. S2B, Supplementary data 2). Differences from the Central/Western European samples are mostly caused by the wide geographical area represented by this region. In the UK, however, the two crops show notably distinct δ13C values (Fig. S2B, Supplementary data 2), which possibly reflects species-specific differences (e.g., the different timing of the vegetation period, which is earlier for barley, while wheat is more impacted by the summer conditions)28,103 or the environmental diversity of the fields used to grow the different crops29. Because human and animal diet is never based on one single crop, the C3 grain sample was kept as one entity for the rest of the analyses. In contrast, the C3 and C4 grains show distinct δ13C values (Fig. S3C, Supplementary data 2), as expected from their different photosynthetic pathways13,14,15. They are thus considered separately in the rest of the study.Temporal isotope variabilityAmong the entire C3 grains dataset, there is hardly any evolution of δ13C over time (linear model [lm] R² = 0.017, p = 1.75e-17; Pearson’s r value = 0.13, p = <2.2e-16; Fig. S3, Supplementary data 2). When distinguishing between the geographical subsets UK (lm R² = 0.044, p = 0.000276; Pearson’s r value = −0.21, p = 0.000276), Southern (lm R² = 0.03, p = 3.43e-17; Pearson’s r value = 0.17, p = <2.2e-16), Central/Western (lm R² = 0.215, p = 1.80e-32; Pearson’s r value = 0.46, p = <2.2e-16) and Northern Europe (including Denmark: lm R² = 0.078, p = 4.13e-18; Pearson’s r value = 0.28, p = <2.2e-16; and excluding Denmark: lm R² = 0.059, p = 1.55e-09; Pearson’s r value = 0.34, p = <2.2e-16), the positive correlation between C3 grain δ13C values and the grain mean date is particularly weak (Fig. S4A–E, Supplementary data 2). In particular, the slightly stronger relationship observed for Central/Western Europe (Fig. S4B) is biased by the youngest samples from Central France, which exhibit particularly enriched δ13C values (Fig. S4C). At the ecozone level, a weak to moderate and significant increase in C3 grains δ13C values can be observed for ecozone 1 (lm R² = 0.359, p = 5.39e-15; Pearson’s r value = 0.60, p = 5.39e-15) and ecozone 17 only (lm R² = 0.223, p = 0.00157; Pearson’s r value = 0.47, p = 1.57e-03) (Figs. 1 and  S5, Supplementary data 2). The C4 grains dataset has a small sample size and each geographical area is represented by short chronologies, which does not enable any proper diachronic analysis (Fig. S6A). The slight decrease in C4 δ13C values over time might thus be only considered statistically significant in Southern Europe despite the chronological gap of nearly a thousand years between the oldest and youngest cluster (Fig. S6B, Supplementary data 2). This implies that the C3 and the C4 grains datasets were not subdivided into different chronological phases for the subsequent analyses.Geographical isotope variabilitySplitting the C3 grain dataset into geographical subsets (UK, Northern, Southern, and Central/Western Europe) shows that the median δ13C value of C3 grains from Northern Europe is approximately 1‰ lower than that from Southern Europe (Fig. 2A, Table 2). This confirms the previous observations made on different types of samples such as faunal remains22,23 and modern plants21. Yet it has to be stressed that the standard deviation (1 SD) is quite large for both regions (±1.35 and ±1.05, respectively), implying some overlap. Despite the high latitude, C3 grains from the UK exhibit among the highest δ13C values across time (Figs. 2A and  S4), which can be related to the oceanic climate22,23. In Denmark, the low δ13C values of the oldest half of the sample (c. 3700–3000 BCE) shift to particularly enriched δ13C values for the most recent half of the sample (c. 1000 BCE–1000 CE) (Fig. S7, Supplementary data 2). This might reflect changes in agricultural practices and soil management following its decrease in quality starting from the Neolithic period30,31.Fig. 2: Geographical isotopic variability in C3 grains.A C3 grains δ13C values in Europe. B C3 grain δ13C variability compared to latitudinal bins within Europe. The middle line of the box represents the median value, the box is delimited by the quartiles Q1 on the left and Q3 on the right and contains the middle half of the sample, the horizontal lines completed by the outlier dots represent the extent of the data. The mean, median, mean absolute deviation (MAD) and standard deviation (1 SD) values for each region and each latitudinal bin are listed in Table 2. The results of the one-way ANOVA tests related to (A) and (B) and of the Pearson’s correlations related to the C3 grain δ13C versus latitude for the whole dataset and for specific subsets are reported in Supplementary data 2. Figure by Margaux L. C. Depaermentier, created using the open source R software.Full size imageTable 2 Statistical summary for the C3 and C4 grain δ13C values over the latitude bins, regions, modern countries and ecozonesFull size tableUsing the whole dataset, we observe a significant decrease in C3 grain δ13C values with increasing latitude (Fig. 2B, Supplementary data 2). From the median values calculated for each latitudinal bin (Table 2, Fig. 2B), the C3 grains δ13C values from sites above 50° latitude are on average 0.54 to 1.72‰ lower than those of grains from sites at latitudes below 50° (Fig. 2B, Table 2). This confirms the mean offset of around 1–2‰ between Southern and Northern Europe. In comparison, there is a mean variation of 0.46‰ among the median δ13C values of the latitudinal bins above 50° and approximately 0.33‰ among the median δ13C values of the bins below 50°. The difference between southern and northern sites is therefore substantial, yet related to an increasing degree of variability towards the north. When excluding the UK and/or Denmark from this dataset due to the overall elevated values in these regions despite their northern latitude, the decrease in δ13C values with increasing latitude is accordingly even stronger and more significant (Pearson’s r value: −0.25 for the whole sample, −0.27 excluding UK and −0.30 excluding UK and Denmark, with a p-value < 2.2e-16 in each case; see Supplementary data 2).At the modern country level, the C3 grains from Lithuania (median −25.18 ± 1.16‰, n = 153) are on average nearly 2.5‰ lower than those from Jordan (median −22.86 ± 0.74‰, n = 46) (Fig. 3A, Table 2) which exceeds the previously defined offset of 1–2‰ between Southern and Northern Europe21,22,23. On the contrary, the Jordan sample is on average 0.80‰ more enriched than those from Greece (median −23.50 ± 0.82‰, n = 383) or Italy (median −23.50 ± 0.97‰, n = 497) despite their shared southern location. Consequently, the North–South-dichotomy is not enough to characterize the different isotopic composition of grains from diverse parts of Europe and does not account for micro-regional environmental diversity. Moreover, the northernmost countries show standard deviations (SD) of 1.29‰ on average (1.12 to 1.67‰ in total), which is sensibly more than in most of the southern (from 0.40 to 1.43‰; mean: 0.92‰) and of the central/western countries (from 0.39 to 1.51‰, mean: 0.90‰) (see the ANOVA test in Supplementary data 2).Fig. 3: Differences in archaeological charred C3 and C4 grain δ13C values in Europe.A C3 grain δ13C values over the modern countries. B C4 grain δ13C values over modern countries. C The comparison between latitude and C4 grain δ13C values shows a weak but significant negative correlation. Boxplots are defined in Fig. 2. The mean, median, MAD and SD values for each modern country are listed in Table 2. The results of the related one-way ANOVA tests and Pearson’s correlations are available from Supplementary data 2. The red labels show non-representative sample sizes (n < 10). Figure by Margaux L. C. Depaermentier, created using the open source R software.Full size imageSimilarly, the C4 grain δ13C values from Lithuania (median: −10.83 ± 0.48‰, n = 20) are on average nearly 1‰ lower than those from Greece (median: −10.19 ± 0.15‰, n = 12), France (median: −10.15‰, n = 16) or Poland (median: −10.19 ± 0.15‰, n = 9) (Fig. 3B, C; Table 2). The sample sizes from Spain (n = 3, median: −10.69 ± 0.15‰) and from the Czech Republic (n = 1) are too low to be considered representative. This trend confirms previous studies from China with depleted C4 grain δ13C values recorded at higher latitudes25,26. The pattern is further supported by the larger C4 grain δ13C dataset resulting from Alpha Magnetic Spectrometer (AMS) in Europe104, showing generally lower δ13C values in Northern compared to Southern or Central Europe (Fig. S8). However, AMS stable isotopic data lack precision due to differences in calibration compared to IRMS and provide particularly wide and unusual δ13C ranges for C4 grains105. Therefore, these data cannot be used to extend the IRMS dataset in this study. Differences in local plant genotypes are considered more likely triggers for the isotopic variability than climate and water availability24. Both C3 and C4 grains exhibit lower δ13C values in some northern regions relative to the rest of Europe, leading to regional variation in the δ¹³C threshold for identifying C4 consumption.Ecological isotopic variabilityThe geographical isotopic variability is related to environmental factors captured in the ecozone model. C3 grain samples from Lithuania (n = 153), Estonia (n = 11), Finland (n = 44), and from parts of Denmark (n = 86) fall into the subhumid temperate lowlands of North-Eastern Europe represented by ecozone 5 (n = 294). Together with ecozone 17 (n = 42)—represented by grains from Norway only—these samples exhibit the lowest median δ13C values (−25.01 ± 1.25‰ and −25.08 ± 1.20‰, respectively) for charred C3 grains in Europe (Fig. 4, Table 2). The high humidity, low temperature and low to moderately elevated topography of these ecozones can explain the depleted δ13C values20. Ecozone 20 (n = 717), represented by balanced temperate plains scattered over Europe and including samples mostly from Denmark, northern Germany and England, exhibit a much higher median δ13C value (−23.00 ± 1.16‰) and its range hardly overlaps with the other northern samples. Beyond the influence of agricultural practices mentioned above30,31, this can be explained by higher temperatures and moderate humidity characterizing this ecozone. In contrast, the 1007 samples from ecozone 17 and the 140 samples from ecozone 1 show the lowest median δ13C values among the sites below 50° latitude (−23.47 ± 0.82‰ and −23.35 ± 1.41‰, respectively). This reflects the mild and moist conditions of these transition zones at a mid-range altitude. In Southern Europe, the ecozones 7 and 13, representing warm highlands in the Mediterranean area, show the highest median δ13C values (−22.93 ± 1.35‰ and −22.56 ± 0.87‰, respectively), deriving from the drier and warmer climatic conditions. The ecozones 2, 8 and 19 are scattered over wide areas of Europe and are not related to extreme temperatures. Their isotopic ratios show intermediate values (Table 2). Ecozones 3, 12, 14, 15 and 18 cannot be included in this isotope investigation due to their small sample size.Fig. 4: C3 grains δ13C values over the European ecozones clusters.The ecozone numbers refer to the numbering in Table 1. Boxplots are defined in Fig. 2. The boxplots are ordered from top to bottom according to decreasing latitude and to increasing temperature within the ecozone. The red labels underline non-representative sample sizes (n < 10). The mean, median, MAD and SD values for each ecozone can be found in Table 2. The results of the one-way ANOVA test are available from Supplementary data 2. Figure by Margaux L. C. Depaermentier, created using the open source R software.Full size imageDiscussionBuilding on the substantial geographical and ecological variation in isotope values within C3 (and to a lesser extent C4) plants, it is essential to revise the commonly used δ13C threshold for identifying C4 consumption, such as −18.0‰ for mammal collagen (for example, ref. 7). At each investigated site, the C3 and C4 grain δ13C values from this dataset (Fig. 5A) were used to create theoretical collagen δ13C values for a diet based exclusively on these crops (Fig. 5B)—which is no realistic diet for humans or animals and was only used for a first theoretical model. This resulted in site-specific estimations for an overall C3 grain-based diet with 10% to 20% C4 grain inputs (Fig. 5B and Supplementary data 1). Our model shows that at several sites from the Baltic and Nordic countries, human or animal collagen δ13C values below −19.0‰ (with a mean SD of 0.59‰)—and up to −19.68 ± 0.94‰ at Bėlis lake, Lithuania (n = 10), for example—already reflect a low C4 input within a primarily C3-grain-based diet. In the Mediterranean, the same C4 input would result in collagen δ13C values above −17.0 ± 0.77 or −16.0 ± 0.57‰ (Fig. 5B and Supplementary data 1).Fig. 5: Point-based approach for baseline C3 grain δ13C values and estimated threshold values for C4 diet identification in mammal collagen at sites with n ≥ 10 grains.A Median C3 grain δ¹³C values (left) and related SD (right). B Median estimated threshold δ¹³C values for mammal collagen (left) and related SD (right) based on a theoretical 100%-grain-based-diet. The mean, median, SD and MAD values for each site are listed in Supplementary data 1. The same maps including even site with n < 10 are in Fig. S9. Figure by Margaux L. C. Depaermentier, created using the open source R software.Full size imageIt is generally accepted that at least 20% of dietary protein from an alternative source (such as C4 compared to C3 crops) is required to be detected in collagen106. However, with a model based on theoretical grain-based diets, a 20% C4 input produces excessively elevated δ13C estimates for mammal collagen (Supplementary data 1). This highlights the limits of a model based on non-realistic grain-based diets, as despite the important role of grains in human diet, both humans and animals have more varied diets in reality—the latter even consuming mostly other parts of the plant and only rarely grains. And with collagen δ13C values reflecting the protein component of the diet16, plant intake contributes less to collagen isotopic composition than animal-derived proteins, which may alter this threshold value in omnivorous diets. Moreover, in regions where the diet includes significant proportions of mushrooms107, forest-derived foods influenced by the canopy effect108,109, and/or freshwater fish110,111,112, consumers are exposed to more depleted δ13C values. These foods have a much stronger influence on collagen δ13C values than any enrichment from C4 plants. And since such conditions are highly plausible in the Baltics and in Scandinavia, the threshold δ13C value for C4 consumption might be even lower than those suggested by the model with 10% C4 input (Figs. 5A, B and  6). On the contrary, potential marine food consumption17 or sea-spray effects in coastal areas113 need to be considered to avoid any over-estimation of C4 plant intake.Fig. 6: Ecozone-based interpolation of the C3 grain δ13C baseline and of the estimated threshold δ13C values for detecting C4 consumers.A Median δ13C (left) and SD values (right) of charred C3 grains interpolated at the ecozone level. B Median δ13C (left) and SD values (right) of the estimated threshold ranges for C4 consumption for each ecozone. Ecozones 3, 12, 14, 15 and 18 are left grey due to their too small sample sizes. Ecozones 4, 9, 10, 11 and 15 are left grey due to the absence of data. The ecozone’s mean, median, SD and MAD values for C3 grains and for C4 consumption are listed in Table 2 and Tab. S1, respectively. Figure by Michael Kempf, created using the open source R and QGIS software.Full size imageA model involving all other food resources would go beyond the scope of this paper. But to accurately estimate the actual local threshold values for C4 consumption, it is essential to use δ13C and δ15N values from as many local and contemporaneous food resources as possible, including crops (for example32). Wild plants would further represent a better baseline for herbivore’s diets. Yet these are much more seldom in archaeological remains and their isotopic ratios are hardly represented in bioarchaeological studies so far. A model based on modern plants21 or on tree δ13C values114 can thus further serve as comparison δ13C baseline for herbivore’s diets. The combination of various isotope systems113,115,116 and/or the application of mixing models117 are further powerful approaches to disentangle the diverse dietary sources.Considerable isotopic variability is observed within each site, region and ecozone (Figs. 2–4 and 5A, Supplementary data 1–2), indicating that a range rather than a sharp threshold value is more appropriate (Figs. 5B and 6B). This isotopic variability is not only related to the gradual and complex variation of environmental settings across Europe, but is also intrinsic to the grains, as grain δ13C may vary for up to 0.5‰ within one ear despite same species and same growing conditions33. Differing growing conditions, in particular various watering practices, can further impact local grain δ13C values by up to 1.7‰24,28. Wild plants δ13C values may be good comparison references to disentangle anthropogenic and natural differences in water regimes21,118,119. Another variability might be induced by diverging analytical uncertainty between datasets, as the data result from different laboratories and protocols and were obtained with different calibrations—which are relevant information to make sure that the datasets are comparable120. Unfortunately, information on calibration, precision and accuracy of the published isotope data were mostly missing, thus it was not possible to assess the impact on the presented results. Yet when available, the error value for replicated samples was very low and still within the range of analytical errors.An evolution of the grain δ13C values might be expected over time as well, especially when considering the various climatic phases that implied a great variation in moisture and temperature throughout Europe over the investigated time frame. The presented regional differences in grain δ13C values have therefore varied between the climatic phases (Fig. S10), yet the isotopic variations within each region over the various climatic phases is overall significantly weak (Supplementary data 2). A larger dataset for the most recent periods might reveal more effective trends. In this study, the impact of the chronological depth and imbalance covered by this dataset can be considered weak (Figs. S2–S5 and S10).The isotopic difference reported in this study between regions and ecozone therefore remains strong enough to highlight environmentally-driven differences (Figs. 2–7, Table 2) even from grains that were possibly undergoing various cultivation practices. Our approach thus demonstrates that the threshold at −18.0‰ for consumer’s collagen is not universal and is mostly valid in the ecozone 2 (−17.51 ± 1.28‰), whereas the northern ecozone 5 and the widespread ecozone 17 require a threshold δ13C value closer to −19.0‰, i.e., −18.59 ± 1.12 and −18.67 ± 1.08‰, respectively. In high-temperature Mediterranean lowlands (ecozone 16), European plains (ecozone 19), and Atlantic regions (ecozone 20), the threshold shifts to −17.0‰ (i.e., −16.88 ± 0.90‰, −17.15 ± 0.74‰, and −16.71 ± 1.05‰, respectively; Figs. 6, 7, Tab. S1). In the arid southern regions (ecozone 13 [−16.31 ± 0.78‰], and partially ecozones 7 and 16), it approaches −16.0‰ (Fig. 7, Tab. S1). Yet in each case, the SD ranges from 0.74 to 1.26‰, stressing again the isotopic variability within ecozones.Fig. 7: Theoretical δ13C threshold ranges for C4 consumption across the European ecozone clusters.The ecozone numbers refer to the numbering in Table 1. Boxplots are defined in Fig. 1. The boxplots are ordered from top to bottom according to decreasing latitude and to increasing temperature within the ecozone. The red labels underline non-representative sample sizes (n < 10). The mean, median, MAD and SD values for each ecozone are listed in Tab. S1. The purple line represents the revised δ13C threshold value for C4 consumption in mammal collagen (i.e., −18.0‰). Figure by Margaux L. C. Depaermentier, created using the open source R software.Full size imageTo conclude, this paper offers both a European-wide δ13C baseline from archaeological charred C3 grains and a threshold-value-model for environmentally adjusted identification of C4 consumption from the site to the ecozone level across Europe (Fig. 7, Tab. S1). The grain baseline offers the advantage to consider a fundamental dietary resource (in particular for humans) as reference data and can be completed by local foods resources at the site level for more holistic and accurate interpretations. A comparison to wild proxies further enhances the results for animal diets. The threshold estimations are particularly suitable in bioarchaeological, ecological or palaeontological studies for which local plants or food resources are unavailable for calculating an isotopic baseline. However, it requires to account for the great degree of isotopic variability within each geographical entity as underlined by the SD values—a variability slightly increasing with latitude. In this context, the point-based approach (Fig. 5) provides more accurate yet geographically discrete data, while the interpolation-based approach (Fig. 6) offers ecologically sensitive estimates over large areas at the ecozone level, both related to some degree of uncertainty or variability. Future datasets could be used to test (and if necessary adjust) the interpolated values in areas of currently low site density. This innovative and context-sensitive ecozone clustering model based on temperature, humidity, and elevation thus enables more accurate interpretations of both animal ecologies and anthropogenic social and agricultural dynamics across Europe by avoiding over- or underestimation of C4 consumption.Methods and materialIsotopic datasetThe material used in this study consists of published δ13C values from charred grains derived from archaeological context, compiled into one single dataset27. The data was collected from 75 publications until September 202528,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,101, also using the open access online repositories IsoArch121, MAIA122, Isotòpia123, IsoMedIta124 and CIMA125. In total, this represents 4,210 δ13C values of C3 and C4 grains derived from 260 sites dated between 8000 BCE and 1800 CE. The represented C3 plants are oat (Avena species, n = 58), rye (Secale species, n = 325), barley (Hordeum species, n = 1843), and wheat (Triticum species, n = 1923). Broomcorn millet (Panicum miliaceum, n = 57) and foxtail millet (Setaria italica, n = 4) represent the C4 crops from this dataset. To facilitate visualization and pattern-recognition, the C3 and C4 grain datasets were considered separately for the various analyses due to their distinct δ13C values and due to the particularly small C4 grain sample size. C3 and C4 grain δ13C values were then combined to address the question of identifying the introduction of C4 crops in C3 plants-based diets. Despite the fact that anthropogenic agricultural practices such as irrigation can impact grain δ13C values24,28, altering the natural and ecological signal, crops remain an important proxy for human diet regardless of agricultural practices, as their isotopic composition would be transferred to human tissues all the same. The ecological differences are considered important enough across Europe to be detectable from grain isotopic composition despite anthropogenic alterations. Wild plants would have represented a better proxy for animal diet, however, this proxy is lacking from archaeological contexts—or could be derived from tree δ13C values for the most recent periods114.Geographically, the research area spans modern Europe and the Mediterranean countries of the Near East, i.e., between 30 and 63° (N) latitude and between −8 and 45° (E) longitude. Yet the data is not evenly distributed and Denmark is over-represented in terms of sites, while Greece is over-represented in terms of number of samples. There are considerable gaps in several regions of Europe (see Fig. 1 and the isotopic dataset27). Chronologically, this dataset covers most archaeological and historical periods and ranges from 8000 BCE to 1800 CE. In this context, it is important to stress that the oldest samples (8000–6000 BCE) exclusively originate from Greece, whereas samples from Northern Europe are predominantly younger than 1000 BCE—except for some larger site samples dated between 4000 BCE and 2000 BCE. Modern data created in the framework of experimental archaeology were not included in the dataset because of the controlled conditions in which they were produced and because of the different present-day atmospheric composition compared to pre-industrial periods126.Only data obtained from Isotope-Ratio Mass Spectrometry (IRMS) were selected for analyses. Despite the fact that these values are not comparable to IRMS isotopic values due to different calibrations105, a European-wide dataset obtained from Accelerator Mass Spectrometry (AMS) in the context of radiocarbon measurements104 was also used as a comparison dataset to verify whether the same trend is visible among both datasets. Importantly, grain δ13C values are sometimes (i.e., in 12 out of 64 publications used in this study) published in form of corrected values to consider the charring effect on the carbon isotope composition of archaeological C3127,128 and C4 grains129,130. But following the approach by Gron et al.30, we are considering here only uncorrected values in order to enhance the comparability between datasets. This is considered to have no significant impact on this study’s results, as the charring effect is not systematic33 particularly low on grain δ13C values (i.e., 0.06 to 0.18‰33,127,128,129,130 on C3 and C4 grains for a heat up to 300 °C) and remains below analytical errors for isotope analyses131.Statistical analysesAll statistical analyses were performed using R software132 and the results and relevant values are summarised in Supplementary data 2. To determine the relationship between grain δ13C values and chronology or geographical location, we applied both Pearson’s correlation tests and linear models—the latter using the lme4 package133. The Pearson correlation coefficient (Pearson’s r value) indicated the strength and direction of a linear relationship between two tested variables, the confidence interval gives the uncertainty range for the true correlation. While in the linear model, the R-squared values show the percentage of the dataset affected by the relationship and the p-values show the significance of the results. To determine the geographical scale at which significant changes in grain values occur, the dataset was divided into various bins. At the largest scale, the main regions of Europe, split into Northern Europe, Southern Europe, Central/Western Europe and the UK. This is not only convenient but also follows expectations based on previous research21,22,23. At the smallest geographical scale, the dataset was binned according to the borders of modern countries. Boxplots were used for data visualization and one-way ANOVA tests were performed to test the difference in grain δ13C values between the investigated clusters. The results and relevant values for the ANOVA tests are summarised in Supplementary data 2). Because our dataset shows a particularly important isotopic variability, the results for each considered bin or entity/group are presented in the text using the median value (since the mean value is more sensitive to extreme outliers) and the related one standard deviation (1 SD). The tables are showing the mean, median, median absolute deviation (MAD) and 1 SD for each category/bin. In order to integrate an ecological dimension to the investigation of grain δ13C variability, the dataset was also binned into newly determined ecozones, as presented in the section below.Environmental clusterFor the cluster analysis, we used these R-packages: terra134, sf 135,136, gtools137, dplyr138 and ggplot2139, ggspatial140, gridExtra141 for plotting, as described in the R-code provided in the open access repository to this paper142. We used an unsupervised k-means clustering approach based on three spatial predictors to differentiate environmental ecozones across Europe. Components include elevation (DEM), mean temperature (T), and the Climatic Moisture Index (CMI). A DEM derived from the USGS (United States Geological Survey, Global Multi-resolution Terrain Elevation Data 2010, https://earthexplorer.usgs.gov/; last accessed 19th of June 2025). Monthly resolved climate variables for the period 1980-2018 were downloaded from CHELSA143. The CMI represents a standardized water availability index, calculated as the ratio of precipitation (P) to potential evapotranspiration (PET) with CMI = P − PET. We used CHELSA v2.1 monthly CMI data based on the Penman-Monteith equation for PET and downscaled from ERA5 reanalysis. The grids were reprojected to a meter-based projection (EPSG:3857) and cropped to the extent of the geographic European landmass. The dataset (n = 468) were aggregated to a 1000 m resolution using bilinear interpolation prior to clustering. Monthly layers were averaged to create a multiannual mean (1980–2018). To allow comparison between variables with different scales and units, all raster layers were standardized using z-score normalization with mean and standard deviation (sd) of the raster (z = cell value − mean_raster / sd_raster). A regular grid of 1 km spacing was generated across the study area and the centroid of each grid cell was calculated for regular point sampling. To reduce edge effects and avoid NA values near the coast, centroids were restricted to land areas using a simplified buffer around the European landmass boundaries. At each centroid location, the values were extracted from the normalized rasters.
K-Means clusteringK-means clustering was performed using the extracted values, with the determined number of clusters k = 20. This dimensionality was chosen to balance regional ecological resolution with model interpretability, including NA values (replaced by −99999 during the cluster analysis to protect correct geographical raster reassignment). Each centroid was assigned a cluster label based on the combined environmental profile of elevation, temperature, and moisture availability. Cluster labels were spatially joined to centroid coordinates and rasterized back into a continuous spatial layer using the DEM grid as a template. The resulting map classifies observational sites into discrete clustered ecozones (Fig. 1).To characterize each of the 20 ecozone clusters based on the data variability, we summarized the environmental properties of each cluster using the mean and standard deviation of the three input variables: Elevation (ELEV), T, and CMI. All values were normalized using z-scores prior to clustering, enabling direct comparison across variables. The clusters were then qualitatively interpreted based on their relative environmental signatures. For example, clusters with low temperatures, moderately dry conditions and high elevations were categorized as Cool|Moderately Dry|Alpine zones. Clusters with high temperatures and low moisture availability in low areas were described as Very Hot|Very Dry|Moderately Low. Intermediate clusters were labeled based on transitional or temperate climate conditions (Table 1).Determining consumer’s δ
13C valuesTo determine the theoretically expected δ13C values of consumer tissues from a mixed C3-C4 diet, we first associated each C3 grain δ13C value to a measured or assumed C4 grain δ13C value. This means that for each site at which δ13C values of both C3 and C4 grains were available, each C3 grain δ13C value got associated with the mean C4 grain δ13C value of the site. However, most of the sites included in this study provided no C4 grain. In this case, a theoretical C4 grain δ13C value was determined for the site based on the observed values from this dataset. In most regions of Europe, the C4 grain δ13C value is thus expected to be −10‰32,72,87. Yet we observed that the C4 grain δ13C values in Lithuania—and by extension presumably in the northernmost latitudes of Europe—were rather around −11‰68. Similarly, C4 grain δ13C values from the western Mediterranean area seem to be closer to −10.5‰35,47. These regional values were thus used as theoretical C4 grain δ13C values associated with the measured C3 grain δ13C values at each site lacking C4 grains (see summary in Supplementary data 1).In a second step, 5‰ was added to each measured or theoretical grain δ13C value to mimic the fractionation offset that applies between the diet and the consumer’s collagenous tissues after consumption16,17. We applied this to both C3 and C4 grains, resulting in theoretical end-members collagen δ13C values for 100% C3 and 100% C4 based diet, respectively. In a third step, we used these end-members values for each grain to create theoretical collagen δ13C values for a C3-grain-based diet including either 10% or 20% C4 grains. These three first steps were done at the grain level to minimize the loss of resolution and information when working with mean values. In a fourth step, we eventually calculated a mean δ13C value for these two types of diet at the site level (Supplementary data 1). It is fundamental to note that these fictitious diets based on 100% grains are not existing in nature and only represent a theoretical model using grains only.Because a 100% grain-based diet does not exist in nature, the model using 10% of C4 input was considered the most reliable basis for estimating the related collagen δ13C values with low C4 input in a normal mixed diet. These results are presented at the site-level to account for local variability in threshold δ13C values for C4 consumption (Fig. 5). The ecozone cluster model was used to create a map of interpolated threshold δ13C values for C4 consumption and hence suggest environmentally adjusted threshold δ13C values for C4 consumption (Fig. 6C). The European background maps used to create Figs. 5 and 6 are vector map data from https://www.naturalearthdata.com/, implemented using the package rnaturalearth in R-Software132,144.Materials & correspondenceThe corresponding authors are MLCD and MK. The isotopic dataset used in this study is available from the open access repository: Depaermentier, M. L. C. (2025). Isotopic Dataset to: Depaermentier, MLC, Kempf, M, Motuzaitė Matuzevičiūtė, G. “Environmentally adjusted δ13C thresholds for accurate detection of C4 plant consumption in Europe” [Data set]. In Communications Earth & Environment. Zenodo. https://doi.org/10.5281/zenodo.17571650 [ref. 27 in this paper]. The data to reproduce the ecozone clusters is available from this open access repository: Kempf, M. (2025): Related files to: Depaermentier, MLC; Kempf, M; Motuzaitė Matuzevičiūtė, G: Environmentally adjusted δ13C thresholds for accurate detection of C4 plant consumption in Europe (2025) [Data set]. Zenodo. https://doi.org/10.5281/zenodo.15695070 [ref. 147 in this paper]. Climate variables used in this article are freely available from Karger et al. (2017): https://chelsa-climate.org/ (last accessed 19th of June 2025) [ref. 143 in this paper]. The Digital Elevation Model (DEM) can be downloaded from the USGS earthexplorer server: https://earthexplorer.usgs.gov/, last accessed 19th of June 2025.Reporting summaryFurther information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

Data availability

The compiled isotopic dataset used in this study is available from the open access repository: Depaermentier, M. L. C. (2025). Isotopic Dataset to: Depaermentier, MLC, Kempf, M, Motuzaitė Matuzevičiūtė, G. “Environmentally adjusted δ13C thresholds for accurate detection of C4 plant consumption in Europe” [Data set]. In Communications Earth & Environment. Zenodo. https://doi.org/10.5281/zenodo.17571650 [ref. 27 in this paper]. The data to reproduce the ecozone clusters is available from this open access repository: Kempf, M. (2025): Related files to: Depaermentier, MLC; Kempf, M; Motuzaitė Matuzevičiūtė, G: Environmentally adjusted δ13C thresholds for accurate detection of C4 plant consumption in Europe (2025) [Data set]. Zenodo. https://doi.org/10.5281/zenodo.15695070 [ref. 147 in this paper]. Climate variables used in this article are freely available from Karger et al. (2017)143: https://chelsa-climate.org/ (last accessed 19th of June 2025). The Digital Elevation Model (DEM) can be downloaded from the USGS earthexplorer server: https://earthexplorer.usgs.gov/, last accessed 19th of June 2025.
Code availability

The code to reproduce the ecozone clusters is available from this open access repository: Kempf, M. (2025): Related files to: Depaermentier, MLC; Kempf, M; Motuzaitė Matuzevičiūtė, G: Environmentally adjusted δ13C thresholds for accurate detection of C4 plant consumption in Europe (2025) [Data set]. Zenodo. https://doi.org/10.5281/zenodo.15695070 [ref. 147 in this paper].
ReferencesVentresca-Miller, A. R. et al. Adaptability of millets and landscapes: ancient cultivation in North-Central Asia. Agronomy 13, 2848 (2023).Article 
CAS 

Google Scholar 
Motuzaitė Matuzevičiūtė, G. Broomcorn millet: from the past to the future. AFF 2, 177–198 (2025).Article 

Google Scholar 
Hakenbeck, S. E., Evans, J., Chapman, H. & Fothi, E. Practising pastoralism in an agricultural environment: an isotopic analysis of the impact of the Hunnic incursions on Pannonian populations. PloS One 12, e0173079 (2017).Article 

Google Scholar 
Kaupová, S. et al. Dukes, elites, and commoners: dietary reconstruction of the early medieval population of Bohemia (9th–11th Century AD, Czech Republic). Archaeol. Anthropol. Sci. 11, 1887–1909 (2019).Article 

Google Scholar 
Martin, L. et al. The place of millet in food globalization during Late Prehistory as evidenced by new bioarchaeological data from the Caucasus. Sci. Rep. 11, 13124 (2021).Article 
CAS 

Google Scholar 
Sneha, M. L. & Arjun, R. Medicinal knowledge In South India (during neolithic to early historic period): an analysis of staple plant dietary nutrition. CMDR J. Soc. Res. 1, 35–48 (2024).
Google Scholar 
Lightfoot, E., Liu, X. & Jones, M. K. Why move starchy cereals? A review of the isotopic evidence for prehistoric millet consumption across Eurasia. World Archaeol. 45, 574–623 (2013).Article 

Google Scholar 
Drucker, D. G. The isotopic ecology of the mammoth steppe. Annu. Rev. Earth Planet. Sci. 50, 395–418 (2022).Article 
CAS 

Google Scholar 
Drucker, D. G. et al. Ecology of large ungulates in the northeastern Iberian Peninsula during the Upper Palaeolithic through stable isotopes and tooth wear analysis. Quat. Environ. Hum. 2, 100011 (2024).
Google Scholar 
Terry, R. C., Guerre, M. E. & Taylor, D. S. How specialized is a diet specialist? Niche flexibility and local persistence through time of the Chisel-toothed kangaroo rat. Funct. Ecol. 31, 1921–1932 (2017).Article 

Google Scholar 
Saarinen, J., Mantzouka, D. & Sakala, J. Aridity, cooling, open vegetation, and the evolution of plants and animals during the cenozoic. In Nature through Time, edited by E. Martinetto, E. Tschopp & R. A. Gastaldo, pp. 83–107 (Springer International Publishing, Cham, 2020).Prasifka, J. & Heinz, K. The use of C3 and C4 plants to study natural enemy movement and ecology, and its application to pest management. Int. J. Pest Manag. 50, 177–181 (2004).Article 

Google Scholar 
Cerling, T. E., Wang, Y. & Quade, J. Expansion of C4 ecosystems as an indicator of global ecological change in the late Miocene. Nature 361, 344–345 (1993).Article 

Google Scholar 
Farquhar, G. D. On the nature of carbon isotope discrimination in C4 species. Funct. Plant Biol. 10, 205 (1983).CAS 

Google Scholar 
O’Leary, M. H. Carbon isotope fractionation in plants. Phytochemistry 20, 553–567 (1981).Article 

Google Scholar 
Ambrose, S. H. Isotopic Analysis of Paleodiets: Methodological and Interpretative Considerations. In Investigations of ancient human tissue. Chemical analyses in anthropology, edited by M. K. Sandford, pp. 59–130 (Gordon and Breach, Philadelphia, 1993).Lee-Thorp, J. A. On isotopes and old bones. Archaeometry 50, 925–950 (2008).Article 
CAS 

Google Scholar 
Kellner, C. M. & Schoeninger, M. J. A simple carbon isotope model for reconstructing prehistoric human diet. Am. J. Phys. Anthropol. 133, 1112–1127 (2007).Article 

Google Scholar 
Froehle, A. W., Kellner, C. M. & Schoeninger, M. J. Multivariate carbon and nitrogen stable isotope model for the reconstruction of prehistoric human diet. Am. J. Phys. Anthropol. 147, 352–369 (2012).Article 
CAS 

Google Scholar 
van Klinken, G. J., Richards, M. P. & Hedges, R. E. M. An Overview of Causes for Stable Isotopic Variations in Past European Human Populations. Environmental, Ecophysiological, and Cultural Effects. In Biogeochemical approaches to paleodietary analysis. Advances in archaeological and museum science, edited by S. H. Ambrose & M. A. Katzenberg, pp. 39–63 (New York, London, 2002).Cooper, C. G., Cooper, M. D., Richards, M. P. & Schmitt, J. Geographic and seasonal variation in δ13C values of C3 plant arabidopsis: Archaeological implications. J. Archaeol. Sci. 149, 105709 (2023).Article 
CAS 

Google Scholar 
van Klinken, G. J., van der Plicht, J. & Hedges, R. E. M. Bone 13C/12C ratios reflect (palaeo) climatic variations. Geophys. Res. Lett. 21, 445–448 (1994).Article 

Google Scholar 
Hedges, R. E., Stevens, R. E. & Richards, M. Bone as a stable isotope archive for local climatic information. Quat. Sci. Rev. 23, 959–965 (2004).Article 

Google Scholar 
Lightfoot, E. et al. Carbon and nitrogen isotopic variability in foxtail millet (Setaria italica) with watering regime. Rapid Commun. Mass Spectrom. 34, e8615 (2020).Article 
CAS 

Google Scholar 
An, C.-B. et al. Variability of the stable carbon isotope ratio in modern and archaeological millets: evidence from northern China. J. Archaeol. Sci. 53, 316–322 (2015).Article 

Google Scholar 
Dong, Y. et al. The potential of stable carbon and nitrogen isotope analysis of foxtail and broomcorn millets for investigating ancient farming systems. Front. plant Sci. 13, 1018312 (2022).Article 
CAS 

Google Scholar 
Depaermentier, M. L. C. Isotopic dataset to: Depaermentier, MLC, Kempf, M, Motuzaitė Matuzevičiūtė, G. “Environmentally adjusted δ13C thresholds for accurate detection of C4 plant consumption in Europe” (2025).Araus, J. L. et al. Identification of Ancient Irrigation Practices based on the Carbon Isotope Discrimination of Plant Seeds: a Case Study from the South-East Iberian Peninsula. J. Archaeol. Sci. 24, 729–740 (1997).Article 

Google Scholar 
Lightfoot, E. & Stevens, R. E. Stable isotope investigations of charred barley (Hordeum vulgare) and wheat (Triticum spelta) grains from Danebury Hillfort: implications for palaeodietary reconstructions. J. Archaeol. Sci. 39, 656–662 (2012).Article 
CAS 

Google Scholar 
Gron, K. J. et al. Archaeological cereals as an isotope record of long-term soil health and anthropogenic amendment in southern Scandinavia. Quat. Sci. Rev. 253, 106762 (2021).Article 

Google Scholar 
Hald, M. M. et al. Farming during turbulent times: agriculture, food crops, and manuring practices in bronze age to viking age Denmark. J. Archaeol. Sci. Rep. 58, 104736 (2024).
Google Scholar 
Nitsch, E. et al. A bottom-up view of food surplus: using stable carbon and nitrogen isotope analysis to investigate agricultural strategies and diet at Bronze Age Archontiko and Thessaloniki Toumba, northern Greece. World Archaeol. 49, 105–137 (2017).Article 

Google Scholar 
Heaton, T. H., Jones, G., Halstead, P. & Tsipropoulos, T. Variations in the 13C/12C ratios of modern wheat grain, and implications for interpreting data from Bronze Age Assiros Toumba, Greece. J. Archaeol. Sci. 36, 2224–2233 (2009).Article 

Google Scholar 
Aguilera, M., Zech-Matterne, V., Lepetz, S. & Balasse, M. Crop fertility conditions in north-eastern Gaul during the la tène and roman periods: a combined stable isotope analysis of archaeobotanical and archaeozoological remains. Environ. Archaeol. 23, 323–337 (2018).Article 

Google Scholar 
Alagich, R., Gardeisen, A., Alonso, N., Rovira, N. & Bogaard, A. Using stable isotopes and functional weed ecology to explore social differences in early urban contexts: the case of Lattara in mediterranean France. J. Archaeol. Sci. 93, 135–149 (2018).Article 

Google Scholar 
Antanaitis, I. & Ogrinc, N. Chemical analysis of bone: stable isotope evidence of the diet of Neolithic and Bronze Age people In Lithuania. Istorija XLV, 3–12 (2000).
Google Scholar 
Araus, J. L. & Buxó, R. Changes in carbon isotope discrimination in grain cereals from the north-western mediterranean basin during the past seven millenia. Funct. Plant Biol. 20, 117 (1993).CAS 

Google Scholar 
Araus, J. L. et al. Changes in carbon isotope discrimination in grain cereals from different regions of the western Mediterranean Basin during the past seven millennia. Palaeoenvironmental evidence of a differential change in aridity during the late Holocene. Glob. Change Biol. 3, 107–118 (1997).Article 

Google Scholar 
Araus, J. L. et al. Isotope and morphometrical evidence reveals the technological package associated with agriculture adoption in western Europe. PNAS 121, e2401065121 (2024).Article 
CAS 

Google Scholar 
Ben Makhad, S. et al. Crop manuring on the Beauce plateau (France) during the second iron age. J. Archaeol. Sci.: Rep. 43, 103463 (2022).
Google Scholar 
Bernardini, S. et al. New multi-proxy isotopic data on the copper age of eastern Liguria. Riv. di Sci. Preistoriche LXXIII S3, 1037–1043 (2023).
Google Scholar 
Bogaard, A. et al. From traditional farming in morocco to early urban agroecology in northern mesopotamia: combining present-day arable weed surveys and crop isotope analysis to reconstruct past agrosystems in (semi-)arid regions. Environ. Archaeol. 23, 303–322 (2018).Article 

Google Scholar 
Bogaard, A. et al. Crop manuring and intensive land management by Europe’s first farmers. Proc. Natl. Acad. Sci. USA 110, 12589–12594 (2013).Article 
CAS 

Google Scholar 
Cortese, F. et al. Isotopic reconstruction of the subsistence strategy for a Central Italian Bronze Age community (Pastena cave, 2 nd millennium BCE) (2022).DiBenedetto, K. E. Investigating Land Use by the Inhabitants of Western Cyprus During the Early Neolithic (2018).Eklund, M. Changing Agriculture. Stable isotope analysis of charred cereals from Iron Age Öland. (Master thesis, Stockholm University, Stockholm, 2019).Fernández-Crespo, T., Ordoño, J., Bogaard, A., Llanos, A. & Schulting, R. A snapshot of subsistence in Iron Age Iberia: the case of La Hoya village. J. Archaeol. Sci.: Rep. 28, 102037 (2019).
Google Scholar 
Fiorentino, G. et al. Third millennium B.C. climate change in Syria highlighted by Carbon stable isotope analysis of 14C-AMS dated plant remains from Ebla. Palaeogeogr. Palaeoclimatol. Palaeoecol. 266, 51–58 (2008).Article 

Google Scholar 
Fiorentino, G., Caracuta, V., Casiello, G., Longobardi, F. & Sacco, A. Studying ancient crop provenance: implications from δ(13)C and δ(15)N values of charred barley in a Middle Bronze Age silo at Ebla(NW Syria). Rapid Commun. Mass Spectrom. 26, 327–335 (2012).Article 
CAS 

Google Scholar 
García-Collado, M. I. et al. First direct evidence of agrarian practices in the alava plateau (northern Iiberia) during the middle ages through carbon and nitrogen stable isotope analyses of charred seeds. Environ. Archaeol. 1–11 (2022).Gavériaux, F. et al. L’alimentation des premières sociétés agropastorales du Sud de la France: premières données isotopiques sur des graines et fruits carbonisés néolithiques et essais de modélisation. Comptes Rendus. Palevol. (2021).Gavériaux, F., Motta, L., Bailey, P., Brilli, M. & Sadori, L. Crop husbandry at gabii during the iron age and archaic period: the archaeobotanical and stable isotope evidence. Environ. Archaeol. 29, 370–383 (2024).Article 

Google Scholar 
Gillis, R. E. et al. Stable isotopic insights into crop cultivation, animal husbandry, and land use at the Linearbandkeramik site of Vráble-Veľké Lehemby (Slovakia). Archaeol. Anthropol. Sci. 12; https://doi.org/10.1007/s12520-020-01210-2 (2020).Gron, K. J. et al. Nitrogen isotope evidence for manuring of early Neolithic Funnel Beaker Culture cereals from Stensborg, Sweden. J. Archaeol. Sci.: Rep. 14, 575–579 (2017).
Google Scholar 
Halvorsen, L. S., Mørkved, P. T. & Hjelle, K. L. Were prehistoric cereal fields in western Norway manured? Evidence from stable isotope values (δ15N) of charred modern and fossil cereals. Veget. Hist. Archaeobot. 32, 583–596 (2023).Article 

Google Scholar 
Fraser, R. A., Bogaard, A., Schäfer, M., Arbogast, R. & Heaton, T. H. E. Integrating botanical, faunal and human stable carbon and nitrogen isotope values to reconstruct land use and palaeodiet at LBK Vaihingen an der Enz, Baden-Württemberg. World Archaeol. 45, 492–517 (2013).Article 

Google Scholar 
Isaakidou, V. et al. Changing land use and political economy at neolithic and bronze Age Knossos, Crete: stable carbon (δ 13 C) and nitrogen (δ 15 N) isotope analysis of charred crop grains and faunal bone collagen. Proc. Prehist. Soc. 88, 155–191 (2022).Article 

Google Scholar 
Kanstrup, M. When δ15N values reveal manuring practice. Empirical evidence from fieldwork, charring experiments and archaeobotanical remains (Aarhus Universitet, Institut for Agroøkologi, Aarhus, 2012).Karakaya, D. Botanical Aspects of the Environment and Economy at Tell Tayinat from the Bronze to Iron Ages (ca. 2.200–600 BCE), In south-central Turkey (Doctoral Dissertation, Universität Tübingen (Germany), Doctoral Dissertation, Universität Tübingen (Germany, 2020).Karaliūtė, R., Motuzaitė Matuzevičiūtė, G., Styring, A. & Stroud, E. 2600 Years of farming in Eastern Lithuania: soil management according to ancient barley isotopic values. Quaternary Environments and Humans (in preparation).Knipper, C. et al. What is on the menu in a Celtic town? Iron age diet reconstructed at Basel-Gasfabrik, Switzerland. Archaeol. Anthropol. Sci. 9, 1307–1326 (2017).Article 

Google Scholar 
Knipper, C. et al. Reconstructing Bronze Age diets and farming strategies at the early Bronze Age sites of La Bastida and Gatas (southeast Iberia) using stable isotope analysis. PloS One 15, e0229398 (2020).Article 
CAS 

Google Scholar 
Lodwick, L. Cultivating villa economies: archaeobotanical and isotopic evidence for iron age to roman agricultural practices on the chalk downlands of Southern Britain. Eur. j. archaeol. 26, 445–466 (2023).Article 

Google Scholar 
Lodwick, L., Campbell, G., Crosby, V. & Müldner, G. Isotopic evidence for changes in cereal production strategies in iron age and Roman Britain. Environ. Archaeol. 26, 13–28 (2021).Article 

Google Scholar 
Maltas, T., Şahoğlu, V., Erkanal, H. & Tuncel, R. From horticulture to agriculture: New data on farming practices in Late Chalcolithic western Anatolia. J. Archaeol. Sci.: Rep. 43, 103482 (2022).
Google Scholar 
Martínez Sánchez, R. M. et al. Archaeology, chronology, and age-diet insights of two late fourth millennium cal BC pit graves from central southern Iberia (Córdoba, Spain). Int. J. Osteoarchaeol. 30, 245–255 (2020).Article 

Google Scholar 
Messager, E. et al. Archaeobotanical and isotopic evidence of Early Bronze Age farming activities and diet in the mountainous environment of the South Caucasus: a pilot study of Chobareti site (Samtskhe–Javakheti region). J. Archaeol. Sci. 53, 214–226 (2015).Article 
CAS 

Google Scholar 
Minkevičius, K. et al. New insights into the subsistence economy of the Late Bronze Age (1100–400 cal BC) communities in the southeastern Baltic. Archaeol. Balt. 30, 58–79 (2023).Article 

Google Scholar 
Mnich, B. et al. Terrestrial diet in prehistoric human groups from southern Poland based on human, faunal and botanical stable isotope evidence. J. Archaeol. Sci.: Rep. 32, 102382 (2020).
Google Scholar 
Mora-González, A. et al. The isotopic footprint of irrigation in the western Mediterranean basin during the Bronze Age: the settlement of Terlinques, southeast Iberian Peninsula. Veget Hist. Archaeobot. 25, 459–468 (2016).Article 

Google Scholar 
Mora-González, A., Teira-Brión, A., Granados-Torres, A., Contreras-Cortés, F. & Delgado-Huertas, A. Agricultural production in the 1st millennium BCE in Northwest Iberia: results of carbon isotope analysis. Archaeol. Anthropol. Sci. 11, 2897–2909 (2019).Article 

Google Scholar 
Mueller-Bieniek, A. et al. Spatial and temporal patterns in Neolithic and Bronze Age agriculture in Poland based on the stable carbon and nitrogen isotopic composition of cereal grains. J. Archaeol. Sci.: Rep. 27, 101993 (2019).
Google Scholar 
Niekamp, A. N. Crop growing conditions and agricultural practices in bronze age greece: a stable isotope analysis of archaeobotanical remains from Tsoungiza (Master’s thesis, University of Cincinnati, 2016).Nitsch, E. K., Jones, G., Sarpaki, A., Hald, M. M. & Bogaard, A. Farming practice and land management at Knossos, Crete: new insights from δ13C and δ15N analysis of Neolithic and Bronze Age crop remains. In Country in the City: Agricultural functions of protohistoric urban settlements (Aegean and Western Mediterranean), edited by D. Garcia, R. Orgeolet, M. Pomadère & J. Zurbach (Archaeopress, Oxford, 2019), pp. 159–173.O’connell, T. C. et al. Living and dying at the Portus Romae. Antiquity 93, 719–734 (2019).Article 

Google Scholar 
Pate, F. D., Henneberg, R. J. & Henneberg, M. Stable carbon and nitrogen isotope evidence for dietary variability at ancient Pompeii, Italy; https://doi.org/10.5281/zenodo.35526 (2015).Pilaar Birch, S. E. et al. Herd management and subsistence practices as inferred from isotopic analysis of animals and plants at Bronze Age Politiko-Troullia, Cyprus. PloS One 17, e0275757 (2022).Article 
CAS 

Google Scholar 
Piličiauskas, G. et al. The earliest evidence for crop cultivation during the Early Bronze Age in the southeastern Baltic. J. Archaeol. Sci.: Rep. 36, 102881 (2021).
Google Scholar 
Riehl, S., Bryson, R. & Pustovoytov, K. Changing growing conditions for crops during the Near Eastern Bronze Age (3000–1200 BC): the stable carbon isotope evidence. J. Archaeol. Sci. 35, 1011–1022 (2008).Article 

Google Scholar 
Speciale, C. et al. The case study of Case Bastione: first analyses of 3rd millennium cal BC paleoenvironmental and subsistence systems in central Sicily. J. Archaeol. Sci.: Rep. 31, 102332 (2020).
Google Scholar 
Styring, A. K. et al. Urban form and scale shaped the agroecology of early ‘cities’ in northern Mesopotamia, the Aegean and Central Europe. J. Agrar. Change 22, 831–854 (2022).Article 

Google Scholar 
Vaiglova, P. et al. An integrated stable isotope study of plants and animals from Kouphovouno, southern Greece: a new look at Neolithic farming. J. Archaeol. Sci. 42, 201–215 (2014).Article 
CAS 

Google Scholar 
Vaiglova, P. et al. Of cattle and feasts: Multi-isotope investigation of animal husbandry and communal feasting at Neolithic Makriyalos, northern Greece. PloS One 13, e0194474 (2018).Article 

Google Scholar 
Vaiglova, P. et al. Exploring diversity in neolithic agropastoral management in mainland greece using stable isotope analysis. Environ. Archaeol. 28, 62–85 (2021).Article 

Google Scholar 
Vanhanen, S. & Ilves, K. Flax use, weeds and manuring in Viking Age Åland: archaeobotanical and stable isotope analysis. Veget. Hist. Archaeobot. 34, 501–517 (2025).Article 

Google Scholar 
Varalli, A. et al. Bronze Age innovations and impact on human diet: a multi-isotopic and multi-proxy study of western Switzerland. PloS One 16, e0245726 (2021).Article 
CAS 

Google Scholar 
Varalli, A. et al. Insights into the frontier zone of Upper Seine Valley (France) during the Bronze Age through subsistence strategies and dietary patterns. Archaeol. Anthropol. Sci. 15; https://doi.org/10.1007/s12520-023-01721-8 (2023).Wallace, M. P. et al. Stable carbon isotope evidence for neolithic and bronze age crop water management in the Eastern Mediterranean and Southwest Asia. PloS One 10, e0127085 (2015).Article 

Google Scholar 
Lempiäinen-Avci, M. et al. New insight into medieval cultivation at the village of Mankby in Espoo, Finland – comparing stable isotopes of carbon δ¹³C and nitrogen δ¹⁵N of Secale and Hordeum from Mankby to 14th century grain materials from Estonia. In Shattered and Scattered Pasts. Festschrift for Professor Georg Haggrén, edited by T. Heinonen, et al., pp. 68–85 (Waasa Graphics, Vaasa, 2025).Ferrio, J. P., Alonso, N., Voltas, J. & Araus, J. L. Grain weight changes over time in ancient cereal crops: Potential roles of climate and genetic improvement. J. Cereal Sci. 44, 323–332 (2006).Article 

Google Scholar 
Della Penna, V. Tradizione e modernità delle pratiche agricole nei Monti dauni: storia e archeologia dei sistemi agroalimentari subappenninici (Unpublished dissertation, Università degli Studi di Foggia, 2022).Dreslerová, D. et al. Maintaining soil productivity as the key factor in European prehistoric and Medieval farming. J. Archaeol. Sci.: Rep. 35, 102633 (2021).
Google Scholar 
Drtikolová Kaupová, S. et al. Stav izotopových výzkumů stravy, rezidenční mobility a zemědělského hospodaření populace Velké Moravy (9.–10. století). Arch. Rozhl. 74, 203–240 (2022).Article 

Google Scholar 
Hamerow, H. et al. An integrated bioarchaeological approach to the medieval ‘agricultural revolution’: a case study from Stafford, England, c.ad 800–1200. Eur. j. archaeol. 23, 585–609 (2020).Article 

Google Scholar 
Herrscher, E. et al. Dietary practices, cultural and social identity in the Early Bronze Age southern Caucasus. Paléorient, 151–174; (2021).Látková, M., Skála, R. & Drtikolová Kaupová, S. Bioarchaeological characteristics of the wheat (triticum aestivum) consumed at different parts of the early medieval settlement agglomeration of Mikulčice-Kopčany (9th–10th Century AD, Czech Republic). Environ. Archaeol. 30, 267–279 (2025).Article 

Google Scholar 
Reed, K. & Wallace, M. To pretreat, or not to pretreat, that is the question. The value of pretreatment protocols in the stable carbon and nitrogen isotope analysis of archaeobotanical cereal grains from Croatia and Serbia. Sci. Technol. Archaeol. Res. 10; https://doi.org/10.1080/20548923.2024.2410092 (2024).Russell, N., Cook, G. T., Ascough, P., Barrett, J. H. & Dugmore, A. Species specific marine radiocarbon reservoir effect: a comparison of ΔR values between Patella vulgata (limpet) shell carbonate and Gadus morhua (Atlantic cod) bone collagen. J. Archaeol. Sci. 38, 1008–1015 (2011).Article 

Google Scholar 
Schlütz, F. et al. Stable isotope analyses (δ15N, δ34S, δ13C) locate early rye cultivation in northern Europe within diverse manuring practices. Philos. Trans. R. Soc. Lond. Ser. B, Biol. Sci. 380, 20240195 (2025).
Google Scholar 
Treasure, E. R. The frontier of Islam: an archaeobotanical study of agriculture in the Iberian Peninsula (c.700 – 1500 CE) (Unpublished dissertation, Durham University, 2020).Vaiglova, P. Neolithic agricultural management in the Eastern Mediterranean: new insight from a multi-isotope approach (Doctoral dissertation, University of Oxford, 2016).Olson, D. M. et al. Terrestrial ecoregions of the world: a new map of life on earth. BioScience 51, 933 (2001).Article 

Google Scholar 
Larsson, M., Bergman, J. & Olsson, P. A. Soil, fertilizer and plant density: Exploring the influence of environmental factors to stable nitrogen and carbon isotope composition in cereal grain. J. Archaeol. Sci. 163, 105935 (2024).Article 
CAS 

Google Scholar 
Filipović, D. et al. New AMS 14C dates track the arrival and spread of broomcorn millet cultivation and agricultural change in prehistoric Europe. Sci. Rep. 10, 13698 (2020).Article 

Google Scholar 
Vaiglova, P., Lazar, N. A., Stroud, E. A., Loftus, E. & Makarewicz, C. A. Best practices for selecting samples, analyzing data, and publishing results in isotope archaeology. Quat. Int.; https://doi.org/10.1016/j.quaint.2022.02.027 (2023).Hedges, R. E. M. Isotopes and red herrings: comments on Milner et al. and Lidén et al. Antiquity 78, 34–37 (2004).Article 

Google Scholar 
O’Regan, H. J., Lamb, A. L. & Wilkinson, D. M. The missing mushrooms: Searching for fungi in ancient human dietary analysis. J. Archaeol. Sci. 75, 139–143 (2016).Article 

Google Scholar 
Drucker, D. G., Bridault, A., Hobson, K. A., Szuma, E. & Bocherens, H. Can carbon-13 in large herbivores reflect the canopy effect in temperate and boreal ecosystems? Evidence from modern and ancient ungulates. Palaeogeogr. Palaeoclimatol. Palaeoecol. 266, 69–82 (2008).Article 

Google Scholar 
Bonafini, M., Pellegrini, M., Ditchfield, P. & Pollard, A. M. Investigation of the ‘canopy effect’ in the isotope ecology of temperate woodlands. J. Archaeol. Sci. 40, 3926–3935 (2013).Article 

Google Scholar 
Häberle, S. et al. Carbon and nitrogen isotopic ratios in archaeological and modern Swiss fish as possible markers for diachronic anthropogenic activity in freshwater ecosystems. J. Archaeol. Sci.: Rep. 10, 411–423 (2016).
Google Scholar 
Guiry, E. Complexities of stable carbon and nitrogen isotope biogeochemistry in ancient freshwater ecosystems: implications for the study of past subsistence and environmental change. Front. Ecol. Evol. 7; https://doi.org/10.3389/fevo.2019.00313 (2019).Robson, H. K. et al. Carbon and nitrogen stable isotope values in freshwater, brackish and marine fish bone collagen from Mesolithic and Neolithic sites in central and northern Europe. Environ. Archaeol. 21, 105–118 (2016).Article 

Google Scholar 
Göhring, A., Hölzl, S., Mayr, C. & Strauss, H. Identification and quantification of the sea spray effect on isotopic systems in α-cellulose (δ13C, δ18O), total sulfur (δ34S), and 87Sr/86Sr of European beach grass (Ammophila arenaria, L.) in a greenhouse experiment. Sci. Total Environ. 856, 158840 (2023).Article 

Google Scholar 
Büntgen, U. Scrutinizing tree-ring parameters for Holocene climate reconstructions. WIREs Clim. Change 13; https://doi.org/10.1002/wcc.778 (2022).Montgomery, J. et al. Strategic and sporadic marine consumption at the onset of the Neolithic: increasing temporal resolution in the isotope evidence. Antiquity 87, 1060–1072 (2013).Article 

Google Scholar 
Göhring, A., Hölzl, S., Mayr, C. & Strauss, H. Multi-isotope fingerprints of recent environmental samples from the Baltic coast and their implications for bioarchaeological studies. Sci. Total Environ. 874, 162513 (2023).Article 

Google Scholar 
Hopkins, J. B. & Ferguson, J. M. Correction: estimating the diets of animals using stable isotopes and a comprehensive bayesian mixing model. PloS One 7; https://doi.org/10.1371/annotation/d222580b-4f36-4403-bb1f-cfd449a5ed74 (2012).Ferrio, J. P., Araus, J. L., Buxó, R., Voltas, J. & Bort, J. Water management practices and climate in ancient agriculture: inferences from the stable isotope composition of archaeobotanical remains. Veget. Hist. Archaeobot. 14, 510–517 (2005).Article 

Google Scholar 
Jones, P. J., O’Connell, T. C., Jones, M. K., Singh, R. N. & Petrie, C. A. Crop water status from plant stable carbon isotope values: a test case for monsoonal climates. Holocene 31, 993–1004 (2021).Article 

Google Scholar 
Szpak, P., Metcalfe, J. Z. & Macdonald, R. A. Best practices for calibrating and reporting stable isotope measurements in archaeology. J. Archaeol. Sci.: Rep. 13, 609–616 (2017).
Google Scholar 
Salesse, K. et al. IsoArcH.eu: an open-access and collaborative isotope database for bioarchaeological samples from the Graeco-Roman world and its margins. J. Archaeol. Sci.: Rep. 19, 1050–1055 (2018).
Google Scholar 
Farese, M. MAIA: mediterranean archive of isotopic dAta. Pandora https://doi.org/10.48493/55v1-xg54 (2023).Article 

Google Scholar 
Formichella, G., Soncin, S. & Cocozza, C. Isotòpia: a stable isotope database for classical antiquity. Pandora v.1 19.09.2023; https://doi.org/10.48493/m0m0-b436 (2023).Mantile, N., Fernandes, R., Lubritto, C. & Cocozza, C. IsoMedIta: a stable isotope database for Medieval Italy. Res. Data J. Humanit. Soc. Sci. 8, 1–13 (2023).Article 

Google Scholar 
Cocozza, C., Cirelli, E., Groß, M., Teegen, W.-R. & Fernandes, R. Presenting the compendium isotoporum medii aevi, a multi-isotope database for Medieval Europe. Sci. Data 9, 354 (2022).Article 

Google Scholar 
Graven, H., Keeling, R. F. & Rogelj, J. Changes to carbon isotopes in atmospheric CO2 over the industrial era and into the future. Glob. Biogeochem. Cycles 34, e2019GB006170 (2020).Article 
CAS 

Google Scholar 
Nitsch, E. K., Charles, M. & Bogaard, A. Calculating a statistically robust δ 13 C and δ 15 N offset for charred cereal and pulse seeds. Sci. Technol. Archaeol. Res. 1, 1–8 (2015).
Google Scholar 
Stroud, E., Charles, M., Bogaard, A. & Hamerow, H. Turning up the heat: Assessing the impact of charring regime on the morphology and stable isotopic values of cereal grains. J. Archaeol. Sci. 153, 105754 (2023).Article 

Google Scholar 
Teira-Brión, A., Stroud, E., Charles, M. & Bogaard, A. The effects of charring on morphology and stable carbon and nitrogen isotope values of common and foxtail millet grains. Front. Environ. Archaeol. 3; https://doi.org/10.3389/fearc.2024.1473593 (2024).Varalli, A., D’Agostini, F., Madella, M., Fiorentino, G. & Lancelotti, C. Charring effects on stable carbon and nitrogen isotope values on C4 plants: Inferences for archaeological investigations. J. Archaeol. Sci. 156, 105821 (2023).Article 
CAS 

Google Scholar 
Styring, A. K. et al. Recommendations for stable isotope analysis of charred archaeological crop remains. Front. Environ. Archaeol. 3; https://doi.org/10.3389/fearc.2024.1470375 (2024).R. Core Team. R: A language and environment for statistical (R Foundation for Statistical Computing, Vienna, 2021).Bates, D., Mächler, M., Bolker, B. & Walker, S. Fitting linear mixed-effects models using lme4. J. Stat. Soft. 67; https://doi.org/10.18637/jss.v067.i01 (2015).Hijmans, R. J. _terra: Spatial Data Analysis_ (2024).Pebesma, E. & Bivand, R. Spatial Data Science (Chapman and Hall/CRC, New York, 2023).Pebesma, E. Simple Features for R: standardized support for spatial vector data. R. J. 10, 439 (2018).Article 

Google Scholar 
Warnes, G. et al. _gtools: Various R Programming Tools_ (2023).Wickham, H., François, R., Henry, L., Müller, K. & Vaughan, D. _dplyr: A Grammar of Data Manipulation_ (2023).Wickham, H. Ggplot2. Elegant graphics for data analysis (Springer Science+Business Media, LLC, New York, NY, 2016).Dunnington, D. _ggspatial: Spatial Data Framework for ggplot2_ (2023).Auguie, B. _gridExtra: Miscellaneous Functions for “Grid” Graphics (2017).Kempf, M. Environmental data to: Depaermentier, MLC; Kempf, M; Motuzaitė Matuzevičiūtė, G: environmentally adjusted δ13C thresholds for accurate detection of C4 plant consumption in Europe (2025).Karger, D. N. et al. Climatologies at high resolution for the earth’s land surface areas. Sci. Data 4, 170122 (2017).Article 

Google Scholar 
South, A., Michael, S. & Massicotte, P. rnaturalearthdata: world vector map data from natural earth used in ‘rnaturalearth’ (2017).Download referencesAcknowledgementsM.L.C.D. and G.M.M. were funded by the European Union with a Consolidator Grant awarded to Giedrė Motuzaitė Matuzevičiūtė (ERC-CoG, MILWAYS, 101087964). Views and opinions expressed are those of the authors only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them. MK’s research is funded by the Swiss National Science Foundation (SNSF/SNF): Project EXOCHAINS − Exploring Holocene Climate Change and Human Innovations across Eurasia (SNSF grant number: TMPFP2_217358).Author informationAuthor notesThese authors contributed equally: Margaux L. C. Depaermentier, Michael Kempf.Authors and AffiliationsFaculty of History, Vilnius University, Vilnius, LithuaniaMargaux L. C. Depaermentier & Giedrė Motuzaitė MatuzevičiūtėDepartment of Environmental Sciences, University of Basel, Basel, SwitzerlandMichael KempfAuthorsMargaux L. C. DepaermentierView author publicationsSearch author on:PubMed Google ScholarMichael KempfView author publicationsSearch author on:PubMed Google ScholarGiedrė Motuzaitė MatuzevičiūtėView author publicationsSearch author on:PubMed Google ScholarContributionsConceptualization: M.L.C.D., M.K. and G.M.M. Isotope data collection and formal analyses: M.L.C.D. Environmental and cluster analyses: M.K. Writing: M.L.C.D and M.K. Editing: M.L.C.D., M.K., G.M.M. Visualisation: M.K. and M.L.C.D. Revision: M.L.C.D., M.K.Corresponding authorsCorrespondence to
Margaux L. C. Depaermentier or Michael Kempf.Ethics declarations

Competing interests
The authors declare no competing interests.

Peer review

Peer review information
Communications Earth & Environment thanks Ashley McCall and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Primary Handling Editors: Nicola Colombo and Aliénor Lavergne. [A peer review file is available].

Additional informationPublisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.Supplementary informationTransparent Peer Review fileSupplementary InformationDescription of Additional Supplementary FilesSupplementary data 1Supplementary data 2Reporting summaryRights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
Reprints and permissionsAbout this articleCite this articleDepaermentier, M.L.C., Kempf, M. & Motuzaitė Matuzevičiūtė, G. Environmentally adjusted δ13C thresholds for accurate detection of C4 plant consumption in Europe.
Commun Earth Environ 6, 1021 (2025). https://doi.org/10.1038/s43247-025-03031-4Download citationReceived: 13 July 2025Accepted: 12 November 2025Published: 18 December 2025Version of record: 18 December 2025DOI: https://doi.org/10.1038/s43247-025-03031-4Share this articleAnyone you share the following link with will be able to read this content:Get shareable linkSorry, a shareable link is not currently available for this article.Copy shareable link to clipboard
Provided by the Springer Nature SharedIt content-sharing initiative More

AbstractAnthropogenic stressors, including ocean acidification (OA), ocean warming (OW), and salinity changes, are rapidly altering marine ecosystems, with Arctic regions being particularly vulnerable. This study investigates the combined effects of these stressors on the fertilization success of the green sea urchin (Strongylocentrotus droebachiensis) from Kongsfjorden, Svalbard. We exposed gametes to various levels of pH, temperature, and salinity to assess their individual and combined impacts on fertilization performance. Our results show that temperature and pH significantly influenced fertilization success, with temperature having the strongest effect, while salinity had no significant impact. A significant statistical interaction between temperature and pH indicated that warming enhanced fertilization more effectively at higher pH levels, while low pH suppressed this increase. To compare the relative influence of each stressor, we used a conceptual model based on standardized slopes, which supported temperature as the dominant driver, followed by pH. These findings highlight the importance of considering the effects of combined stressors when assessing marine organism responses to climate change, especially in polar ecosystems. Our study underscores the need for further research into the mechanisms driving these combined effects, given that Arctic ecosystems face accelerated environmental changes.

Similar content being viewed by others

Regional sensitivity patterns of Arctic Ocean acidification revealed with machine learning

Article
Open access
19 April 2022

Direct and latent effects of ocean acidification on the transition of a sea urchin from planktonic larva to benthic juvenile

Article
Open access
01 April 2022

Impacts of hypoxic events surpass those of future ocean warming and acidification

Article

11 January 2021

IntroductionAnthropogenic stressors such as ocean acidification (OA) and ocean warming (OW) are negatively impacting marine ecosystems worldwide1. Polar ecosystems, and particularly Arctic ecosystems, are among the most sensitive ecosystems to these stressors2,3,4. Arctic marine ecosystems are experiencing some of the fastest rates of warming and acidification globally, along with a continued decline in sea-ice volume5. These adverse effects are being exacerbated by additional stressors already present in coastal ecosystems, such as pollution or freshening6.Due to their sensitivity to environmental stressors, early-life stages and their associated processes (e.g. fertilization, metamorphosis) are considered critical bottlenecks in the life cycle of marine benthic organisms7. Understanding how future environmental changes will impact them is therefore crucial. Extensive research has been conducted on fertilization success in broadcast spawners under future scenarios of global change, particularly focusing on acidification conditions (reviewed in8). The findings suggest that sea urchin fertilization is robust to small changes in pH9,10,11,12. However, potential negative impacts are anticipated at more extreme levels13.The bulk of this research has concentrated on tropical and temperate species, leaving polar species – especially Arctic ones – relatively understudied. Furthermore, the majority of Arctic research has focused on single-stressor responses such as OA, while it has become evident that a multiple-stressor approach is more suitable to uncover complex combined effects14. Multiple-stressor studies are therefore essential to fully understand the extent to which future scenarios of global change will impact marine ecosystems.Future environmental changes in the Arctic are expected to alter the dynamics of benthic and coastal communities, particularly around Svalbard. Understanding how these communities will fare is vital for informing management practices and policy development. In this context, we investigated the effects of pH, temperature and salinity on the fertilization success of the green sea urchin (Strongylocentrotus droebachiensis) collected from Kongsfjorden, a glacial fjord in the Svalbard archipelago. This species has a circumboreal distribution, extending from temperate to Arctic regions, and plays a central role in the benthic ecosystem of Kongsfjorden15. The population in this fjord seems to be increasing, particularly in the outer reaches, where sea urchins have grazed down significant portions of kelp biomass16. Known for its utility in developmental and environmental biology, S. droebachiensis is well-suited for laboratory studies dues to its ease of collection, maintenance, and spawning.Kongsfjorden is one of the most extensively studied Arctic fjord systems and serves as a reference site for marine science and monitoring17. The unique diversity and abundance of its fauna make it an important early indicator of changes associated to global-change stressors, such as acidification or warming17. Furthermore, ongoing monitoring efforts in the area document temporal changes in environmental parameters (e.g.18,19).In this study, we examined the fertilization of S. droebachiensis gametes under current and extreme levels of seawater pH, temperature, and salinity. Additionally, we assessed the relative contributions of each stressor by developing an index based on performance curves for individual stressors, aiming to discern the mechanisms underlying the influence of these environmental stressors on fertilization in this sea urchin population.Materials and methodsAnimal collection and acclimationAdult green sea urchins were collected in June 2021 by scuba divers from Kongsfjorden, Svalbard (78°59’05.0”N, 11°57’55.0”E; Fig. 1) at depths of 3 to 10 m. Although biometric measurements were not taken, all individuals were sexually mature, as indicated by the presence of ripe gametes and successful spawning following KCl injection. The urchins were transported to the Ny-Ålesund marine laboratory (Ny-Ålesund, Svalbard) acclimated in flow-through tanks supplied with ambient Kongsfjorden seawater (temperature ~ 2 °C, salinity ~ 34, and pH ~ 8.1) and fed ad libitum with seaweed for one week before the experiments commenced.Fig. 1Location of the sea urchin collection site in Kongsfjorden (diving site at Hansneset, Blomstrand), Ny-Ålesund with location of the marine laboratory (experimental site), and the Kb3 station (used for reference environmental data; Fransson et al., 2016). Figure created with ArcMap 10.8.2. and Adobe Illustrator 2022. Map basis and cartography by the Norwegian Polar Institute, hydrographic data from the Norwegian Mapping Authority.Full size imageFertilization assays and experimental designSpawning was induced by means of intracoelomic injection of 0.1 M KCl and gametes were collected separately for quality check. Eggs were collected in glass vials with filtered seawater and inspected for viability (shape and colour). Sperm was collected dry and kept on ice, and quickly inspected for viability (activation when in contact with seawater). Gametes of two females and four males were selected and pooled prior to fertilization assays to reduce individual variability and represent the population-level response.A pilot fertilization assay (five concentrations, four replicates) was carried out to establish the sperm concentration that would lead to a fertilization success of approximately 50% to optimize the chances of observing effects20. Details and results of this pilot study can be found in supplementary material (S1).The main fertilization assays took place in small 24 multi-well plates and involved testing four pH levels (nominal pH 7.1, 7.4, 7.7 and 8.1), four salinities (28, 30, 32 and 34) and three temperatures (1.5 °C, 2.2 °C and 5.3 °C). The selected levels of pH, salinity, and temperature were chosen to reflect both current and projected variability in Kongsfjorden. Reference conditions (pH 8.1, salinity 34, temperature 2.2 °C) correspond to average field measurements18, while altered conditions represent extreme but plausible future scenarios under ongoing climate change19. The experimental design included a factorial combination of the three temperatures and four salinities at pH ~ 8.1, as well as a factorial combination of the same three temperatures and four pHs at ambient salinity of 34, resulting in a total of 21 treatments. Each treatment was replicated six times, and the controls were replicated 12 times, adding up to a total of 144 experimental units (Table 1).Table 1 Treatments used for the fertilization assay.Full size tableFor the assay, each individual well was filled with 2 ml of filtered seawater (FSW, 0.2 μm) from the corresponding treatment. Subsequently, approximately 200 eggs were added and left to sink to the bottom of the well. After 5 min, 1 µl of sperm dilution (at a concentration of 3.9 × 106 sperm/ml) was added in each well and the plates were closed and incubated in a temperature-controlled water bath at the target temperature for 15 min. After that time, one drop of 4% paraformaldehyde (PFA) was added to each well to stop the fertilization process and fix the eggs. The fertilization success in each well was subsequently assessed under a microscope, with the appearance of the perivitelline membrane indicating successful fertilization (see supporting data in21).Preparation of treatments and measurementsThe pH treatments were prepared by bubbling pure CO2 into FSW and mixing manually until reaching target values. For pH, the target values were 7.1, 7.4, 7.7 and 8.1. The pH measurements were taken manually with a pH probe (Hanna Instruments, model HI-98190) calibrated with NIST buffers (LabChem, Zelie-nople, PA, USA), and the final pH values for the treatments were pHNBS 8.04 ± 0.01, pHNBS 7.62 ± 0.04, pHNBS 7.36 ± 0.06 and pHNBS 7.14 ± 0.04. The salinity treatments were prepared in advance by diluting FSW with tap water, until reaching four final salinities: 28, 30, 32 and 34. The temperature conditions were achieved and maintained by means of water baths inside controlled temperature rooms (5.28 °C ± 0.02, 2.18 °C ± 0.05 and 1.55 °C ± 0.06). All physical parameters for the different treatments were measured and recorded immediately before the start of each assay (Table 2). The carbonate chemistry parameters of the experimental water were calculated using the CO2sys package22 and are reported in Table 2, using as input the pH measurements of the experimental water and values of total alkalinity (TA) obtained from18 for Kongsfjorden. Table 2 Mean seawater parameters in the different treatments during the fertilization experiment.Full size tableStatistical analysesAll statistical analyses were performed with R v. 4.0.3. and RStudio v.1.4.110323. A linear model (package “stats”) was used to test for effects of pH, temperature, salinity and their interactions on the fertilization success. The model and final data obtained are reported in the results section. The exploratory models and validation plots are detailed in Supplementary 2.Calculation of relative contributions and stress indexTo quantify the relative contribution of each driver to fertilization success, we fitted simple linear regressions for each individual variable (pH, salinity, and temperature), using subsets of the data where the other two variables were held constant at reference levels (pH 8.1, temperature 2.2 °C, and salinity 34; see Supplementary 3). These reference conditions were chosen to represent typical environmental values at Kb3, the closest oceanographic monitoring station to our sea urchin sampling site in Kongsfjorden18 (Fig. 1).The slope of each regression was used to calculate a comparative index of contribution, standardized relative to a one-unit change in pH (Supplementary 3 and supporting data in24). Specifically, the contribution of temperature (CT) relative to pH was calculated as the ratio between the slope of the pH model and the slope of the temperature model (slope_pH / slope_temp). Similarly, the contribution of salinity (CS) was derived from the ratio between the slope of the pH model and the slope of the salinity model (slope_pH / slope_sal). This approach allowed us to express each stressor’s contribution in standardized units for conceptual comparison. The relative stress index for each variable was calculated individually and used as input for the total stress index (TS). In this context, ‘stress’ is defined as the departure from reference conditions for a given environmental driver (e.g., deviation from pH 8.1), weighted by the relative effect size (slope) estimated from the simple linear models. The stress related to pH (SPH) was calculated by subtracting the reference pH value (pHREF = 8.1) from the pH of each observation (pHO). The stress related to temperature (ST) was calculated by subtracting the reference temperature value (TREF = 2.2°C) from the temperature for each observation and dividing the result by the calculated contribution value for the temperature, according to the following formula:ST = [(TO – TREF)/ CTWhere: ST = stress related to temperature, TO = temperature in each observation,TREF = reference temperature (2.2 °C) and CT = contribution due to the temperature.The salinity-related stress was calculated by subtracting the reference salinity value (SREF = 34) from the salinity for each observation and dividing the result by the calculated contribution value for the salinity, according to the following formula:SS = – [(SO – SREF)/ CSWhere: SS = stress related to salinity, SO = salinity in each observation, SREF = reference salinity (34) and CS = contribution due to the salinity.Finally, the total stress (TS) was calculated additively, according to the following formula:TS = SpH + ST + SS.Where: TS = Total stress; SpH = stress related to pH; ST = Stress related to temperature and SS = stress related to salinity.The fitted regressions and subsequent calculations are reported in the result section (Fig. 3), the model details and verification plots are in Supplementary 3.ResultsEffects of temperature, salinity and pH on the fertilization successFertilization rates across all treatments ranged from 22.5% (at pH 7.1, salinity 34 and temperature 1.5 °C) to 94.3% (at pH 8.1, salinity 28 and temperature 5.3 °C; Fig. 2; and supplementary data in21).The analysis revealed that fertilization success was significantly influenced by pH (negative effect) and temperature (positive effect), as well as their interactions (p < 0.05 and p < 0.001, respectively, Fig. 2, Supplementary 2). In contrast, salinity alone and its interaction with temperature did not significantly affect fertilization (p = 0.96 and p = 0.52, respectively, Fig. 2, Supplementary 2).Fig. 2Fertilization success (%) per temperature (vertical facets), salinity (horizontal facets) and pH (x-axis).Full size imageRelative contribution of environmental drivers to the fertilization successAt reference conditions of pH 8.1 and salinity 34 (Fig. 3a), fertilization success at the three tested temperatures ranged from 37.9% to 88.2% at 1.5 °C, which was similar but less than the upper value at 5.3 °C. At the reference conditions of pH 8.1 and temperature 2.2 °C (Fig. 3b), fertilization rates across different salinities varied from 23.2% (at salinity 28) to 84.2% (at salinity 32). Under reference salinity (34) and temperature (2.2 °C), fertilization success increased with pH, ranging from 30.5% (at pH 7.1) to 79.6% (at pH 8.1) (Fig. 3c). These modelled relationships are derived from simple linear regression models fitted separately for each driver under fixed reference conditions. They are intended to illustrate the relative influence of temperature, salinity, and pH on fertilization under standardized scenarios (Fig. 3). Full model coefficients and diagnostics are provided in Supplementary 3.Fig. 3Fertilization success (%) at the reference conditions according to (a) temperature, (b) salinity and (c) pH.Full size imageFor each pH unit, the calculated relative contribution of temperature to fertilization success was CT = 8.29 °C, while the relative contribution of salinity was CS = 10.71 salinity units. These stress values reflect the standardized magnitude of deviation from reference conditions for each variable, scaled by its relative contribution to fertilization success as estimated by slope values from the individual linear regressions. These values were incorporated into the calculation of the total stress (TS) index, which is shown for each fertilization point in relation to fertilization success (Fig. 4). Further details on the total stress values for each observation are provided in Supplementary 3 and in24.Fig. 4Calculated stress (CS) index based on the relative contribution of all environmental variables.Full size imageDiscussionOur investigations revealed that temperature and pH, both individually and in combination, significantly influenced fertilization success in S. droebachiensis. Specifically, fertilization success generally increased with higher temperatures and was reduced under lower pH conditions. Importantly, the interaction between pH and temperature was significant: fertilization increased more strongly with temperature at high pH, but remained low under low pH conditions even when temperature was elevated. This suggests that warming alone may not fully offset the negative effects of ocean acidification. In contrast, salinity did not have a statistically significant effect on fertilization success, either alone or in combination with temperature. These trends are supported by the fitted interaction model (Supplementary 2) as visualized in Fig. 2.Fertilization success in sea urchins is influenced by temperature through enhanced sperm motility and by gamete sensitivity to pH. Higher temperatures have been shown to increase sperm swimming speed and metabolic activity, which can improve sperm–egg encounter rates and fertilization success25. In contrast, reduced pH can impair intracellular pH regulation in sperm, disrupt the acrosome reaction, and reduce motility by affecting dynein ATPase activity, which is essential for flagellar movement26,27. Further, dynein ATPase activity has been observed to increase linearly with intracellular alkalinity in Strongylocentrotus purpuratus, clearly linking internal pH (pHi) control and motility28. More recent evidence from ”kina” Evechinus chloroticus demonstrates that ocean acidification causes a failure of sperm to maintain pHi, leading to declines in the proportion of motile sperm and swimming speed, which likely reduce fertilization success under near-future acidification scenarios29. These physiological impairments under low pH conditions can lower the proportion of motile sperm and slow swimming velocity, ultimately decreasing fertilization success. Such mechanisms are consistent with our observed results, where warming alone could not fully compensate for the negative effects of low pH on fertilization outcomes.Simplified performance curves and combined stressorsPerformance curves for temperature, pH and salinity are rarely linear and are rather unimodal. However, these curves have a linear part within their tolerance range. We assumed simplified linear performance curves for fertilization success within the tested range of pH, salinity and temperature used in our experiments (Fig. 3). The significant interaction observed between pH and temperature suggests that their combined impact, while generally following the direction of their individual effect, is not mathematically additive consistent from non-linear effects expected from a unimodal performance curve and/or stressors interactions. These findings underscore the importance of experiments covering a wide range of variability to capture the full performance curve.Relative contribution of individual driversOur results are in line with previous studies that investigated fertilization success in sea urchins. Our regression analysis showed that temperature has a more substantial impact on fertilization success compared to pH. This interpretation is supported by the fitted interaction model (Supplementary 2) and the visualized patterns in Fig. 2, which show differential temperature effects across pH levels.The main bulk of research on this topic points towards a general consensus in that fertilization in sea urchins is resilient to low pH in the range expected for OA (reviewed in8), although some studies have shown contrasting results, even though these seem to be due to differences in methodology. For example, in the tropical shortspined sea urchin (Heliocidaris erythrogramma), exposure of gametes to acidification (− 0.4 pH units) negatively affected fertilization success (observed reductions of 24%30). Similarly, acidification (800 and 1800 ppm CO2) was found to decrease fertilization success at lower sperm concentrations in the temperate red sea urchin (Strongylocentrotus franciscanus)31. In the sea urchins Hemicentrotus pulcherrimus and Echinometra mathaei, fertilization rate decreased with increasing pCO2 concentrations, although the impact differed between individual females32.To standardize the relative effect size of each environmental driver, we extracted slope coefficients from simple linear regressions fitted separately for temperature, pH, and salinity under reference conditions (Supplementary 3). The regressions indicated that fertilization success increased by approximately 2.8% per 1 °C increase in temperature and by 2.3% per 0.1-unit decrease in pH. Expressed differently, a 10% change in fertilization success would correspond to roughly a 3.5 °C increase in temperature or a 0.4-unit decrease in pH. These values provide a more comparable basis for interpreting the relative strength of each driver. Nevertheless, they should not be taken as predictive relationships, since changes of this magnitude may not occur with equal ecological probability. Instead, they offer an approximate measure of relative contributions within the limits of our linear models.Our results indicate that the coefficient for temperature is higher, indicating its greater influence on reproductive outcomes. This finding is consistent with projections suggesting that rising temperatures will have a more pronounced effects on sea urchin fertilization success than changes in pH alone11,33. Previous studies across a range of sea urchin species and thermal environments support the critical role of temperature in determining fertilization success. For example, warming by 2–4 °C enhanced fertilization success in the Antarctic sea urchin Sterechinus neumayeri, particularly at low sperm concentrations34. Similarly, fertilization in tropical species such as Echinometra mathaei and Toxopneustes roseus showed strong sensitivity to warming, with success declining sharply above 30°C35. In Lytechinus variegatus, fertilization was negatively affected by a 3 °C temperature increase, particularly when combined with lowered pH36. Conversely, Heliocidaris erythrogramma, a temperate species, exhibited resilience to warming up to 6 °C, although developmental impacts became evident at higher thresholds9. These results align with our findings that temperature had the strongest effect on fertilization in S. droebachiensis and suggest that sensitivity to warming may reflect species-specific thermal tolerances and local adaptation. The contrasting fertilization responses observed across sea urchin species likely reflect differences in physiological adaptation to their native thermal and chemical environments. Polar species such as S. droebachiensis and Sterechinus neumayeri are adapted to cold, stable conditions and may be more susceptible to rapid warming. In contrast, tropical and temperate species like Heliocidaris erythrogramma often tolerate higher temperatures but can be more sensitive to reduced pH, potentially due to species-specific differences in gamete structure and function. For instance, ocean acidification has been shown to impair sperm motility and reduce fertilization success in several echinoid species34, while also disrupting intracellular pH regulation and acrosomal function26. Conversely, some species show fertilization success to be more strongly influenced by temperature than by pH9. These divergent patterns emphasize the role of environmental history, local adaptation, and physiological plasticity in shaping species-specific vulnerabilities to climate change. Sensitivity analysis further supports this, showing that temperature’s coefficient indicates greater sensitivity compared to pH. This underscores the importance of focusing on temperature management in future conservation and mitigation strategies.Combined stressor effectsWhen comparing our results with other studies investigating the effects of multiple stressors (mostly pH, temperature and salinity, although some studies tested contaminants too) on the fertilization success in sea urchins, we observe contrasting results. Similarly to our observations, fertilization in the tropical sea urchin Pseudoboletia indiana, significantly decreased with acidification (− 0.3 to 0.5 pH units) while increasing with warming33. Other studies have shown no effect of pH or temperature on fertilization success. For example, fertilization success in the tropical sea urchin H. erythrogramma was highly dependent on sperm density, and not on temperature (2–4 °C above ambient), pH (0.4–0.6 pH units below ambient), pCO2 (367–1892 ppm) or their interactions10. Similarly, no significant effect of warming and acidification was found on the percentage of fertilization in a suite of tropical echinoids (H. erythrogramma, Heliocidaris tuberculata, Tripneustes gratilla and Centrostephanus rodgersii) when exposed to all combinations of three temperatures and three pHs10. On the other hand, fertilization in the subtropical sea urchin Heliocidaris crassispina was affected by warming (28 °C to 43 °C) and freshening37. In another species of subtropical rock-boring sea urchin, Echinometra lucunter, fertilization rates were only negatively affected by temperature increase (2 °C), and pH decrease or presence of lead contamination (alone and in combination) did not seem to affect them38. In the same species, pH decrease alone negatively affected the fertilization success, even at optimal temperatures39. These apparent differences between studies may likely be a consequence of local adaptation to the present range of pH variability40,41. In many of the cited studies, ocean acidification scenarios were selected based on IPCC scenarios for open ocean, neglecting the high range of pH variability in coastal areas. As a consequence, some of the tested pH levels fell within the present range, and thus did not represent a true stress condition nor a realistic ocean acidification scenario per se42.The combined effects of temperature, pH, and salinity on fertilization success were explored using a conceptual framework based on standardized slope coefficients from simple linear regressions (Supplementary 3). These ran under reference conditions for each driver, allowed us to compare relative contributions and visualize how individual stressors scale in relation to one another. The results suggest that temperature had the strongest influence on fertilization, followed by pH, while salinity played a minor role. This approach oversimplifies the modelling of the performance curve for these three drivers and only considers its linear part. The observed statistical interactions between temperature and pH support non-linear effects through unimodal performance curves and/or potential stressors interactions. New experiments expanding the pH and temperature range as well as additional information on the drivers’ mode of actions would be needed to fully resolve the combined effect.Lessons-learned applied to Arctic regionsIt is important to note that this study focuses on a single population of S. droebachiensis collected from Kongsfjorden, Svalbard, which may not fully represent the species’ diversity. The green sea urchin has a broad distribution across the North Atlantic and Arctic, including populations in Norway, Greenland, and North America, which experience different environmental conditions. It is therefore possible that other populations may exhibit different sensitivities to temperature and pH based on local adaptation. Additionally, the fertilization assays were conducted using gametes from two females and four males, which limits the genetic diversity represented in the experimental crosses. Fieldwork in polar environments poses substantial logistical and ethical challenges, and the opportunity to collect and experimentally test high numbers of individuals and multiple populations is often constrained by weather, accessibility, and seasonal availability. Despite these limitations, our findings provide critical data points for Arctic populations and can serve as a baseline for future comparative studies aimed at understanding intraspecific variability across broader spatial and environmental gradients.Research in polar sea urchins is rather limited, and the few available studies show contrasting results. Fertilization in the Antarctic sea urchin Sterechinus neumayeri was resilient to acidification (pH ~ 7.7 and 7.5) at ambient temperature (0 °C), but not under elevated temperatures (1.5 °C and 3 °C), where it was negatively impacted by both drivers independently as well as by their interactions43. However, when tested at a range of sperm concentrations, warming (1 °C, 3 °C and 5 °C; 2–4 °C above ambient) enhanced fertilization in the same sea urchin S. neumayeri at the lowest sperm concentrations, whereas decreased pH (pH 8.0, pH 7.8 and pH 7.6; 0.2–0.4 pH units below ambient) did not have any effect34. Another study of the same species found reduced fertilization success under environmental-relevant OA scenarios although the responses observed indicated a high degree of individual variability44. Even less information is available on Arctic species of sea urchins, and the only available study focused on the effects of acidification alone. Fertilization in S. droebachiensis collected from Kongsfjorden (Svalbard) was impaired by acidification, decreasing significantly at extreme levels of acidification (~ 2000 µatm pCO2) but not at the other levels of acidification tested45, results that contrast with our observations. Although some Antarctic sea urchins have shown resilience to OA when tested in isolation, their fertilization success may still be negatively affected under combined stressors such as warming and acidification, suggesting limited overall resilience under realistic environmental conditions.ConclusionsOur investigation provides insights into the complexity of the relation between environmental stressors and the responses in marine organisms. While our observations align with previous research, it also contrasts with other studies. The variability in responses underscores the difficulties in providing definitive conclusions about the impacts of multiple stressors. This variability is not only present among populations but also among individuals, reflecting the complex nature of multiple stressors and their effects on marine organisms. Our results corroborate the notion that multiple stressor investigations are inherently complicated, and the relationships among stressors are often difficult to discern.Despite these complexities, our study emphasizes the critical need for continued research on organismal and ecosystem responses to multiple stressors, particularly in polar ecosystems. Given that environmental changes are occurring at an accelerated pace in polar regions compared to other parts of the world, understanding these mechanisms is crucial. This research is timely and necessary to develop effective management and conservation strategies for polar marine ecosystems, ensuring their resilience in the face of rapid environmental changes.

Data availability

**Supporting data for this study is available at 10.21334/NPOLAR.2024.9A45CB17 (fertilization data; 21) and 10.21334/NPOLAR.2025.9BF189E2 (calculations of the stress index and individual contributions; 24).**.
ReferencesPörtner, H. O. et al. IPCC, 2019: Summary for policymakers. In: Special Report on the Ocean and Cryosphere in a Changing Climate. (2019).Hoegh-Guldberg, O. & Bruno, J. F. The impact of climate change on the world’s marine ecosystems. Science 328, 1523–1528 (2010).Article 
ADS 
CAS 
PubMed 

Google Scholar 
Fabry, V. J., Seibel, B. A., Feely, R. A. & Orr, J. C. Impacts of ocean acidification on marine fauna and ecosystem processes. ICES J. Mar. Sci. 414–432 (2008).Doney, S. C., Fabry, V. J., Feely, R. A. & Kleypas, J. A. Ocean acidification: the other CO2 problem. Annual Rev. Mar. Sci. 1, 169–192 (2009).Article 
ADS 

Google Scholar 
AMAP. AMAP Assessment 2018: Arctic Ocean Acidification. Arctic Monitoring and Assessment Programme (AMAP), Tromsø, Norway. (2018).AMAP. Arctic climate change update 2021: Key trends and impacts. Summary for policy-makers. (2021).Byrne, M. Impact of ocean warming and ocean acidification on marine invertebrate life history stages: vulnerabilities and potential for persistence in a changing ocean. Oceanogr. Mar. Biol. Annu. Rev. 49, 1–42 (2011).
Google Scholar 
Ross, P. M., Parker, L., O’Connor, W. A. & Bailey, E. A. The impact of ocean acidification on reproduction, early development and settlement of marine organisms. Water 3, 1005–1030 (2011).Article 
CAS 

Google Scholar 
Byrne, M. et al. Temperature, but not pH, compromises sea urchin fertilization and early development under near-future climate change scenarios. Proc. Royal Soc. B: Biol. Sci. 276, 1883–1888 (2009).Article 

Google Scholar 
Byrne, M. et al. Fertilization in a suite of coastal marine invertebrates from SE Australia is robust to near-future ocean warming and acidification. Mar. Biol. 157, 2061–2069 (2010).Article 

Google Scholar 
Byrne, M., Soars, N., Selvakumaraswamy, P., Dworjanyn, S. A. & Davis, A. R. Sea urchin fertilization in a warm, acidified and high pCO2 ocean across a range of sperm densities. Mar. Environ. Res. 69, 234–239 (2010).Article 
CAS 
PubMed 

Google Scholar 
Dupont, S., Ortega-Martinez, O. & Thorndyke, M. Impact of near-future ocean acidification on echinoderms. Ecotoxicology 19, 449–462 (2010).Article 
CAS 
PubMed 

Google Scholar 
Kurihara, H. & Shirayama, Y. Effects of increased atmospheric CO2 on sea urchin early development. Mar. Ecol. Prog Ser. 274, 161–169 (2004).Article 
ADS 

Google Scholar 
Kroeker, K. J., Kordas, R. L. & Harley, C. D. G. Embracing interactions in ocean acidification research: confronting multiple stressor scenarios and context dependence. Biol Lett 13, (2017).Voronkov, A., Hop, H. & Gulliksen, B. Diversity of hard-bottom fauna relative to environmental gradients in Kongsfjorden, Svalbard. Polar Res. 32, 1–27 (2013).Article 

Google Scholar 
Hop, H. et al. Scientific diving in Arctic Kelp forest to detect climate-related changes. Fram Forum. 14, 46–51 (2025).
Google Scholar 
Hop, H. et al. Zooplankton in Kongsfjorden (1996–2016) in relation to climate change. In The Ecosystem of Kongsfjorden, Svalbard. Advances in Polar Ecology, vol 2 (eds Hop, H. & Wiencke, C.) 185–219 (Springer, Cham, doi:https://doi.org/10.1007/978-3-319-46425-1_7. (2019).Chapter 

Google Scholar 
Fransson, A. et al. Late winter-to-summer change in ocean acidification state in Kongsfjorden, with implications for calcifying organisms. Polar Biol. 39, 1841–1857 (2016).Article 

Google Scholar 
Tverberg, V. et al. The Kongsfjorden transect: seasonal and Inter-annual variability in hydrography. In The Ecosystem of Kongsfjorden, Svalbard. Advances in Polar Ecology, vol 2 (eds Hop, H. & Wiencke, C.) 49–104 (Springer, Cham, doi:https://doi.org/10.1007/978-3-319-46425-1_3. (2019).Chapter 

Google Scholar 
Levitan, D. R. The relationship between conspecific fertilization success and reproductive isolation among three congeneric sea urchins. Evolution 56, 1599–1689 (2002).Article 
PubMed 

Google Scholar 
Espinel-Velasco, N., Kvernvik, A. C., Hop, H. & Dupont, S. Fertilization Rates in the Green Sea Urchin Stronglyocentrotus Droebachiensis from Kongsfjorden when Exposed To Multiple Environmental Stressors [Dataset] (Norwegian Polar Institute, 2025).Pierrot, D., Epitalon, J. M., Orr, J. C., Lewis, E. & Wallace, D. W. R. MS Excel program developed for CO2 system calculations – version 3.0. (2021).R Core Team. R: A language and environment for statistical computing. (2021).Espinel-Velasco, N., Kvernvik, A. C., Hop, H. & Dupont, S. Calculations of the Stress Index and Relative Contributions of Multiple Environmental Stressors on the Fertilization Rates in the Green Sea Urchin Stronglyocentrotus Droebachiensis from Kongsfjorden [Dataset] (Norwegian Polar Institute, 2025).Leuchtenberger, S. G. et al. The effects of temperature and pH on the reproductive ecology of sand dollars and sea urchins: impacts on sperm swimming and fertilization. PLoS ONE 17, (2022).Christen, R., Schackmann, R. W. & Shapiro, B. M. Metabolism of sea urchin sperm. Interrelationships between intracellular pH, ATPase activity, and mitochondrial respiration. J. Biol. Chem. 258, 5392–5399 (1983).Article 
CAS 
PubMed 

Google Scholar 
Mita, M. & Nakamura, M. Energy metabolism of sea urchin spermatozoa: An approach based on echinoid phylogeny. Zool. Sci. 15 1–10. https://doi.org/10.2108/zsj.15.1 (1998). Kapsenberg, L., Okamoto, D. K., Dutton, J. M. & Hofmann, G. E. Sensitivity of sea urchin fertilization to pH varies across a natural pH mosaic. Ecol. Evol. 7, 1737–1750 (2017).Article 
PubMed 
PubMed Central 

Google Scholar 
Hudson, M. E. & Sewell, M. A. Ocean acidification impacts sperm swimming performance and pHi in the new Zealand sea urchin Evechinus chloroticus. J. Exp. Biol. 225, (2022).Havenhand, J. N., Buttler, F. R., Thorndyke, M. C. & Williamson, J. E. Near future levels of ocean acidification reduce fertilization success in a sea urchin. Curr. Biol. 18, R651–R652 (2008).Article 
CAS 
PubMed 

Google Scholar 
Reuter, K. E., Lotterhos, K. E., Crim, R. N., Thompson, C. A. & Harley, C. D. G. Elevated pCO2 increases sperm limitation and risk of polyspermy in the red sea urchin Strongylocentrotus franciscanus. Glob Chang. Biol. 17, 163–171 (2011).Article 
ADS 

Google Scholar 
Kurihara, H. Effects of CO2-driven ocean acidification on the early developmental stages of invertebrates. Mar. Ecol. Prog Ser. 373, 275–284 (2008).Article 
ADS 
CAS 

Google Scholar 
Foo, S. A., Dworjanyn, S. A., Khatkar, M. S., Poore, A. G. B. & Byrne, M. Increased temperature, but not acidification, enhances fertilization and development in a tropical urchin: potential for adaptation to a tropicalized Eastern Australia. Evol. Appl. (2014).Ho, M. A., Price, C., King, C. K., Virtue, P. & Byrne, M. Effects of ocean warming and acidification on fertilization in the Antarctic echinoid Sterechinus neumayeri across a range of sperm concentrations. Mar. Environ. Res. 90, 136–141 (2013).Article 
CAS 
PubMed 

Google Scholar 
Mejía-Gutiérrez, L. M., Benítez-Villalobos, F. & Díaz-martínez, J. P. Effect of temperature increase on fertilization, embryonic development and larval survival of the sea urchin Toxopneustes roseus in the Mexican South Pacific. J. Therm. Biol. 83, 157–164 (2019).Article 
PubMed 

Google Scholar 
Lenz, B., Fogarty, N. D. & Figueiredo, J. Effects of ocean warming and acidification on fertilization success and early larval development in the green sea urchin Lytechinus variegatus. Mar. Pollut Bull. 141, 70–78 (2019).Article 
CAS 
PubMed 

Google Scholar 
Mak, K. K. Y. & Chan, K. Y. K. Interactive effects of temperature and salinity on early life stages of the sea urchin Heliocidaris Crassispina. Mar. Biol. 165, 1–11 (2018).Article 

Google Scholar 
Sperandio Caetano, L. et al. Impact on fertility rate and embryo larval development due to the association acidification, ocean warming and lead contamination of a sea urchin Echinometra lucunter (Echinodermata: Echinoidea). Bull. Environ. Contam. Toxicol. 106, 923–928 (2021).Article 

Google Scholar 
Pereira, T. M. et al. The success of the fertilization and early larval development of the tropical sea urchin Echinometra lucunter (Echinodermata: Echinoidea) is affected by the pH decrease and temperature increase. Mar. Environ. Res. 161, 105106 (2020).Article 
CAS 
PubMed 

Google Scholar 
Vargas, C. A. et al. Species-specific responses to ocean acidification should account for local adaptation and adaptive plasticity. Nat. Ecol. Evol. 1, 1–7 (2017).Article 
CAS 

Google Scholar 
Vargas, C. A. et al. Upper environmental pCO2 drives sensitivity to ocean acidification in marine invertebrates. Nat. Clim. Chang. 12, 200–207 (2022).Article 
ADS 
CAS 

Google Scholar 
Boyd, P. W. et al. Experimental strategies to assess the biological ramifications of multiple drivers of global ocean change—a review. Glob. Chang. Biol. 24 (2018).Ericson, J. A. et al. Combined effects of two ocean change stressors, warming and acidification, on fertilization and early development of the Antarctic echinoid Sterechinus neumayeri. Polar Biol. 35, 1027–1034 (2012).Article 

Google Scholar 
Sewell, M. A., Millar, R. B., Yu, P. C., Kapsenberg, L. & Hofmann, G. E. Ocean acidification and fertilization in the Antarctic sea urchin Sterechinus neumayeri: the importance of polyspermy. Environ. Sci. Technol. 48, 713–722 (2014).Article 
ADS 
CAS 
PubMed 

Google Scholar 
Bögner, D., Bickmeyer, U. & Köhler, A. CO2-induced fertilization impairment in Strongylocentrotus droebachiensis collected in the Arctic. Helgol. Mar. Res. 68, 341–356 (2014).Article 
ADS 

Google Scholar 
Download referencesAcknowledgementsSpecial thanks to Markus Brand (AWI) and the divers from the AWIPEV station in Ny-Ålesund for the collection of the specimens.FundingOpen access funding provided by University of Gothenburg.Author informationAuthor notesNadjejda Espinel-VelascoPresent address: Department of Marine Science, Tjärnö Marine Laboratory, University of Gothenburg, Gothenburg, SwedenAuthors and AffiliationsNorwegian Polar Institute, Fram Centre, Tromsø, 9296, NorwayNadjejda Espinel-Velasco, Ane Cecilie Kvernvik & Haakon HopDepartment of Biological and Environmental Sciences, University of Gothenburg, Kristineberg, Gothenburg, SwedenSam DupontMarine and Freshwater Research Institute, Fornubudir 4, Hafnarfjörður, 220, IcelandSam DupontAuthorsNadjejda Espinel-VelascoView author publicationsSearch author on:PubMed Google ScholarAne Cecilie KvernvikView author publicationsSearch author on:PubMed Google ScholarHaakon HopView author publicationsSearch author on:PubMed Google ScholarSam DupontView author publicationsSearch author on:PubMed Google ScholarContributionsNEV conceptualized and designed the study, conducted the experiments, analysed the data and led the manuscript writing. ACK contributed to the experimental work and reviewed the manuscript. HH provided critical insights into Arctic ecosystems, contributed to the interpretation of results, and assisted in manuscript revisions. SD contributed to the conceptual framework, provided guidance on experimental methodology and data analysis, and critically revised the manuscript. All authors reviewed and approved the final manuscript.Corresponding authorCorrespondence to
Nadjejda Espinel-Velasco.Ethics declarations

Competing interests
The authors declare no competing interests.

Additional informationPublisher’s noteSpringer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.Supplementary InformationBelow is the link to the electronic supplementary material.Supplementary Material 1Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
Reprints and permissionsAbout this articleCite this articleEspinel-Velasco, N., Kvernvik, A.C., Hop, H. et al. Combined effects of ocean acidification, warming, and salinity on the fertilization success in an Arctic population of sea urchins.
Sci Rep 15, 44090 (2025). https://doi.org/10.1038/s41598-025-27725-zDownload citationReceived: 01 April 2025Accepted: 05 November 2025Published: 18 December 2025Version of record: 18 December 2025DOI: https://doi.org/10.1038/s41598-025-27725-zShare this articleAnyone you share the following link with will be able to read this content:Get shareable linkSorry, a shareable link is not currently available for this article.Copy shareable link to clipboard
Provided by the Springer Nature SharedIt content-sharing initiative
Keywords
Strongylocentrotus droebachiensis
Global change related stressorsAnthropogenic pressuresMultiple stressorsBenthic communitiesEarly-life processesKongsfjordenSvalbard More

AbstractUnderstanding the response of biofilms to climatic conditions and their effects on cultural heritage sites is crucial for developing effective conservation strategies. Previous studies have primarily focused on microbial community responses to environmental factors, and little is known about how climatic conditions influence biofilm-induced deterioration and protection. Here we analyzed genomic data from stone heritage sites in East and South Asia and found that biofilm roles shifted from causing deterioration to offering protection along the transition from temperate to tropical climates. This shift was likely regulated by climate-driven variations in functional genes associated with dissimilatory nitrate reduction (napAB, narGHI, nrfAH, and nirBD) and assimilatory sulfate reduction (cysJI and sir). The expression of genes related to these pathways inhibits the accumulation of soluble salts and biogenic acids, leading to protective effects. This study elucidates the dynamic role of microbes in cultural heritage conservation and lays the foundation for sustainable preservation strategies.

Similar content being viewed by others

Climate-driven succession in marine microbiome biodiversity and biogeochemical function

Article
Open access
25 April 2025

Spatial and functional differentiation of microbial biofilms in a traditional cheese ripening environment

Article
Open access
28 November 2025

Soil microbiome engineering for sustainability in a changing environment

Article

30 October 2023

IntroductionAs of September 2023, the World Heritage List recorded a total of 1154 cultural World Heritage Sites worldwide, including 897 cultural heritage sites (https://whc.unesco.org/en/list/). Over time, immovable cultural heritage exposed to the natural environment conditions experiences deterioration due to physical, chemical, and biological weathering processes1,2,3,4. The Angkor monuments in Cambodia5, the Pharaonic sandstone monuments in Luxor, Egypt6, the statues of the Dazu Rock Carvings in China7, and many other stone cultural heritage sites exhibit such damage. This type of damage not only threatens the integrity of these ancient monuments, but also leads to irreparable losses in inherited culture. Therefore, mitigating the destruction of cultural heritage sites and developing protective strategies have garnered widespread global attention.Biofilms are aggregates of microorganisms and their extracellular polymeric substances (polysaccharides, proteins, and environmental DNA)8. Studies have shown that microorganisms, particularly bacteria, play an important role in the formation of biofilms on the stone surfaces of heritage sites9. However, biofilms have been reported to exert two distinct effects on stone: biodeterioration and bioprotection8. Biodeterioration, ranging from material discoloration to structural damage, can result in a loss of the historical and aesthetic value of cultural heritage sites10,11,12. This process can damage the rock structures through physical penetration, organic and inorganic acid corrosion, and the redox reactions of mineral cations13,14,15. Conversely, bioprotection refers to the positive effects of biofilms on cultural heritage by which biofilms form protective layers on stone surfaces through biomineralization, such as calcium carbonate and calcium oxalate precipitation, which enhances resistance to environmental stress16,17. Acting as physical barriers, biofilms can help protect stone surfaces from weathering, acid rain, and UV radiation18. However, microbial deterioration and protection are influenced by a combination of multidimensional factors, such as the microbial community, climatic conditions, and rock properties. The mechanisms underlying the functional transition of biofilms is not well understood, which limits the design and development of biointervention strategies.A growing body of research suggests that climatic gradients may drive the transition of microbial functions between biodeterioration and bioprotection19,20,21. The structure of microbial communities varies greatly across different heritage sites and is influenced not only by the type of stone but also by environmental factors such as temperature, precipitation, sunlight, and salinity22,23,24. Generally, warm and humid climates provide favorable environmental conditions for the growth of most organisms9,25. Precipitation influences the diversity and abundance of microbial communities26, while an increase in salinity reduces microbial and enzyme activity27. An increase in pollutant concentrations lowers the complexity and stability of microbial ecological networks28. Although several studies have examined the effect of climatic conditions on the structure of microbial communities, cross-climatic investigations on microbial deterioration and protection remain scarce. Consequently, the effects of climate-driven transitions on these microbial functions are largely unexplored, and their underlying mechanisms are yet to be elucidated.This study evaluated the dual roles of biofilms at cultural heritage sites by examining the relationship between the relative bioprotection rate of biofilms and climatic conditions. Specifically, 16S rRNA data from 91 sampling points from 10 World Heritage Sites across different climatic environments in East and South Asia were analyzed. By assessing the microbial community structures under varying climatic conditions, we explored the impact of climate on the dominant species, diversity, and network structure of microbial communities. Additionally, functional genes associated with biofilm protection and metabolic processes involved in deterioration were predicted to determine the biofilm protection rate. This rate was subsequently used to evaluate the impact of biological protection and deterioration on cultural heritage. Finally, the relationship between biofilm protection rate and climatic conditions was analyzed to elucidate the key mechanisms through which climatic gradients affect biofilm function. The results of this study provide new insights into the climate-driven transition of biofilm functions, promoting the sustainable conservation of stone cultural heritage sites by mitigating microbial deterioration and designing innovative biointervention strategies.ResultsClimate-driven variations of microbial communities at stone cultural heritage sitesTo determine how climatic conditions affect microbial communities on stone surfaces, the microbial community compositions of 91 sampling points from ten heritage sites in East and South Asia under different climatic conditions were analyzed. Based on the sampling location, these microbial communities were classified into three groups: N, ST, and T, which corresponded to temperate, subtropical, and tropical climates, respectively. Non-metric multidimensional scaling (NMDS) based on Bray–Curtis distances reflected the β-diversity of microbial community composition across different climatic environments (stress = 0.184) (Fig. 1a). The results showed significant differences in microbial structures among the different climatic conditions (PERMANOVA: F = 4.308, P < 0.001). Notably, the ST group overlapped with both the N and T groups, whereas the N and T groups were distinctly separated. This indicated that bacterial communities different notably between temperate and tropical climates.Fig. 1: Microbial community structure of heritage sites under different climatic environments.a Non-metric multidimensional scaling (NMDS) analysis of microbial communities based on species. b Relationships between mean annual precipitation (MAP) and Shannon index of bacterial communities. c Relationships between MAP and Chao index of bacterial communities. d Relationships between Minimum Temperature of the Coldest Month (MINTCM) and Shannon index of bacterial communities. e Relationships between MINTCM and Chao index of bacterial communities. f The top 10 bacterial phyla in microbial communities at heritage sites across different climatic environments.Full size imageChao and Shannon indices were used to measure the alpha diversity of different bacterial communities. Both indices exhibited an increasing trend with increasing precipitation (Figs. 1b, c). The increase in precipitation can provide sufficient water resources for microbial communities, thereby promoting the growth in the total species number and enhancing species richness within the community, ultimately increasing the overall species diversity29. However, the relationship between bacterial alpha diversity and minimum temperature of the coldest month (MINTCM) was not linear (Figs. 1d, e). The highest Shannon and Chao indices were observed in samples with minimum temperature between 0 and 10 °C (approximately 25–30°N latitude), likely because temperature extremes increase interspecies competition and favor the growth of strains better adapted to the environment. In contrast, the abundance of other species with low-temperature adaptability decreased or even ceased to exist30. These results suggest that the alpha diversity of the microbial community was strongly influenced by precipitation. The ranking of climatic factors based on their importance to the alpha diversity of the microbial community further supports this, as precipitation-related factors, including precipitation of the wettest quarter (PWETQ), precipitation of the warmest quarter (PWARQ), and mean annual precipitation (MAP), received the highest importance scores among all climatic factors (Supplementary Fig. 1).In terms of community composition, Actinomycetota, Pseudomonadota, and Cyanobacteriota were the most widely distributed bacterial phyla on the stone surfaces (Fig. 1f). Their combined relative abundances exceeded 50%, making them the dominant bacterial phyla across all studied heritage sites. Climatic conditions significantly influenced the abundance of the dominant bacterial phyla (Supplementary Fig. 2). The relative abundances of Actinomycetota and Cyanobacteriota decreased with increasing precipitation and temperature, whereas Pseudomonadota accumulated along the precipitation and temperature gradients. Actinomycetota and Cyanobacteriota are typically adapted to drier environments31,32. In particular, the abundance of Cyanobacteriota decreased as water availability decreases32. Conversely, Pseudomonadota are sensitive to water, and their abundance tends to show strong and consistent variation with precipitation33.To reveal potential bacterial interactions at the heritage sites under different climatic conditions, an operational taxonomic units (OTU)-based co-occurrence network was constructed (Fig. 2; Supplementary Table 1). Among these, microbial communities in tropical climates (T) exhibited the most complex network structures. Their core networks, composed of nodes (OTUs) and edges (correlations between OTUs), were the largest, with N = 1024 and L = 10,214. This indicated that bacterial communities in tropical climates have stronger interconnections than those in temperate climates. Studies have shown that in warm and humid environments, microbial communities tend to be more complex because they can sustain stronger ecological interactions34. These stable environmental conditions provide ample water and suitable temperatures for bacterial growth and reproduction, thereby promoting interconnectivity. However, the Venn diagram showed that the subtropical bacterial community had the highest number of unique OTUs (Supplementary Fig. 3), which is likely related to the broader distribution range of the subtropical samples.Fig. 2: Co-occurrence network analysis of microbial communities in heritage sites.The co-occurrence network was based on Spearman’s correlation between OTUs, which were used to construct the network exists in at least 60 % of the samples. All the connections have correlation coefficients r > |0.9| and a P < 0.01. Nodes were colored according to different phylum levels. N denotes the node and L denotes the edge.Full size imageClimate-driven microbial functions of biofilmsBiofilms on stone surfaces may exert both protective and detrimental effects, which are controlled by microbial metabolic functions. To assess the influence of climatic conditions on microbial metabolic function, we analyzed the functional genes related to negative or positive metabolic processes in biofilms (Supplementary Table 2)35. Microbial metabolic functions were closely correlated with climatic factors (Fig. 3a). Correlations between positive and negative metabolic functions and climatic factors were consistent. These functions were significantly positively correlated with temperature or precipitation extremes (such as minimum temperature (MINTCM) and precipitation of the wettest quarter (PWETQ)), and negatively correlated with temperature seasonality (TSEA) (P < 0.001). This relationship may occur because higher temperatures and humid environments typically enhance microbial growth and metabolic activity, thereby promoting metabolic processes36. However, when temperature fluctuations are large, they may affect the stability of microbial metabolism, inhibiting the normal progression of certain metabolic functions37. Notably, the maximum temperature of the warmest month (MAXTWM) was negatively correlated with assimilatory nitrate reduction and denitrification but positively correlated with oxalate biosynthesis. Although these three metabolic processes are considered to have positive functions, they exhibit different trends under the influence of different climatic factors. This suggests that the climatic factors dynamically affect biofilm functions.Fig. 3: Climate-driven microbial functions of biofilms.a Pearson correlation heatmap between climatic factors and microbial metabolic function. *P < 0.05, **P < 0.01, ***P < 0.001. The abbreviation for climatic condictions was detailed at Supplementary Table 3. b The relationship between minimum temperature (MINTCM) and biofilm protection rate. c The relationship between mean annual precipitation (MAP) and biofilm protection rate.Full size imageTo quantitatively assess the influence of bacterial communities on cultural heritage materials, the “biofilm protection ratio” that characterizes the dual role of biofilm was analysed8,38. The “biofilm protection ratio” is the ratio between the sum of natural bioprotection effect genes and the sum of biodeterioration effect genes. A value greater than 1 indicates that the biofilm plays a protective role, meaning that the non-colonized areas exhibit a more severe deterioration pattern than biofilm-colonized areas. According to this study, the bioprotection ratio was closely related to climatic factors (Figs. 3b, c). Specifically, MINTCM and MAP showed a positive correlation with the bioprotection ratio, with correlation coefficients of R = 0.67 and R = 0.44, respectively (Figs. 3b and 3c, P < 0.001). This indicated that the protective role of biofilms is driven by climatic conditions. Biofilms exhibit strong protective effects at heritage sites in warm and humid tropical and subtropical environments.Key functional genes that control the protective role of biofilmsClimate drives the acceleration of the microbial geochemical cycle affecting the deterioration or protection of stone, as evidenced by the significant correlation between the N, S and C cycles and climatic factors (Supplementary Fig. 4). To explore the mechanisms by which microbial communities contribute to the protection of stone cultural heritage sites, we analyzed the variations in functional genes related to nitrogen (N), sulfur (S), and carbon (C) cycles under different climatic conditions (Fig. 4), focusing on their effects on the biofilm protection ratio. These functional genes regulate key metabolic pathways and drive biofilms formation to protect or degrade stone surfaces under different climatic conditions.Fig. 4: The relationship between microbial community metabolic capacity and MINTCM in heritage sites.a The relationship between microbial N cycle-related genes and MINTCM. b The relationship between microbial S cycle-related genes and MINTCM. Among them, the red arrow represents the potential positive metabolism, and the blue arrow represents the potential negative metabolism.Full size imageThe abundance of nitrogen cycle genes was significantly correlated with the minimum annual temperature (MINTCM) (P < 0.001, Fig. 4a). The abundance of dissimilatory nitrate-reduction genes (napAB, narGHI, nrfAH, and nirBD) increased with increasing temperature (R = 0.39, P < 0.001). These genes facilitate the conversion of nitrate into ammonium (NO₃⁻ → NO2⁻ → NH₄⁺)39, which reduced the risk of acid corrosion of stone caused by nitrate accumulation and enhanced the protective effect of biofilm on heritage40. In addition, the abundance of assimilated nitrate reduction genes (nasAB and nirA) and denitrification-related genes (nirKS, norBC, and norZ) decreased with increasing temperature. However, its abundance is much lower than that of the dissimilatory nitrate reduction gene, which makes biofilms in tropical environments more likely to exert protective effects. This finding shows that the effect of biofilms on stone is a dynamic proces dominated by dissimilatory nitrate reduction.Moreover, assimilatory sulfate reduction and dissimilatory sulfate reduction play positive roles in the S cycle of microorganisms, whereas sulfur oxidation plays a negative role. The abundance of genes involved in assimilatory sulfate reduction (cysJI and sir) increased with increasing minimum annual temperatures (R = 0.24, P < 0.05). In warmer environments, microorganisms preferentially utilize sulfate (SO₄²⁻) to synthesize sulfur-containing organic compounds, thus minimizing the corrosive effects of sulfuric acid on stone surfaces. Conversely, genes associated with negative effects, such as sulfur oxidation, were negatively correlated with temperature (Fig. 4b). In hot and humid environments, the detrimental metabolic processes of biofilms were weaker.Chemolithoautotrophic organisms can fix CO₂ through the Calvin–Benson–Bassham (CBB) and Wood–Ljungdahl (WL) pathways41,42,43. Genes associated with these cycles account for 0.25%–0.90% of the entire dataset. The genes associated with Carbon fixation were not significantly correlated with MINTCM (P > 0.05). Microorganisms involved in the C cycle can thrive in both cold and warm environments. On nutrient-poor rock surfaces, carbon fixation plays an important role in the initial formation of biofilms by providing bioavailable organic carbon. However, from the perspective of cultural heritage preservation, this process may have potential negative impacts. The organic matter produced by carbon fixation can create favorable conditions for other microorganisms, such as algae and lichens, which may accelerate the degradation of heritage materials.Furthermore, several metabolic pathways with potential protective functions on the stone surface, including biomineralization, oxalate biosynthesis and metal resistance-related functional gene abundance, were positively correlated with MINTCM (P < 0.01). Biomineralization-related genes may promote the precipitation of carbonate minerals, enhancing the connection of micro-pores in the stone18. The oxalate biosynthetic pathway can form an insoluble calcium oxalate layer, providing a protective film that resists further chemical erosion44. Genes related to metal resistance help mitigate microbial stress caused by heavy metals, thereby stabilizing the biofilm structure and alleviating the deterioration process45. Therefore, in warmer areas, the protective effect of biofilms may be further enhanced.DiscussionClimate is vital for microbial communities and functions of stone heritageThe microbial communities on stone cultural heritage sites under different climatic conditions were predominantly composed of Actinomycetota, Pseudomonadota, and Cyanobacteriota, each playing a distinct ecological role. Actinomycetota and Pseudomonadota, primarily associated with filamentous bacteria46, were observed on sandstone surfaces such as those at Dazu Rock Carvings. Studies have shown that these microbial groups are influenced by rainfall47, which can modify the conditions on stone surfaces and impact microbial colonization. Cyanobacteriota were present in all samples as photosynthetic autotrophs48. Notably, the genus Chroococcidiopsis within Cyanobacteriota is commonly found in extreme and arid habitats and was relatively abundant in samples from Beishiku Grottoes12. Their widespread distribution across various habitat types suggests strong adaptability to arid conditions. Network analysis revealed the potential bacterial interactions at heritage sites under different climatic conditions. Pechlivanis et al. found that microbial communities in tropical regions exhibited lower diversity but maintained a denser network, which is consistent with the conclusion of this study49.In recent years, the beneficial role of biofilms on stone heritage, that is, bioprotection, has garnered increasing attention. Our findings suggest that biofilms in tropical environments may offer greater protection (Fig. 5). This observation could be linked to the high expression levels of genes related to positive metabolic pathways in microbial communities under tropical climate. In particular, genes associated with nitrate reduction (napAB and narGHI) and sulfate reduction (cysJI and sir) contribute to mitigating material degradation by consuming soluble salts35. The dissolution and crystallization of salt can exert physical stress to stone, resulting in erosion, cracking, or pulverization50. Nitrate reduction converts nitrate to nitrite or ammonium (NH4+), while sulfate reduction converts sulfate to hydrogen sulfide. Both processes help to alleviate the damage to stone by altering the salt concentration in the surrounding environment.Fig. 5: Protection and deterioration of cultural heritage by biofilms in different climate environments.Compared to temperate climates, microorganisms in tropical climates exhibit more complex network structures and higher diversity. The strong expression of protection-related genes makes biofilms in tropical climates more inclined to have a protective effect on cultural heritage.Full size imageIt should be emphasized that, although this study indicates that these genes may be related to heritage bioprotection or biodeterioration, the underlying mechanisms still needs further study and verification. For example, while assimilatory sulfate reduction primarily supports biosynthetic processes51, the hydrogen sulfide (H₂S) produced during dissimilatory sulfate reduction—a volatile and highly reactive gas— has the potential to exacerbate the deterioration of stone surfaces52. When specific metal ions, such as iron, are present, H₂S can be oxidized to sulfuric acid (H₂SO₄), which then reacts with the calcium in the rock to form gypsum (CaSO₄·2H₂O). This leads to the accumulation of acidic substances, ultimately resulting in surface crust formation, granular disintegration, increased porosity, and loss of cohesion53. Therefore, although sulfate reduction consumes sulfate ions, its products, especially H₂S produced during dissimilatory sulfate reduction, may still react with minerals and accumulate in poorly ventilated or humid environments, thereby triggering localized stone deterioration.Oxalate biosynthesis and biomineralization are the crucial processes to protect the stone heritage from erosion35. Genes associated with these processes were identified in each sample. The biosynthetic pathway of oxalate is particularly important in the formation of insoluble calcium oxalate crystals. These crystals can form a protective barrier on the surface of the stone and effectively resist the chemical erosion of environmental factors such as acid rain and ultraviolet light54. In addition, genes related to biomineralization can promote the precipitation of carbonate calcium, which can fill pores to strengthen the microstructure of stone55. Gene markers related to metal ion absorption and regulation can be used as an indicator of ion extraction in rocks56,57. Metabolism genes of siderophores indicate that bacteria promote the absorption of iron released during rock dissolution by generating siderophores, thereby enhancing the leaching of mineral elements56. Metal resistance genes also help reduce heavy metal stress, stabilize biofilm structure, and slow down biofilm deterioration45,57.Sustainable conservation strategies for cultural heritageOur study reveals the dual role of biofilms in the protection and deterioration of cultural heritage sites, with climate being a key factor influencing biofilm function. Specifically, biofilms in tropical climates may offer protective effects, but they can also contribute to deterioration through biochemical processes (such as sulfate and nitrate reduction) and physical processes, including the retention of capillary water, which accelerates abiotic weathering. In contrast, biofilms in cold and dry temperate climates tend to exhibit stronger deteriorative effects. Therefore, biointervention strategies should be tailored to the specific climatic conditions of heritage sites. In hot and humid areas, while biofilms may have protective benefits, it is essential to monitor and maintain them in their original state to balance their protective and deteriorative effects. In these regions, microbial community-directed regulators, such as the inhibition of acid-producing bacteria or the inoculation of bacteria with nitrate-reducing functions, could be developed using environmental DNA detection. Additionally, the use of nanocoating materials with strong antibacterial properties and good compatibility may offer effective protection for stone heritage sites in these climates. Therefore, dynamically managing biofilms based on environmental conditions and adopting customized protection strategies are considered critical for future cultural heritage protection.It is important to evaluate the metabolic processes that may affect the stone heritage to assess the impact of bacterial communities on stone deterioration. According to the calculation results of biofilm protection rate, it was found that the bacterial community may tend to protect in the tropical environment. However, even so, the bacterial community that inhabits the surface of the rock matrix retains the possibility of a negative impact on the stone. This study attempts to quantify the dual role of biofilms by predicting the metabolic capacity of bacterial communities. This method enables us to determine whether biofilms are more inclined to protect or deteriorate. However, due to the complexity of the natural environment and microbial communities, specific bioprotection or biodeterioration effects need to be further studied.ConclusionThis study investigates the impacts of climatic environments on microbial communities in stone cultural heritage, with a focus on the structural and functional differences across climatic zones. The composition of bacterial communities in temperate, subtropical and tropical regions was analyzed using 16S rRNA data from stone surfaces of East Asian and South Asian cultural heritage sites. The findings demonstrate variations in microbial diversity driven by climatic factors such as temperature and precipitation.Additionally, this study predicted functional genes associated with biofilm protective and degradative metabolic processes, and employed the relative biofilm protective rate to assess the relationship between biodeterioration and bioprotection at cultural heritage sites under different climatic environments. As the climate shifts from temperate to tropical zones, the role of biofilms shifts towards the protective function for cultural heritage. This transition may be modulated by climate-driven variations in functional genes linked to assimilatory nitrate reduction (napAB, narGHI, nrfAH, and nirBD) and assimilatory sulfate reduction (cysJI and sir). The expression of these pathway-related genes inhibits the accumulation of soluble salts and bioacids, thereby generating protective effects.This study provides an exploration to assess the dual role of stone heritage biofilms. Although DNA-based studies mainly reveal potential rather than actual functional capabilities, our data provide valuable insights into the metabolism of bacterial communities on stone heritage. It forms the foundation for the development of sustainable preservation strategies that account for the interaction between microbial activity and climatic conditions.MethodsBioinformatic data of cultural heritageIn this study, high-throughput sequencing data of 91 samples from 10 stone cultural heritage sites in Asia (mainly East Asia and South Asia) were compiled and analyzed (Fig. 6). Each sampling point is from an approximately open environment. Each of heritage site has a unique climate: Beishiku Grottoes (Sample: BSK1-9), Tiantishan Grottoes (Sample: TT1-3), and Maijishan Grottoes (Sample: MJ1-5) in China have temperate climates; Dazu Rock Carvings (Sample: DZ1-14), Leshan Giant Buddha (Sample: LS1-10), Feilai Peak in Hangzhou (Sample: XH1-20), and Leizhou Stone Dog sculptures (Sample: LZ1-5) in China have subtropical climates; the Preah Vihear Temple (Sample: BWX1-8), the Royal Palace of Angkor Thom in Cambodia (Sample: WG1-7), and the Historic Stone Ruins of Tamil Nadu (Sample: India1-10) in India have tropical climates. These sampling points belonging to the ‘Grottoes’ are close to the outdoor environment, so their environmental conditions can be considered similar to those of the open-air sites. Among them, the samples from Dazu were collected and tested by our team, while the data from other heritage sites were from NCBI database. Detailed sample information is shown in Supplementary Table 4.Fig. 6: Global distribution maps of the cultural heritage sampling sites considered in this study.These examples belong to temperate zone (N), subtropical (ST) and tropical (T).Full size imageWater and temperature play an important role in the colonization and deterioration of microbial communities25,26. Appropriate levels of water and temperature contribute to the growth and metabolism of the microbial community. Therefore, 21 bioclimatic variables related to temperature and precipitation, with a resolution of 2.5 km, were obtained from the WorldClim database (www.worldclim.org), following references from other literature24. Temperature-related variables included the mean annual temperature (MAT), temperature seasonality (TSEA), mean diurnal temperature range (MDTR), isothermality (ISO), maximum temperature of the warmest month (MAXTWM), minimum temperature of the coldest month (MINTCM), annual temperature range (TRANGE), mean temperature of the wettest quarter (TWETQ), mean temperature of the driest quarter (TDQ), mean temperature of the warmest quarter (TWARQ), and mean temperature of the coldest quarter (TCQ). The precipitation-related variables include mean annual precipitation (MAP), precipitation of the wettest month (PWETM), precipitation of the driest month (PDM), precipitation seasonality (PSEA), precipitation of the wettest quarter (PWETQ), precipitation of the driest quarter (PDQ), precipitation of the warmest quarter (PWARQ), and precipitation of the coldest quarter (PCQ). The human influence index (HII) and altitude (ALT) were also included. Detailed environmental information is provided in Supplementary Table 5.The 14 samples from Dazu were taken from the statues of Dazu Rock Carvings in Chongqing, China. Biofilm samples containing a small amount of weathered sandstone were carefully collected from the biofilm-sandstone interface using a sterile surgical knife. Each sample weighed approximately 2 g and was stored in a sterile centrifuge tubes. Following the manufacturer’s instructions, DNA was extracted from all samples using the YH-Soil FastPure soil DNA extraction kit (product number: T09-96; mJYH Biotech, Shanghai, China). The quantity and quality of the extracted DNA was assessed using a Nanodrop instrument, and the extraction quality was verified through 1.2 % agarose gel electrophoresis. PCR amplification was conducted using universal primers: 338 F (5′-ACTCCTACGGGAGGCAG-3′) and 806 R (5′-GGACTACHVGGGTWTCTAAT-3′) to target the V3–V4 region of the 16S rRNA gene58. The PCR product was extracted from a 2% agarose gel and purified using the gel extraction kit (AXYGEN Co., China), following the manufacturer’s instructions. The purified amplicons were quantified using Qubit 4.0 (Thermo Fisher Scientific, USA). Amplicons were then pooled in equimolar amounts and paired-end sequenced on an Illumina MiSeq PE300 platform (Illumina, USA), according to standard protocols provided by Majorbio Bio-Pharm Technology Co. Ltd. (Shanghai, China).High-throughput amplicon sequencing analysisThe raw amplicon sequencing dataset was primarily analyzed using the open-source microbiome ecological quantitative insight tool, QIIME. Fastp software was used to remove primer sequences, after which Trimmomatic was used to filter out low-quality sequences with length <150 bases or a 20-base wide moving window average quality score <2059. The original paired-ended sequencing data from each study were first combined using FLASH60. Subsequently, the USEARCH software (v11)61 was applied with the ‘fastx_uniques’ and ‘unoise3’ functions (default parameters) for dereplication and denoising (error correction) of the sequence.Since the data set contains sequencing sequences from different regions of the 16S rRNA gene, analysis at the single nucleotide difference level, such as zero-radius operational taxonomic units (zOTU), was not possible and the sequences after quality control cannot be directly compared. Therefore, these fragments were aligned to the full-length sequence of the 16S rRNA gene (SILVA 138 database) using a closed reference process, with the ‘closed_ref ‘function to set a 97 % identity threshold62. Sequences that did not match any entries in the database were assigned to the “unclassified” or “norank” categories. The matched full-length sequences and their annotations were used as representative sequences and taxonomic classifications for subsequent studies. However, the use of different primers for different taxa introduces an objective bias in detection efficiency, which may affect the detection results of diversity, relative abundance, or specific groups.”Statistical analysisAlpha diversity analysis is used to investigate the diversity of species in a localized, homogeneous habitat and is an important component of microbiome diversity studies. This approach included measures of species richness (Chao and Ace) and diversity (Shannon and Simpson). In this study, one-way analysis of variance was used to compare significant differences in microbial diversity metrics between different groups. NMDS is an indirect gradient analysis method based on dissimilarity or distance matrices and is commonly used in microbiome research to display community beta diversity. This was performed using the “Vega” package in R. Spearman correlation between species and environmental factors was calculated using the “psych” package in R. The importance ranking of environmental factors was conducted using the random forest algorithm, with analysis completed in Python. Spearman correlation coefficients between OTUs were calculated using the “psych” package in R to obtain a correlation matrix. OTUs with correlation coefficients (r) greater than |0.9| and P-values less than 0.01 were retained, and a co-occurrence network were generated. The network was imported into Gephi for visualization using the Fruchterman-Reingold layout, and the nodes were colored according to the modules and phyla. The topological parameters of the network, including average node degree, clustering coefficient, average path length, and modularity, were calculated.Microbial community functional predictionThe metagenomic function was predicted using PICRUSt2 (Phylogenetic Investigation of Communities by Reconstruction of Unobserved States) based on OTU representative sequences63. PICRUSt2 is a software suites that includes a series of tools, including HMMER, which aligns OTU representative sequences with reference sequences. While the microbiome functions predicted by the PICRUSt2 should be considered as potential, it can still provide important insights for research64. Functional genes associated with carbon fixation, and nitrogen and sulfur metabolism were identified using the Kyoto Encyclopedia of Genes and Genomes (KEGG) database. Based on the current literature18,65,66, functional genes related to rock interaction were selected, including those associated with metal uptake/ resistance, oxalate biosynthesis, biomineralization, inorganic P dissolution, and organic P mineralization. The metabolic processes potentially involved in the negative and positive effects on the stone heritage are primarily summarized in Table 1. A complete list of the surveyed genes is provided in Supplementary Table 2.Table 1 Metabolic processes that may have positive and negative effects on stone heritageFull size tableThe “biofilm relative bioprotective ratio” is defined as the ratio of the sum of the biofilm-protective effects to the sum of the biodeterioration effects. A value greater than 1 indicated that the biofilm had a protective effect. The following formula was used for the calculations:$$r=frac{{sum }_{i=1}^{n}{x}_{i}}{{sum }_{j=1}^{m}{y}_{j}}$$
(1)
where ({sum }_{i=1}^{n}{x}_{i}) represents the sum of the relative abundances of genes with positive effects in the biofilm, and ({sum }_{j=1}^{m}{y}_{j}) represents the sum of the relative abundances of genes with negative effects in the biofilm; xi and yi denote the relative abundances of individual genes.Reporting summaryFurther information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

Data availability

The original sequence data can be obtained in NCBI, and the corresponding references of the BioProject number is shown in Supplementary Table 4. Meteorological data used in this study are from the WorldClim database (www.worldclim.org) and can be found in Supplementary Table 5. Supplementary Table 2, Supplementary Table 4, Supplementary Table 5 and other data are available at https://doi.org/10.6084/m9.figshare.28738748.
ReferencesYang, H. Q., Chen, C. W., Ni, J. H. & Karekal, S. A hyperspectral evaluation approach for quantifying salt-induced weathering of sandstone. Sci. Total Environ. 885, 163386 (2023).Article 

Google Scholar 
Li, X., Yang, H., Chen, C., Zhao, G. & Ni, J. Deterioration identification of stone cultural heritage based on hyperspectral image texture features. J. Cult. Herit. 69, 57–66 (2024).Article 

Google Scholar 
Peng, L. X., Bo, W., Yang, H. Q. & Li, X. Y. Deep learning-based image compression for enhanced hyperspectral processing in the protection of stone cultural relics. Expert Syst. Appl. 271, 126691 (2025).Article 

Google Scholar 
Yang, H., Li, X. & Cappitelli, F. Interplay between geological materials and the environment at the Dazu Rock Carvings, China. J. Cult. Herit. 75, 74–83 (2025).Article 

Google Scholar 
Zhang, X., Ge, Q., Zhu, Z., Deng, Y. & Gu, J.-D. Microbiological community of the Royal Palace in Angkor Thom and Beng Mealea of Cambodia by Illumina sequencing based on 16S rRNA gene. Int. Biodeter. Biodegrad. 134, 127–135 (2018).Article 
CAS 

Google Scholar 
Fitzner, B., Heinrichs, K. & Bouchardiere, D. L. Weathering damage on Pharaonic sandstone monuments in Luxor-Egypt. Build. Environ. 38, 1089–1103 (2003).Article 

Google Scholar 
Li, X. et al. A Non-Contact Disease Identification Method for Geoheritage Based on the Spectra and Image Fusion Technology. Geoheritage 17, 100 (2025).Article 

Google Scholar 
Berti, L., Villa, F., Toniolo, L., Cappitelli, F. & Goidanich, S. Methodological challenges for the investigation of the dual role of biofilms on outdoor heritage. Sci. Total Environ. 954, 176450 (2024).Article 
CAS 

Google Scholar 
Zhu, C. et al. Analysis of the Microbiomes on Two Cultural Heritage Sites. Geomicrobiol. J. 40, 203–212 (2023).Article 
CAS 

Google Scholar 
Silva, I. et al. Microbial i nduced stone discoloration in alcobaça monastery: A comprehensive study. J. Cult. Herit. 67, 248–257 (2024).Article 

Google Scholar 
Santo, A. P. et al. Black on White: Microbial Growth Darkens the External Marble of Florence Cathedral. Appl. Sci. 11, 6163 (2021).Article 
CAS 

Google Scholar 
Wu, F. S. et al. Community structures of bacteria and archaea associated with the biodeterioration of sandstone sculptures at the Beishiku Temple. Int. Biodeter. Biodegrad. 164, 105290 (2021).Article 
CAS 

Google Scholar 
Ding, X. et al. Microbiome and nitrate removal processes by microorganisms on the ancient Preah Vihear temple of Cambodia revealed by metagenomics and N-15 isotope analyses. Appl. Microbiol. Biotechnol. 104, 9823–9837 (2020).Article 
CAS 

Google Scholar 
Li, J. et al. The active microbes and biochemical processes contributing to deterioration of Angkor sandstone monuments under the tropical climate in Cambodia – A review. J. Cult. Herit. 47, 218–226 (2021).Article 

Google Scholar 
Yang, H., Cappitelli, F. & Li, X. Pollution gradients shape structure and functions of stone heritage bacterial communities at global scale. Sci. Total Environ. 971, 179087 (2025).Article 
CAS 

Google Scholar 
Tian, J., Ge, F., Zhang, D. Y., Deng, S. Q. & Liu, X. W. Roles of Phosphate Solubilizing Microorganisms from Managing Soil Phosphorus Deficiency to Mediating Biogeochemical P Cycle. Biology 10, 158 (2021).Article 
CAS 

Google Scholar 
Dreyfuss, T. Interactions on site between powdering porous limestone, natural salt mixtures and applied ammonium oxalate. Herit. Sci. 7, 5 (2019).Article 

Google Scholar 
Liu, X. B., Koestler, R. J., Warscheid, T., Katayama, Y. & Gu, J. D. Microbial deterioration and sustainable conservation of stone monuments and buildings. Nat. Sustain. 3, 991–1004 (2020).Article 

Google Scholar 
Lanzén, A. et al. Multi-targeted metagenetic analysis of the influence of climate and environmental parameters on soil microbial communities along an elevational gradient. Sci. Rep. 6, 28257 (2016).Article 

Google Scholar 
Bui, A. et al. Soil fungal community composition and functional similarity shift across distinct climatic conditions. FEMS Microbiol. Ecol. 96, fiaa193 (2020).Article 
CAS 

Google Scholar 
Donhauser, J. & Frey, B. Alpine soil microbial ecology in a changing world. FEMS Microbiol. Ecol. 94, fiy099 (2018).Article 
CAS 

Google Scholar 
Viles, H. A. & Cutler, N. A. Global environmental change and the biology of heritage structures. Global Change Biol. 18, 2406–2418 (2012).Article 

Google Scholar 
Qian, Y., Liu, X., Hu, P., Gao, L. & Gu, J.-D. Identifying the major metabolic potentials of microbial-driven carbon, nitrogen and sulfur cycling on stone cultural heritage worldwide. Sci. Total Environ. 954, 176757 (2024).Article 
CAS 

Google Scholar 
Yu, Y. J. et al. Unearthing the global patterns of cultural heritage microbiome for conservation. Int. Biodeter. Biodegrad. 190, 105784 (2024).Article 
CAS 

Google Scholar 
Yuan, M. M. et al. Climate warming enhances microbial network complexity and stability. Nature Clim. Change 11, 343–U100 (2021).Article 

Google Scholar 
Liu, X., Meng, H., Wang, Y., Katayama, Y. & Gu, J.-D. Water is a critical factor in evaluating and assessing microbial colonization and destruction of Angkor sandstone monuments. Int. Biodeter. Biodegrad. 133, 9–16 (2018).Article 
CAS 

Google Scholar 
Xiao, W., Chen, X., Jing, X. & Zhu, B. A. A meta-analysis of soil extracellular enzyme activities in response to global change. Soil Biol. Biochem. 123, 21–32 (2018).Article 
CAS 

Google Scholar 
Zhang, H. et al. Pollution gradients shape the co-occurrence networks and interactions of sedimentary bacterial communities in Taihu Lake, a shallow eutrophic lake. J. Environ. Manage. 305, 114380 (2022).Article 
CAS 

Google Scholar 
He, Y., Yu, M., Ding, G., Wang, C. & Zhang, F. Precipitation amount and event intervals interact to change plant diversity during dry years in a desert shrubland. Ecol. Indic. 145, 109701 (2022).Article 
CAS 

Google Scholar 
Shen, Z. et al. Warming reduces bacterial diversity and stability in Lake Bosten. J. Environ. Manage. 375, 124352 (2025).Article 

Google Scholar 
Wang, Y. et al. Changes in bacterial community composition and soil properties altered the response of soil respiration to rain addition in desert biological soil crusts. Geoderma 409, 115635 (2022).Article 
CAS 

Google Scholar 
Jung, P. et al. Water availability shapes edaphic and lithic cyanobacterial communities in the Atacama Desert. J. Phycol. 55, 1306–1318 (2019).Article 
CAS 

Google Scholar 
Li, X., Yan, Y., Lu, X., Fu, L. & Liu, Y. Responses of soil bacterial communities to precipitation change in the semi-arid alpine grassland of Northern Tibet. Front. Plant Sci. 13, 1036369 (2022).Article 

Google Scholar 
Yuan, M. M. et al. Climate warming enhances microbial network complexity and stability. Nature Clim. Change 11, 343–348 (2021).Article 

Google Scholar 
Mugnai, G. et al. Ecological strategies of bacterial communities in prehistoric stone wall paintings across weathering gradients: A case study from the Borana zone in southern Ethiopia. Sci. Total Environ. 907, 168026 (2024).Article 
CAS 

Google Scholar 
Hu, Y. L. et al. Effect of increasing precipitation and warming on microbial community in Tibetan alpine steppe. Environ Res 189, 109917 (2020).Article 
CAS 

Google Scholar 
Zhang, Y. et al. Temperature fluctuation promotes the thermal adaptation of soil microbial respiration. Nat. Ecol. Evol. 7, 205–213 (2023).Article 

Google Scholar 
Liu, X. et al. Biofilms on stone monuments: biodeterioration or bioprotection? Trends Microbiol 30, 816–819 (2022).Article 
CAS 

Google Scholar 
Chen, X. et al. Spatiotemporal successions of N, S, C, Fe, and As cycling genes in groundwater of a wetland ecosystem: Enhanced heterogeneity in wet season. Water Res. 251, 121105 (2024).Article 
CAS 

Google Scholar 
Ding, X. et al. An internal recycling mechanism between ammonia/ammonium and nitrate driven by ammonia-oxidizing archaea and bacteria (AOA, AOB, and Comammox) and DNRA on Angkor sandstone monuments. Int. Biodeter. Biodegrad. 165, 105328 (2021).Article 
CAS 

Google Scholar 
Greening, C. & Grinter, R. Microbial oxidation of atmospheric trace gases. Nat. Rev. Microbiol. 20, 513–528 (2022).Article 
CAS 

Google Scholar 
Imhoff, J. F., Rahn, T., Künzel, S. & Neulinger, S. C. Phylogeny of Anoxygenic Photosynthesis Based on Sequences of Photosynthetic Reaction Center Proteins and a Key Enzyme in Bacteriochlorophyll Biosynthesis, the Chlorophyllide Reductase. Microorganisms 7, 576 (2019).Article 

Google Scholar 
Ray, A. E. et al. Atmospheric chemosynthesis is phylogenetically and geographically widespread and contributes significantly to carbon fixation throughout cold deserts. ISME J 16, 2547–2560 (2022).Article 
CAS 

Google Scholar 
Karche, T. & Singh, M. R. Biologically induced calcium oxalate mineralization on 15th century lime mortar, Murud Sea fort, India. J. Archaeol. Sci. Rep. 39, 103178 (2021).
Google Scholar 
Ren, X. K. et al. Engineered microbial platform confers resistance against heavy metals via phosphomelanin biosynthesis. Nat. Commun. 16, 4836 (2025).Article 
CAS 

Google Scholar 
López-Archilla, A. I., Gérard, E., Moreira, D. & López-García, P. Macrofilamentous microbial communities in the metal-rich and acidic River Tinto, Spain. FEMS Microbiol. Lett. 235, 221–228 (2004).Article 

Google Scholar 
Lee, E. S. et al. Distribution and characteristics of geosmin and 2-MIB-producing actinobacteria in the Han River, Korea. Water Supply 20, 1975–1987 (2020).Article 
CAS 

Google Scholar 
Filippidou, S. et al. A Combination of Extreme Environmental Conditions Favor the Prevalence of Endospore-Forming Firmicutes. Front. Microbiol. 7, 1707 (2016).Pechlivanis, N. et al. Microbial co-occurrence network demonstrates spatial and climatic trends for global soil diversity. Sci. Data 11, 672 (2024).Article 

Google Scholar 
Sun, B. et al. Experimental study on the effects of hydrochemistry and periodic changes in temperature and humidity on sandstone weathering in the Longshan Grottoes. Herit. Sci. 11, 173 (2023).Article 
CAS 

Google Scholar 
Zhou, L. J. et al. Assimilatory and dissimilatory sulfate reduction in the bacterial diversity of biofoulant from a full-scale biofilm-membrane bioreactor for textile wastewater treatment. Sci. Total Environ. 772, 145414 (2021).Article 

Google Scholar 
Huber, B., Drewes, J. E., Lin, K. C., König, R. & Müller, E. Revealing biogenic sulfuric acid corrosion in sludge digesters: detection of sulfur-oxidizing bacteria within full-scale digesters. Water Sci. Technol. 70, 1405–1411 (2014).Article 
CAS 

Google Scholar 
Yang, H., Ni, J., Chen, C. & Chen, Y. Weathering assessment approach for building sandstone using hyperspectral imaging technique. Herit. Sci. 11, 70 (2023).Article 

Google Scholar 
Doherty, B. et al. Durability of the artificial calcium oxalate protective on two Florentine monuments. J. Cult. Herit. 8, 186–192 (2007).Article 

Google Scholar 
Daskalakis, M. I. et al. Pseudomonas, Pantoea and Cupriavidus isolates induce calcium carbonate precipitation for biorestoration of ornamental stone. J. Appl. Microbiol. 115, 409–423 (2013).Article 
CAS 

Google Scholar 
Nir, I., Barak, H., Kramarsky-Winter, E., Kushmaro, A. & de Los Ríos, A. Microscopic and biomolecular complementary approaches to characterize bioweathering processes at petroglyph sites from the Negev Desert, Israel. Environ. Microb. 24, 967–980 (2022).Article 

Google Scholar 
Potysz, A. & Bartz, W. Bioweathering of minerals and dissolution assessment by experimental simulations-Implications for sandstone rocks: A review. Const. Build. Mater. 316, 125862 (2022).Article 
CAS 

Google Scholar 
Zeng, G. M. et al. Response of compost maturity and microbial community composition to pentachlorophenol (PCP)-contaminated soil during composting. Bioresour. Technol. 102, 5905–5911 (2011).Article 
CAS 

Google Scholar 
Bolger, A. M., Lohse, M. & Usadel, B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics 30, 2114–2120 (2014).Article 
CAS 

Google Scholar 
Magoc, T. & Salzberg, S. L. FLASH: fast length adjustment of short reads to improve genome assemblies. Bioinformatics 27, 2957–2963 (2011).Article 
CAS 

Google Scholar 
Edgar, R. C. UPARSE: highly accurate OTU sequences from microbial amplicon reads. Nat. Methods 10, 996 (2013).Article 
CAS 

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

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

Google Scholar 
Chu, Y., Zhang, X., Tang, X., Jiang, L. & He, R. Uncovering anaerobic oxidation of methane and active microorganisms in landfills by using stable isotope probing. Environ. Res. 271, 121139 (2025).Article 
CAS 

Google Scholar 
Liang, J. L. et al. Novel phosphate-solubilizing bacteria enhance soil phosphorus cycling following ecological restoration of land degraded by mining. ISME J 14, 1600–1613 (2020).Article 
CAS 

Google Scholar 
Luo, G. et al. Soil Carbon, Nitrogen, and Phosphorus Cycling Microbial Populations and Their Resistance to Global Change Depend on Soil C:N:P Stoichiometry. mSystems 5, e00162–00120 (2020).Article 
CAS 

Google Scholar 
Zerboni, A. Holocene rock varnish on the Messak plateau (Libyan Sahara): Chronology of weathering processes. Geomorphology 102, 640–651 (2008).Article 

Google Scholar 
Wild, B., Gerrits, R. & Bonneville, S. The contribution of living organisms to rock weathering in the critical zone. npj Mater. Degrad. 6, 98 (2022).Article 

Google Scholar 
Ding, X. H. et al. Microbiome characteristics and the key biochemical reactions identified on stone world cultural heritage under different climate conditions. J. Environ. Manage. 302, 114041 (2022).Article 
CAS 

Google Scholar 
Zeng, Q., Hao, T. W., Mackey, H. R., van Loosdrecht, M. C. M. & Chen, G. H. Recent advances in dissimilatory sulfate reduction: From metabolic study to application. Water Res 150, 162–181 (2019).Article 

Google Scholar 
Michael, G. et al. Oxalate production by fungi: significance in geomycology, biodeterioration and bioremediation. Fungal Biol. Rev. 28, 36–55 (2014).Article 

Google Scholar 
Burgos-Cara, A., Ruiz-Agudo, E. & Rodriguez-Navarro, C. Effectiveness of oxalic acid treatments for the protection of marble surfaces. Mater. Des. 115, 82–92 (2017).Article 
CAS 

Google Scholar 
Ivanov, A. I., Gavriushkin, A. V., Siunova, T. V., Khasanova, L. A. & Khasanova, Z. M. Resistance of certain strains of Pseudomonas bacteria to toxic effect of heavy metal ions. Mikrobiologiia 68, 366–374 (1999).
Google Scholar 
Download referencesAcknowledgementsThis work was supported by National Natural Science Foundation of China (No. 52179096, 22376221), Natural Science Foundation of Hunan Province, China (No. 2024JJ2074), and Young Elite Scientists Sponsorship Program by CAST (No. 2023QNRC001). This work was partly supported by the High Performance Computing Center of Central South University.Author informationAuthors and AffiliationsSchool of Civil Engineering, Chongqing University, Chongqing, ChinaHaiqing Yang & Xingyue LiSchool of Metallurgy and Environment, Central South University, Changsha, ChinaLiyuan ChaiSchool of Civil Engineering, Tianjin University, Tianjin, ChinaLe WangSchool of Resources and Safety Engineering, Central South University, Changsha, ChinaChongchong QiAuthorsHaiqing YangView author publicationsSearch author on:PubMed Google ScholarXingyue LiView author publicationsSearch author on:PubMed Google ScholarLiyuan ChaiView author publicationsSearch author on:PubMed Google ScholarLe WangView author publicationsSearch author on:PubMed Google ScholarChongchong QiView author publicationsSearch author on:PubMed Google ScholarContributionsC.C.Q. conceived and directed the project; H.Q.Y. and X.Y.L. performed the simulations; X.Y.L., L.Y.C. and L.W. analyzed the experimental results, H.Q.Y., X.Y.L. and C.C.Q. and wrote the manuscript. All authors edited the manuscript before submission.Corresponding authorCorrespondence to
Chongchong Qi.Ethics declarations

Competing interests
The authors declare no competing interests.

Peer review

Peer review information
Communications Earth & Environment thanks Fadwa Jroundi and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Primary Handling Editors: Somaparna Ghosh. A peer review file is available.

Additional informationPublisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.Supplementary informationTransparent Peer Review fileSupplementary InformationReporting SummaryRights and permissions
Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.
Reprints and permissionsAbout this articleCite this articleYang, H., Li, X., Chai, L. et al. Climate-driven transition in microbial deterioration and protection of stone surfaces at cultural heritage sites.
Commun Earth Environ 6, 1019 (2025). https://doi.org/10.1038/s43247-025-02993-9Download citationReceived: 17 June 2025Accepted: 31 October 2025Published: 18 December 2025Version of record: 18 December 2025DOI: https://doi.org/10.1038/s43247-025-02993-9Share this articleAnyone you share the following link with will be able to read this content:Get shareable linkSorry, a shareable link is not currently available for this article.Copy shareable link to clipboard
Provided by the Springer Nature SharedIt content-sharing initiative More

AbstractMicro-watershed prioritization using morphometric parameter analysis is a systematic approach to identifying and ranking micro-watersheds based on their susceptibility to soil erosion. This helps in implementing effective soil and water conservation measures. Soil erosion is one of the earth’s surface environmental problems. The objectives of this study were to (1) analyze 27 morphometric parameters in kiltie watershed. (2) Prioritize micro watersheds using weighted sum priority techniques. A morphometric analysis of 10 micro-watersheds was conducted, assessing characteristics such as linear, areal, shape, and relief, features to estimate erosion vulnerability. Shuttle Radar Topography Mission (SRTM) digital elevation model (DEM) with 30*30 m spatial resolution was used to delineate the micro-watersheds and drainage networks through ArcGIS 10.7.1 software. Micro watersheds were ranked based on their soil erosion susceptibility by using a weighted sum priority (WSP) value derived from various morphometric parameters. The morphometric analysis and erosion evaluation results showed that six micro watersheds (MW1, MW2, MW3, MW4, MW5, and MW8) contributed very high soil erosion in the study area with 78.07% of the total area of the watershed (17.74 km2), three micro watersheds (MW6, MW7, and MW10) contributed high soil erosion with 20.64% of the total area of the watershed(17.74 km2), and one micro watershed (MW9) contributed medium soil erosion with 1.66% of the total area of the watershed (17.74 km2). Therefore, in the study area where soil erosion is significantly very high and we recommend the implementation of effective erosion control techniques that can contribute to long-term soil and water conservation efforts.

Similar content being viewed by others

Application of smart technologies for predicting soil erosion patterns

Article
Open access
21 July 2025

Impact of soil erosion on agricultural sustainability based on crop water productivity in semi-arid Iran

Article
Open access
18 November 2025

Assessing soil erosion and farmers’ decision of reducing erosion for sustainable soil and water conservation in Burji woreda, southern Ethiopia

Article
Open access
15 April 2024

Data availability

The datasets utilized and analyzed in this study are available from the corresponding author on reasonable request.
ReferencesAbdeta, G. C., Tesemma, A. B., Tura, A. L. & Atlabachew, G. H. Morphometric analysis for prioritizing sub-watersheds and management planning and practices in Gidabo Basin, Southern rift Valley of Ethiopia. Appl. Water Sci. 10 (7). https://doi.org/10.1007/s13201-020-01239-7 (2020).Kudnar, N. S. & Rajasekhar, M. A. Study of the morphometric analysis and cycle of erosion in Wainganga Basin, India. Model. Earth Syst. Environ. 6 (1), 311–327. https://doi.org/10.1007/s40808-019-00680-1 (2020).
Google Scholar 
Asfaw, D. & Neka, M. Factors affecting adoption of soil and water conservation practices: the case of wereillu woreda (District), South Wollo Zone, Amhara Region, Ethiopia. Int. Soil. Water Conser. Res. 5 (4), 273–279. https://doi.org/10.1016/j.iswcr.2017.10.002 (2017).
Google Scholar 
Woldemariam, G. W., Iguala, A. D., Tekalign, S. & Reddy, R. U. Spatial modeling of soil erosion risk and its implication for conservation planning: the case of the Gobele Watershed, East Hararghe Zone, Ethiopia. Land 7 (1), 25. https://doi.org/10.3390/land7010025 (2018).
Google Scholar 
Hossain, A. et al. Agricultural land degradation: processes and problems undermining future food security. Environ. Clim. Plant. Veg. Growth. https://doi.org/10.1007/978-3-030-49732-3-2 (2020).
Google Scholar 
Muralitharan, J., Francis, L. & Zerihun, D. Morphometric analysis and prioritization of Sub-watersheds for soil erosion using geomatics technologies in Megech river Catchment, lake Tana Basin, North Western Ethiopia. E J. S S D 8(1). https://doi.org/10.20372/ejssdastu:v8.i1.2021.225 (2021).Abera, A. D., Pritam, C. & Rinku, K. Prioritization of soil erosion-prone sub-watersheds using geomorphometric and statistical-based weighted sum priority approach in the middle Omo-Gibe river basin, Southern Ethiopia. Int. J. Digit. Earth 17(1).https://doi.org/10.1080/17538947.2024.2350198(2024).Okoli, J. et al. LIDAR-Derived DEM for landslide susceptibility assessment using AHP and fuzzy logic in Serdang, Malaysia. Geoscience 13 (2), 34. https://doi.org/10.3390/geosciences13020034 (2023).
Google Scholar 
Valkanou, K. et al. Assessment of neotectonic landscape deformation in evia Island, Greece, using GIS-based multi-criteria analysis. ISPRS Int. J. Geoinf. 10 (3), 118. https://doi.org/10.3390/ijgi10030118 (2021).
Google Scholar 
Meshram, S. G. et al. Assessing erosion prone areas in a watershed using interval rough-analytical hierarchy process(IR-AHP) and fuzzy logic (FL). Stoch. Environ. Res. Risk Assess. 36 (2), 297–312. https://doi.org/10.1007/s00477-021-02134-6 (2021).
Google Scholar 
Altaf, S., Meraj, G. & Romshoo, S. A. Morphometry and land cover based multi-criteria analysis for assessing the soil erosion susceptibility of the Western Himalayan watershed. Environ. Monit. Assess. 186 (12), 8391–8412. https://doi.org/10.1007/s10661-014-4012-2 (2014).
Google Scholar 
Haokip, P., Khan, M. A., Choudhari, P., Kulimushi, L. C. & Qaraev, I. Identification of erosion-prone areas using morphometric parameters, land use land cover and multi-criteria decision-making method: geo-informatics approach. Environ. Dev. Sustain. 24 (1), 527–557. https://doi.org/10.1007/s10668-021-01452-7 (2021).
Google Scholar 
Kushwaha, N. L. & Yousuf, A. Soil erosion risk mapping of watersheds using RUSLE, remote sensing and GIS: a review. Res. J. Agric. Sci. 8 (2), 269–277 (2017).
Google Scholar 
Duressa, A. A. et al. Identification of soil erosion prone areas for effective mitigation measures using combined approach of morphometric analysis and geographical information system. Results Eng. https://doi.org/10.1016/j.rineng.2023.101712 (2024).
Google Scholar 
Smith, J. & Johnson, A. A novel approach for prioritizing erosion-prone areas using land cover analysis. Environ. Manag. 45 (3), 321–335. https://doi.org/10.1007/s40899-024-01134-y (2020).
Google Scholar 
Mohamed, E. Watershed delineation and morphometric analysis using remote sensing and GIS mapping techniques in Qena-Safaga-Bir Queh, central Eastern desert. Int. J. Water Res. Environ. Eng. 12 (2), 22–46. https://doi.org/10.5897/IJWREE2019.0896 (2020).
Google Scholar 
Tamir, A., Aramde, F. & Million, G. Morphometric analysis: sub-watershed prioritization in the Temcha watershed, upper blue nile Basin, Ethiopia. Sustain. Water Resour. Manag. https://doi.org/10.1007/s40899-024-01134-y (2024). 10,155.
Google Scholar 
Negash, D. A., Moisa, M. B., Merga, B. B., Sedeta, F. & Gemeda, D. O. Soil erosion risk assessment for prioritization of sub watershed: the case of Chogo watershed, Horo Guduru Wollega Ethiopia. Environ. Earth Sci. 80:589.https://doi.org/10.1007/s12665-021-09901-2(2021).Gutema, D., Kassa, T. & Sifan, A. K. Morphometric analysis to identify erosion prone areas on the upper blue nile using GIS (case study of Didessa and Jema sub-basin, Ethiopia). Int. Res. J. Eng. Technol. 4(8). (2017).Aher, P. D., Adinarayana, J. & Gorantiwar, S. D. Quantifcation of morphometric characterization and prioritization for management planning in semi-arid tropics of india: a remote sensing and GIS approach. J. Hydrol. 511, 850–860. https://doi.org/10.1016/j.jhydrol.2014.02.028 (2014).
Google Scholar 
Peel, M. C., Finlayson, B. L. & Mahon, T. A. Updated world map of the Koppen-Geiger climate classification. Hydrol. Earth Syst. Sci. Discuss. 4, 439–473. https://doi.org/10.5194/hessd-4-439-2007 (2007).
Google Scholar 
NMSA (National metrological Service Agency) National metrological service agency at Dangila metrological station, Ethiopia annual report (2023).Food and Agriculture Organization of the United Nations (FAO). World Reference Base for Soil resources, 2006: a Framework for International classification, correlation, and Communication (Food and Agriculture Organization of the United Nations, 2006).Working Group, W. R. B. & IUSS World Reference Base for Soil Resources 2014, Update 2015 International Soil Classification System for Naming Soils and Creating Legends for Soil Maps (World Soil Resources Reports No. 106. FAO, 2015).Benzougagh, B. et al. Identification of critical watershed at risk of soil erosion using morphometric and geographic information system analysis. Appl. Water Sci. 12 (1). https://doi.org/10.1007/s13201-021-01532-z (2022).Kulimushi, L. C., Bashagaluke, J. B., Choudhari, P., Masroor, M. & Sajjad, H. Novel combination of analytical hierarchy process and weighted sum analysis for watersheds prioritization. A study of ulindi catchment, congo river basin. Geocarto Int. 37 (25), 8456–8494. https://doi.org/10.1080/10106049.2021.2002426 (2021).
Google Scholar 
Kadam, A. K., Jaweed, T. H., Kale, S. S., Umrikar, B. N. & Sankhua, R. N. Identification of erosion-prone areas using modified morphometric prioritization method and sediment production rate: a remote sensing and GIS approach. Geomat Nat. Hazards Risk 10(1), 986–1006. https://doi.org/10.1080/19475705.2018.1555189 (2019).
Google Scholar 
Nookaratnam, K., Srivastava, Y. K., Rao, V. V., Amminedu, E. & Murthy, K. S. Check dam positioning by prioritization micro-watersheds using SYI model and morphometric analysis—remote sensing and GIS perspective. J. Indian Soc. Remote Sens. 33, 25–38. https://doi.org/10.1007/BF02989988 (2005).
Google Scholar 
Das, D. Identification of erosion prone areas by morphometric analysis using GIS. J. Inst. Eng. (India). 95 (1), 61–74. https://doi.org/10.1007/s40030-014-0069-8 (2014).
Google Scholar 
Malik, A., Kumar, A. & Kandpal, H. Morphometric analysis and prioritization of sub-watersheds in a hilly watershed using weighted sum approach. Arab. J. Geosci. 12(4), 118 .https://doi.org/10.1007/s12517-019-4310-7(2019).Shekar, P. R. et al. Prioritizing sub-watersheds for soil erosion using Geospatial techniques based on morphometric and hypsometric analysis: a case study of the Indian Wyra river basin. Appl. Water Sci. 13 (7). https://doi.org/10.1007/s13201-023-01963-w (2023).Strahler, A. N. Quantitative Geomorphology of Drainage Basins and Channel Networks, In: Chow, V., Ed., Handbook of Applied Hydrology 439–476 (McGraw Hill, New York, 1964).Horton, R. E. Erosional development of streams and their drainage basins: hydropysical approach to quantitative morphology. Bull. Geol. Soc. Am. 56(2), 75–370 (1945). (Progress in Physical Geography: Earth and Environment, 19(4), 533–554). https://doi.org/10.1177/030913339501900406Strahler, A. N. Quantitative analysis of watershed geomorphology, eco transects. AGU Union. 38 (6), 913–920. https://doi.org/TR0381006P00913 (1957).
Google Scholar 
Schumm, S. A. Evolution of drainage systems and slopes in badlands at Perth Amboy New Jersey. Geol. Soc. Am. Bull. 67, 597–646. https://doi.org/10.1130/0016-7606 (1956).
Google Scholar 
Faniran A. The index of drainage intensity provisional new drainage factor. Aust. J. Sci. 31, 328–330 (1968).
Google Scholar 
Suresh, M., Sudhakar, S., Tiwari, K. N. & Chowdary, V. M. Prioritization of watersheds using morphometric parameters and assessment of surface water potential using remote sensing. J. Ind. Soc. Remote Sens. 32, 249–259. https://doi.org/10.1007/BF03030885 (2004).
Google Scholar 
Miller, V. C. A Quantitative Geomorphic Study of Drainage Basin Characteristics in the Clinch Mountain Area, Virginia and Tennessee 389–402 (Department of Geology Columbia University, 1953).Arefin, R., Mohir, M. M. I. & Alam, J. Watershed prioritization for soil and water conservation aspect using GIS and remote sensing: PCA-Based approach at Northern elevated tract Bangladesh. Appl. Water Sci. 10 (4), 1–19. https://doi.org/10.1007/s13201-020-1176-5 (2020).
Google Scholar 
Wischmeier, W. H. & Smith, D. D. Predicting Rainfall Erosion Losses.A Guide to Conservation Planning the USDA Agricultural Handbook No. 537, Maryland (1978).Prabhakaran, A. & Jawahar, R. N. Drainage morphometric analysis for assessing form and processes of the watersheds of Pachamalai hills and its Adjoinings, central Tamil Nadu, India. Appl. Water Sci. 8 (1), 1. https://doi.org/10.1007/s13201-018-0646-5 (2018).
Google Scholar 
Iqbal, M. & Sajjad, H. Watershed prioritization using morphometric and land use/land cover parameters of Dudhganga catchment Kashmir Valley India using Spatial technology. J. Geophy Remote Sens. 3 (1), 1–12. https://doi.org/10.4172/2169-0049.1000115 (2014).
Google Scholar 
Brahim, B. et al. Identification of critical watershed at risk of soil erosion using morphometric and geographic information system analysis. Appl. Water Sci. https://doi.org/10.1007/s13201-021-01532-z (2022).
Google Scholar 
Shekar, P. R. & Mathew, A. Morphometric analysis of watersheds: a comprehensive review of data sources, quality, and geospatial techniques. Watershed Ecol. Environ. https://doi.org/10.1016/j.wsee.2023.12.001 (2023).
Google Scholar 
Mukaka, M. M. Statistics corner: a guide to appropriate use of correlation coefficient in medical research. Malawi Med. J. 24, 69–71 (2012).
Google Scholar 
Santosh, W. & Vivek, M. A GIS-based morphometric prioritization of watersheds for soil erosion planning: a case study. Environ. Earth Sci. https://doi.org/10.1007/s12665-023-11155-z (2023). 82,443.
Google Scholar 
Bashir, B. & Alsalman, A. Geospatial analysis for tectonic assessment and soil erosion prioritization: a case study of Wadi Al-Lith, red sea Coast, Saudi Arabia. Appl. Sci. 13 (22), 12523. https://doi.org/10.3390/app132212523 (2023).
Google Scholar 
Khalifa, A., Bashir, B., Alsalman, A. & Bachir, H. Morphometric-Hydro characterization of the coastal line between El-Qussier and Marsa-Alam, Egypt: preliminary flood risk signatures. Appl. Sci. 12 (12), 6264. https://doi.org/10.3390/app12126264 (2023).
Google Scholar 
Shekar, P. R. & Mathew, A. Morphometric analysis for prioritizing sub-watersheds of Murredu river basin, Telangana State, India, using a geographical information system. J. Eng. Appl. Sci. 69 https://doi.org/10.1186/s44147-022-00094-4 (2022). :44.Gajbhiye, S., Mishra, S. K. & Pandey, A. Prioritizing erosion-prone area through morphometric analysis: an RS and GIS perspective. Appl. Water Sci. 4 (1), 51–61. https://doi.org/10.1007/s13201-013-0129-7 (2013).
Google Scholar 
Arabameri, A., Pradhan, B., Pourghasemi, H. R. & Rezaei, K. Identification of erosion-prone areas using different multicriteria decision-making techniques and GIS. Geomat. Nat. Hazards. 9 (1), 1129–1155. https://doi.org/10.1080/19475705.2018.1513084 (2018).
Google Scholar 
Mohammed, J. A., Gashaw, T. & Yimam, Z. A. Identification of erosion-prone watersheds for prioritizing soil and water conservation in a changing climate using morphometric analysis and GIS. Nat. Hazards. 121, 4171–4189. https://doi.org/10.1007/s11069-024-06952-z (2024).
Google Scholar 
Klingenberg, C. P. & Size Shape, and form: concepts of allometry in geometric morphometrics. Dev. Genes Evol. 226 (3), 113–137. https://doi.org/10.1007/s00427-016-0539-2 (2016).
Google Scholar 
Sahu, N. et al. Morphometric analysis in basaltic terrain of central India using GIS techniques: a case study. Appl. Water Sci. https://doi.org/10.1007/s13201-016-0442-z (2016).
Google Scholar 
Magesh, N. S., Jitheshlal, K. V., Chandrasekar, N. & Jini, K. Geographical information system-based morphometric analysis of Bharathapuzha river basin Kerala. India Appl. Water Sci. 3, 467–477. https://doi.org/10.1007/s13201-013-0095-0 (2013).
Google Scholar 
Asode, A. N., Sreenivasa, A. & Lakkundi, T. K. Quantitative morphometric analysis in the hard rock Hirehalla sub-basin, Bellary and Davanagere Districts, Karnataka, India, using RS and GIS. Arab. J. Geosci. 9 (6), 1–14. https://doi.org/10.1007/s12517-016-2414-x (2016).
Google Scholar 
Fenta, A. A., Hiroshi, Y., Katsuyuki, S., Haregeweyn, N. & Woldearegay, K. Quantitative analysis and implications of drainage morphometry of the Agula watershed in the semi-arid Northern Ethiopia. Appl. Water Sci. 7, 3825–3840. https://doi.org/10.1007/s13201-017-0534-4 (2017).
Google Scholar 
Javarayigowda, N. H., Basavaraju, G. K. S. & Jayaram, S. H. Morphometric analysis of Karadya micro watershed: a case study of Mandya district. Am. J. Remote Sens. 6 (1), 15–22. https://doi.org/10.11648/j.ajrs.20180601.13 (2018).
Google Scholar 
Harsha, J. & Shivakumar, A. Evaluation of morphometric parameters and hypsometric curve of Arkavathy river basin using RS and GIS techniques. Appl. Water Sci. 10 (3), 1–15. https://doi.org/10.1007/s13201-020-1164-9 (2020).
Google Scholar 
Download referencesAcknowledgementsThe authors acknowledge Bahir Dar University and Dilla University for financial support to conduct this study.FundingOpen access funding provided by Bahir Dar University and Dilla University.Author informationAuthors and AffiliationsDepartment of Natural Resources Management, College of Agriculture and Environmental Science, Bahir Dar University, P.O. Box 79, Bahir Dar, EthiopiaGetu Abey Denekewu & Enyew Adgo TsegayeDepartment of Natural Resources Management, College of Agriculture and Natural Resources, Dilla University, P.O. Box 419, Dilla, EthiopiaGetu Abey DenekewuArid Land Research Center, Tottori University, Tottori, JapanDerege Tsegaye MesheshaAuthorsGetu Abey DenekewuView author publicationsSearch author on:PubMed Google ScholarDerege Tsegaye MesheshaView author publicationsSearch author on:PubMed Google ScholarEnyew Adgo TsegayeView author publicationsSearch author on:PubMed Google ScholarContributionsGetu Abey: Conceptualization, Methodology, Investigation, Formal analysis, writing (original draft), writing (review & editing) and validation. Derege Tsegaye Meshesha and Enyew Adgo Tsegaye: Supervision, methodology, editing, writing, review, validation and visualization.Corresponding authorCorrespondence to
Getu Abey Denekewu.Ethics declarations

Competing interests
The authors declare no competing interests.

Additional informationPublisher’s noteSpringer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.
Reprints and permissionsAbout this articleCite this articleDenekewu, G.A., Meshesha, D.T. & Tsegaye, E.A. Morphometric analysis and weighted sum priority based prioritization of micro-watershed, kiltie watershed, Upper Blue Nile basin, Ethiopia.
Sci Rep (2025). https://doi.org/10.1038/s41598-025-32595-6Download citationReceived: 13 July 2025Accepted: 11 December 2025Published: 18 December 2025DOI: https://doi.org/10.1038/s41598-025-32595-6Share this articleAnyone you share the following link with will be able to read this content:Get shareable linkSorry, a shareable link is not currently available for this article.Copy shareable link to clipboard
Provided by the Springer Nature SharedIt content-sharing initiative
KeywordsMorphometric parametersWeighted sum prioritySRTM 30 mMicro watershed prioritization More

AbstractPolychlorinated biphenyls (PCBs) remain a global concern due to their environmental persistence and toxicity. However, their distribution in industrial and residential areas in Nigeria is insufficiently documented. This study investigated PCB concentrations in soils and commonly consumed plants from five sites around industrial areas in Koko Town, Delta State. Soil and plant samples were extracted using a Soxhlet extraction method and analyzed via GC-MS following the USEPA method 3540 C. Mean PCB concentrations were significantly higher in plants (20.75 mg kg− 1) than in soils (10.32 mg kg− 1), with both matrices exceeding the WHO recommended limits. PCB accumulation was highest in Pueraria phaseoloides, followed by Ceiba pentandra, Chromolaena odorata, Vermonia amygdalina, and Musa sapientum. The estimated daily intake (EDI) for adults and children exceeded the USEPA reference dose (0.000007 mg kg− 1) with a hazard ratio > 1, indicating a notable health risk, particularly for children. The findings underscore the need for regular monitoring and mitigative strategies for communities vulnerable to PCB pollution.

Similar content being viewed by others

Association between polychlorinated biphenyl (PCB) and dioxin with metabolic syndrome (METS): a systematic review and meta-analysis

Article
Open access
02 August 2024

Organochlorine pesticide residues in plants and their possible ecotoxicological and agri food impacts

Article
Open access
08 September 2021

Decrypting bacterial polyphenol metabolism in an anoxic wetland soil

Article
Open access
29 April 2021

Data availability

Data are available upon request from the corresponding author.
ReferencesAziza, A. E., Iwegbue, C. M. A., Tesi, G. O., Nwajei, G. E. & Martincigh, B. S. Concentrations, sources and exposure risk of polychlorinated biphenyls in soil profiles of the floodplain of the lower reaches of the River Niger, Nigeria. Environ. Monit. Assess. 193(9), 579. https://doi.org/10.1007/s10661-021-09310-9 (2021).
Google Scholar 
Jing, R., Fusi, S. & Kjellerup, B. V. Remediation of polychlorinated biphenyls (PCBs) in contaminated soils and sediment: state of knowledge and perspectives. Front. Environ. Sci. 6 https://doi.org/10.3389/fenvs.2018.00079 (2018).Zhang, C. et al. Uptake and translocation of organic pollutants in plants: A review. J. Integr. Agric. 16(8), 1659–1668. https://doi.org/10.1016/S2095-3119(16)61590-3 (2017).
Google Scholar 
Fernandez-Gonzalez, R., Yebra-Pimental, I., Martinez-Carballo, E. & Simal-Gandara, J. A critical review about human exposure to polychlorinated dibenzo-p-dioxins (PCDDs); polychlorinated dibenzofurans (PCDFs) and polychlorinated biphenyls (PCBs) through food. Crit. Rev. Food Sci. Nutr. 55(11), 1590–1617. https://doi.org/10.1080/10408398.2012.710279 (2015).
Google Scholar 
Iwegbue, C. M. A., Bebenimibo, E., Tesi, G. O., Egobueze, F. E. & Martincigh, B. S. Spatial characteristics and risk assessment of PCBs in surficial sediments around crude oil production facilities in the Escravos River Basin, Niger Delta, Nigeria. Mar. Pollut Bull. 159, 111462. https://doi.org/10.1016/j.marpolbul.2020.111462 (2020).
Google Scholar 
Arshad, M. et al. Monitoring of level of mean concentration and toxicity equivalence (TEQ) of polychlorinated biphenyls (PCBs) to selected vegetables, beans and grains to Khanewal and Multan, Pakistan. Saudi J. Biol. Sci. 29(4), 2787–2793. https://doi.org/10.1016/j.sjbs.2022.01.009 (2022).
Google Scholar 
Olatunji, O. S. Evaluation of selected polychlorinated biphenyls (PCBs) congeners and dichlorodiphenyltrichloroethane (DDT) in fresh root and leafy vegetables using GC-MS. Sci. Rep. 9, 538. https://doi.org/10.1038/s41598-018-36996-8 (2019).
Google Scholar 
Eze, T. C. & Eze, A. G. Control of pollution arising from oil and gas industry: appraising the scope of provisions under the 1999 Nigerian Constitution. Nnamdi Azikiwe University Journal of International Law and Jurisprudence. 8(2), 58–60 (2017).Michael, A., Ekperusi, O. A., Okeke, P. N., Ihejirika, C. E. & Ejiogu, C. C. Assessment of dioxin concentration in soil and edible plants around industrial sites in Koko Town, Southern Nigeria. In SPE Nigeria Annual International Conference and Exhibition D031S020R006 https://doi.org/10.2118/228672-MS (2025).Plank, V. W. & Plant Sampling Mineral Nutrition-TNAU Agritech portal. https://agritech.tnau.ac.in/agriculture/agri_min_nutri_plantsampling.html (1979).USEPA. SW-846 Test Method 3540 C. Soxhlet Extraction part of test methods for evaluating solid waste, physical/chemical methods. https://www.epa.gov/hw-sw846/sw-846-test-method-3540c-soxhlet-extraction (2022).FAO. FAOSTAT – Food Balances (2010–2022). Accessed July 2024. https://www.fao.org/faostat/en/#data/FBS (2022).Walpole, S. C. et al. The weight of nations: an estimation of adult human biomass. Public. Health. 12, 439. https://doi.org/10.1186/1471-2458-12-439 (2012).
Google Scholar 
Goon, D. T. et al. Anthropometrically determined nutritional status of urban primary schoolchildren in Makurdi, Nigeria. BMC Public. Health. 11, 769–776. https://doi.org/10.1186/1471-2458-11-769 (2011).
Google Scholar 
USEPA. Regional Screening Level (RSL) Summary Table (TR = 1E-06, HQ = 1) (2024).DeVito, M. et al. The 2022 world health organization re-evaluation of human and mammalian toxic equivalency factors for polychlorinated dioxins, dibenzofurans, and biphenyls. Regul. Toxicol. Pharmacol. 146, 105525. https://doi.org/10.1016/j.yrtph.2023.105525 (2024).
Google Scholar 
Lacomba, I. et al. Levels and risk assessment of dl-PCBs and dioxins in soils surrounded by cement plants from industrial areas of Colombia and Spain. Emerg. Contaminants. 11(1), 100427. https://doi.org/10.1016/j.emcon.2024.100427 (2025).
Google Scholar 
Sandu, M. A. et al. Trends in polychlorinated biphenyl contamination in Bucharest’s urban soils: A two-decade perspective (2002–2022). Processes 13(5), 1357. https://doi.org/10.3390/pr13051357 (2025).
Google Scholar 
Košnář, Z. & Tlustoš, P. Translocation and dissipation of seven indicator polychlorinated biphenyls from contrast soils cultivated with different root vegetables. Environ. Sci. Eur. 36, 177. https://doi.org/10.1186/s12302-024-01006-4 (2024).
Google Scholar 
Esposito, M. et al. Occurrence of polychlorinated dibenzo-p-dioxin and dibenzofurans and polychlorinated biphenyls in fruit and vegetables from the land of fires. Area South. Italy Toxics. 5(4), 33. https://doi.org/10.3390/toxics5040033 (2017).
Google Scholar 
Tuomisto, J. Dioxins and dioxin-like compounds: toxicity in humans and animals, sources and behaviour in the environment. Wiki J. Med. 6(1), 8. https://doi.org/10.15347/wjm/2019.008 (2019).
Google Scholar 
Irerhievwie, G. O. et al. Spatial characteristics, sources and ecological and human health risks of polychlorinated biphenyls in sediments from some river systems in the Niger Delta, Nigeria. Mar. Pollut Bull. 160, 111605. https://doi.org/10.1016/j.marpolbul.2020.111605 (2020).
Google Scholar 
Wang, Z. et al. Removal of cadmium and polychlorinated biphenyls by clover and the associated microbial community in a Long-term co-contaminated soil. Sci. Total Environ. 871. https://doi.org/10.1016/j.scitotenv.2023.161983 (2023).Michael, A., Okeke, P. N., Ihejirika, C. E. & Ejiogu, C. C. Assumptions on health risks in consuming vermonia amygdalina and fruits (Musa sp.) in Koko, Nigeria. EJFOOD 5(5), 12–15. https://doi.org/10.24018/ejfood.2023.5.5.711 (2023).
Google Scholar 
Bantum, J., Dodoo, D., Kwakye, P., Essumang, D. & Adjei, G. Spatial and temporal distribution of polychlorinated biphenyl residues in tropical soils. Open. J. Appl. Sci. 6, 234–247. https://doi.org/10.4236/ojapps.2016.64024 (2016).
Google Scholar 
Eghbaljoo, H. et al. Analysis of polychlorinated biphenyls (PCBs) in edible oils using the QuEChERS/GC-MS method: A health risk assessment study. Heliyon 9(11), e21317. https://doi.org/10.1016/j.heliyon.2023.e21317 (2023).
Google Scholar 
Kumar, B. et al. Polychlorinated biphenyls in residential soils and their health risk and hazard in an industrial city in India. J. Public. Health Res. 3(2), 252. https://doi.org/10.4081/jphr.2014.252 (2014).
Google Scholar 
Li, X. & Su, X. Assessment of the polychlorinated biphenyl (PCB) occurrence in copper sulfates and the influential role of PCB levels on grapes. PLoS ONE. 10(12), e0144896. https://doi.org/10.1371/journal.pone.0144896 (2015).
Google Scholar 
WHO/EURO. PCBs, PCDDs, and PCDFs. Prevention and Control of Accidental and Environmental exposures. Copenhagen 227 Environmental Health Series No. 23 https://apps.who.int/iris/bitstream/10665/39892/1/9241510684_eng.pdf (World Health Organization Regional Office for Europe, 1987).Boakye-Yiadom, M., Kumadoh, D., Adase, E. & Woode, E. Medicinal plants with prospective benefits in the management of peptic ulcer diseases in Ghana. Biomed. Res. Int. 2021, 5574041. https://doi.org/10.1155/2021/5574041 (2021).Habtamu, A. & Melaku, Y. Antibacterial and Antioxidant Compounds from the Flower Extracts of Vermonia amygdalina. Adv. Pharmacol. Pharm. Sci. 2018(1), 4083736. https://doi.org/10.1155/2018/4083736 (2018).Download referencesAcknowledgementsThe authors are grateful to the Herbarium Unit, University of Benin where identification of plant species was carried out by Amaka Michael and Akinnibosun H. Adewale.Author informationAuthors and AffiliationsDepartment of Environmental Management and Pollution, Faculty of Environmental Management, Nigeria Maritime University, Okerenkoko, Delta State, NigeriaAmaka Michael, Abraham O. Ekperusi & Anthonia E. GbuvboroDepartment of Environmental Management, School of Environmental Sciences, Federal University of Technology, Owerri, Imo State, NigeriaNdu P. OkekeAuthorsAmaka MichaelView author publicationsSearch author on:PubMed Google ScholarAbraham O. EkperusiView author publicationsSearch author on:PubMed Google ScholarNdu P. OkekeView author publicationsSearch author on:PubMed Google ScholarAnthonia E. GbuvboroView author publicationsSearch author on:PubMed Google ScholarContributionsAmaka Michael, prepared the manuscript, interpreted the results and wrote the paper, Abraham O. Ekperusi revised the work, Ndu P. Okeke initiated the study, and Anthonia E. Gbuvboro wrote the summary. All authors read and approve the final manuscript.Corresponding authorCorrespondence to
Amaka Michael.Ethics declarations

Competing interests
The authors declare no competing interests.

Additional informationPublisher’s noteSpringer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.Supplementary InformationBelow is the link to the electronic supplementary material.Supplementary Material 1Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.
Reprints and permissionsAbout this articleCite this articleMichael, A., Ekperusi, A.O., Okeke, N.P. et al. Distribution and human health risk of polychlorinated biphenyls in soil and plants in Koko Town, Delta State, Nigeria.
Sci Rep (2025). https://doi.org/10.1038/s41598-025-33241-xDownload citationReceived: 07 July 2025Accepted: 17 December 2025Published: 18 December 2025DOI: https://doi.org/10.1038/s41598-025-33241-xShare this articleAnyone you share the following link with will be able to read this content:Get shareable linkSorry, a shareable link is not currently available for this article.Copy shareable link to clipboard
Provided by the Springer Nature SharedIt content-sharing initiative
KeywordsPolychlorinated biphenylsBioaccumulationNiger deltaHuman exposureRisk assessment More

AbstractSophora moorcroftiana (Benth.) Baker is a crucial shrub species for windbreak and soil conservation in Tibet (Xizang), playing an important role in plateau ecological protection and the biopharmaceutical industry. However, there have been no studies reporting on the phenological characteristics or effective accumulated temperature of this species. This research systematically describes the developmental processes of S. moorcroftiana across various phenological growth stages using the BBCH scale (Biologische Bundesanstalt, Bundessortenamt, and Chemische Industrie). Through long-term photography, tracking, and observation, eight main phenological stages were identified: bud development, leaf development, shoot development, inflorescence emergence, flowering, fruit development, seed maturation, and senescence and beginning of dormancy. Additionally, 41 secondary growth stages were detailed, accompanied by characteristic images, standardizing morphological features and phenological observation criteria for S. moorcroftiana. This study provides a scientific reference for research on the biological characteristics and breeding of superior varieties of S. moorcroftiana.

Similar content being viewed by others

Impacts of climate change on reproductive phenology in tropical rainforests of Southeast Asia

Article
Open access
21 April 2022

Persistent microbiome members in the common bean rhizosphere: an integrated analysis of space, time, and plant genotype

Article
Open access
26 March 2021

Divergent phenological responses of soil microorganisms and plants to climate warming

Article

29 July 2025

IntroductionSophora moorcroftiana (Benth.) Baker (Leguminosae: Sophora) is a deciduous dwarf shrub widely distributed along riverbanks, valleys, and hillsides in the Yarlung Tsangpo, Nyenchu, and Lhasa Rivers in Tibet1. This species serves as an important pioneer species for windbreak, sand fixation and soil conservation in the region. It exhibits drought tolerance, soil adaptation and robust environmental adaptability, playing a key role in the plateau ecosystem2,3. Current research predominantly focuses on its medicinal properties and drought resistance mechanisms. Yin, X. et al.4,5 found that S. moorcroftiana seeds are rich in flavonoids, matrine-type alkaloids, terpenoids, steroids, and their derivatives. Matrine-type alkaloids exhibit a wide spectrum of pharmacological activities, including antiviral, antifibrotic, antitumor, leukocytosis-promoting, and immunomodulatory effects6,7. The tender branches, leaves, and seeds of S. moorcroftiana also have high nutritional value, with the seeds containing substantial protein. This makes the species an important source of nutrition in the protein-deficient regions of Tibet and a high-quality forage resource for plateau livestock8.Phenology is the study of periodic biological events in relation to environmental changes. The BBCH scale (Biologische Bundesanstalt, Bundessortenamt, and Chemische Industrie) provides a standardized coding system that describes comparable phenological stages in both monocotyledonous and dicotyledonous plants. It was developed through the collaboration among scientists from the German Federal Biological Research Centre for Agriculture and Forestry (BBA), the Federal Plant Variety Office (BSA), and the German Agricultural Chemical Industry (IVA)9,10. BBCH scale is now widely used for crops, forestry species, and medicinal plants, and has standardized descriptions of phenological periods and their time periods for all flowering plants11. Many woody plants are described by the BBCH, such as mango12, longan11, sugar apple13, sweet cherry14, and persimmon15. The BBCH scale can be used to characterize the different phases in a hierarchical manner, which makes it more suitable for reflecting the phenology of shrubs at high altitudes and facilitates direct comparisons with other species or taxonomic groups. In 1997, the BBCH scale was officially recommended for use by the European and Mediterranean Plant Protection Organization (EPPO), and the Global Phenological Monitoring Programme also adopt this scale as a standard for phenological observations9. To complete certain developmental stages, plants require a specific amount of thermal energy16, with effective accumulated temperature being the most effective method for estimating plant maturity17. Predicting plant maturity can reduce agricultural costs and guide optimal cultivation timing18.Although substantial progress in understanding the ecological, medicinal and forage values of S. moorcroftiana has been done, the species continues to face several challenges, including declining populations, underutilization of its resources, outdated cultivation techniques, and absence of a standardized framework for phenological characterization. To address these issues, we conducted a two-year systematic observation of the growth stages of S. moorcroftiana. A comprehensive phenology framework was also established for the plateau shrub S. moorcroftiana. Each developmental stage was described in detail, and the corresponding temperature requirements for each phenophases were analyzed. These analyses enabled the identification of the optimal seeds harvesting period, the growth stages most susceptible to insect infestation, as well as the appropriate temperature conditions for the cultivation of S. moorcroftiana. Collectively, these findings provide a scientific basis for conservation, breeding, introductions and biopharmaceutical exploration of S. moorcroftiana, and establish a robust phenology model for this ecologically and economically important species.Materials and methodsExperimental materialsThe experiment was conducted at the nursery base of the Xizang Agricultural and Animal Husbandry University in Nyingchi, Tibet Autonomous Region (94°20′40″E; 29°40′24″N). Nyingchi, located in southeastern Tibet along the middle reaches of the Yarlung Tsangpo River, has an average altitude of 2800 m. It has unique topographic features such as high mountains and valleys, and great diversity of landforms with huge vertical drop. Nyingchi is one of the main distribution areas of S. moorcroftiana. Climate data from the past 30 years, sourced from the Xizang Meteorological Administration, indicates that the region exhibits a plateau temperate humid-subhumid monsoon climate. The annual average temperature is 9.48 °C, with the highest temperature recorded at 16.49 °C (July) and the lowest at 1.15 °C (January). Its soil is sandy loam with high year-round surface temperatures, regular watering and irrigation, no towering irrigation or trees around, and excellent sun exposure conditions, providing favorable conditions for S. moorcroftiana growth.The plant material of S. moorcroftiana used in this study was identified by Prof. Fumei Xin from Xizang Agricultural and Animal Husbandry University. The identification was based on the herbarium (PE 02330018) provided by Xiangyun Zhu and deposited in the Chinese Virtual Herbarium (https://www.cvh.ac.cn). Our S. moorcroftiana plants (id:19944579) have been preserved and shared in the Plant Photo Bank of China (http://ppbc.iplant.cn/). This study used all S. moorcroftiana plants at the nursery to measure effective cumulative temperature, and selected twelve 6-year-old plants with the same height and crown width for phenological observations. Maintained by irrigation and weeding. Traits were recorded and photographed weekly or bi-weekly, and dormancy lasted for two years, from bud dormancy to pre-germination to flowering and fruiting. Climate data is sourced from the Xizang Meteorological Administration, with daily temperature monitoring conducted between 2023 and 2024. During this period, no abnormal climate conditions were observed, ensuring the data’s strong representativeness. Representative photographs of each phenological stage were selected, describing traits from bud dormancy to flowering, fruiting, and the next dormancy period, spanning two years. During the study period, the plants were kept healthy with consistent water and fertilizer management and no use of pesticides, external hormones.BBCH scaleThe BBCH scale was employed to classify the phenological stages of S. moorcroftiana. It uses a two-digit code: one digit for the primary growth stage (0–9) and another for secondary stages (0–9)12. This study identified eight of the ten primary BBCH stages, starting from bud development (0), leaf development (1), stem elongation (3), inflorescence emergence (5), flowering (6), fruit development (7), seed maturation (8), to senescence and dormancy (9). The coding system facilitates comparisons, with higher numbers indicating further progression within the same primary stage. Separate overlapping stages with diagonal strokes9.Effective accumulated temperaturePlant threshold temperature varies according to species and phenological period, and the threshold temperature of S. moorcroftiana has not been reported yet. For S. moorcroftiana, a threshold of 7 °C was established based on the average temperature of initial developmental stages and literature references19. We calculated the effective accumulated temperature using the meteorological dataset recorded by Weather Station No. 1 at the experimental site from 1991 to 2020. Phenological observations and photography documented the morphology at each stage, calculating the heat accumulation needed for each phase.$${text{GDD = }}sumlimits_{i = 1}^{n} {left( {T_{{text{i}}} – B} right)}$$
(1)
where: GDD ( °C·d) denotes growing cumulative temperature; n denotes the number of days elapsed during the development period; Ti denotes the average air temperature on day i; B = 7 ℃ denotes the biological zero degree of the development stage20.ResultsDescription of S. moorcroftiana phenological stages using BBCH CodesThrough observations and recordings of the phenological stages of S. moorcroftiana, this study identified eight major growth stages and their durations (Table 1, Fig. 1). The plant morphology of S. moorcroftiana is not amenable to description in terms of the second stage (formation of lateral buds or tillers) and the fourth stage (development of harvestable vegetative parts or vegetative reproductive organs/bud formation). The study documented changes in S. moorcroftiana phenology through photography, irrigation, and weeding, recording the status at regular intervals from bud dormancy until the next dormancy, lasting two years. We described the morphological characteristics of individual plants using the BBCH scale and found no phenological differences between them.Table 1 Description of S. moorcroftiana phenological stages according to the BBCH Scale.Full size tableFig. 1Phenological stages of S. moorcroftiana. Vertical axis represents temperature and precipitation, monthly average, maximum, and minimum temperatures and average precipitation in Nyingchi over 30 years (1991–2020). Horizontal axis indicating the duration of each phenological stage of S. moorcroftiana. Light blue (00–09) indicates bud development stage, green (10–19) indicates leaf development stage, blue (30–39) indicates shoot development stage, mauve (50–59) indicates inflorescence emergence stage, purple (60–69) indicates flowering stage, yellow (70–79) indicates fruit development stage, brown (80–89) indicates seed maturation stage, and gray (90–99) indicates senescence and dormancy stage.Full size imagePrimary growth stage 0: bud developmentBud dormancy at the observation site lasted until early March, with bud swelling beginning in mid-March, and flowering occurring at the end of bud development (Fig. 2).Fig. 2Vegetative growth stages of S. moorcroftiana according to BBCH classification.Full size image00. Bud dormancy: Leaf buds are grayish-white and closed in a dormant state.01. Bud swelling starts: The leaf buds begin to expand again in mid-March, and the buds are ready to unfold.03. Bud swelling ends: Buds reach full swelling and slightly open.07. Bud scale separation: Bud scales separate, revealing purple leaf tips.09. Visible purple bud tip: Purple tips become visible, protruding 5–10 mm above the grayish-white scales.Primary growth stage 1: leaf developmentLeaf development occurs from mid-April to mid-August, taking about four months for complete leaf maturation (Fig. 2).10. First leaf separation: Leaves turn from purple to gray-green, with the outermost leaf separating.11. First leaf unfolds: The leaflets of the first compound leaf open, with petioles extending to 10% of their final length.14. 40% of juvenile leaves unfolded: Juvenile leaves spread out, turning light green, and petioles lengthen.17. 70% of leaves unfolded: Leaf unfolds to final 70%, leaf is light green, leaf thickens, petiole hardens.19. All leaves fully unfolded: Leaf blades reach full size, petioles are elongated to their final size, and leaflets are mostly thickened to a dark green color.Primary growth stage 3: shoot developmentShoot development often occurs concurrently with leaf development. S. moorcroftiana shoot development over an extended period, ceasing during fruit development. Non-reproductive branch tips form spines (Fig. 2).30. Shoot begin to elongate: The bud axis becomes visible and starts developing as the first compound leaf separates. This stage coincides with stage 10.31. 10% of the final shoot length: The stem thickens and appears dark green, covered with white trichomes.32. 20% of the final shoot length: The stem thickens further, coinciding with stage 11.34. 40% of the final shoot length: The stem lengthens and turns gray-green as trichomes decrease.36. 60% of the final shoot length: The shoot continues to extend, with most branches forming spikes at the end, non-reproductive branch tips form spines.39. 90% or more of the final shoot length: The shoot hardens, completing development, turning green with almost no trichomes.Primary growth stage 5: inflorescence emergenceThe inflorescence development stage occurs from early April to mid-May, with flower buds and leaves developing simultaneously (Fig. 3).Fig. 3Reproductive growth stages of S. moorcroftiana according to BBCH classification.Full size image50. Flower buds begin to swell: The bracts are close to the buds, and the surface is densely covered with gray-white tomentum.51. Continued bud swelling: Bracts unfold, and purple flower buds reach 50% of their final size.52. Bud axis starts elongating: The bud axis extends to 20% of its final length, with small flowers densely arranged and dark purple flower buds covered with gray-white trichomes.54. Bud axis reaches 40% of final length: Small flowers slightly separate, and the color lightens.56. Bud axis reaches 60% of final length: Small flowers become purple and increase in size.59. Bud axis reaches full length: The inflorescence develops to its final size, with some small flowers emerging from the bracts and revealing purple petals.Primary growth stage 6: floweringThe flowering stage occurs from mid-May to early June. The flowers are purple, and the pods appear as the petals begin to fade (Fig. 3).62. Anthesis: Approximately 20% of flowers are open.64. Further flowering: More flowers open, reaching 40% of the final flowering stage.65. Full flowering: All flowers are fully open.67. Start of wilting: Petals turn yellow and begin to drop or dry out.69. End of flowering: Most of the petals fall off and dry out, the undersides of most petals are grayish-purple, and fruit pods appear.Primary growth stage 7: fruit developmentFruit development for S. moorcroftiana occurs from early June to early August. Each pod typically contains 1–5 seeds (Fig. 3).70. Seeds begin to form: Wilted petals encase 1–2 cm long, slightly flattened pods covered in white trichomes.74. Pods elongate: Pods grow longer, with wilted petals falling off completely, reaching 30% of their final length. The pods are grayish-white, and trichomes decrease.75. Pods further elongate: Pods reach 90% of their final length, turning light green, with seeds swelling and the remaining areas still covered with white trichomes.77. Seeds swell further: Pods reach full length, and seeds are dark green, swollen and oval, with reduced trichome density.79. Seed development complete: Pods turn gray-green, with brown spots on the swollen parts, and slightly wrinkle.Primary growth stage 8: seed maturationSeed maturation of S. moorcroftiana occurs from early August to mid-September (Fig. 3).81. Early seed maturation: Pods are gray-green, with fewer trichomes on the swollen seed areas, and seeds are tender and light green.82. Pods begin to change color: Pods turn yellow-green with brown spots appearing on the surface, and peduncles turn brown.85. Seeds mature further: Pods become yellow–brown, with reduced trichomes, and seeds turn light yellow.89. Full development and maturity: Pods dry out, turn brown, wrinkle slightly, and the tips of the pods crack, revealing yellow seeds.Primary growth stage 9: senescence and beginning of dormancyThe onset of senescence and dormancy occurs in tandem with the process of fruit ripening, and October to December is the stage of defoliation of S. moorcroftiana, with little variation among individual plants (Fig. 4).Fig. 4Senescence and dormancy stages of S. moorcroftiana.Full size image90. Bud and leaf development cease: The buds stop growing and the leaves remain dark green in color.91. Leaves begin to discolor: Leaflet tips yellowed with dark brown dots.95. Half of the leaves fall off: 50% of the leaves fall off, most of the leaves are yellow and some are still green.97. Easily or completely dislodged: Leaf blades are almost entirely lost, petioles remain on the tree but are easily shed.99. Dormancy: The trees completely enter into dormant period.Primary growth stage 2 and 4The BBCH standard is predicated on the developmental stage of the main stem, and given the high degree of uniformity in the development of the plant’s main stem, the growth stage of lateral branches is not usually described specifically. The temporal dynamics, quantitative development, and spatial distribution of lateral branches are found to be profoundly influenced by light, nutrients, and hormones. Furthermore, it is important to note that different lateral branches of the same plant may be at different stages of development. The incorporation of the developmental stage of lateral branches may serve to reduce the generalizability of the scale. Stage 4 of the BBCH generally describes asexual reproductive organs and harvestable nutrient parts such as root tillers, stolons and root systems, whereas asexual reproduction in S. moorcroftiana tends to be at the root tiller level, and harvestable root systems belong to the stage of nutrient growth that is not applicable to Stage 4 of the BBCH. Similar to Xanthoceras sorbifolium11, Sapindus mukorossi21, Spondias dulcis22 and shea tree23.Effective accumulated temperatureThe primary objective of phenological research is to link climatic events with specific phenological stages, providing an approximate prediction of these stages24. Temperature is a critical factor affecting plant development. Different growth stages of plants require varying effective cumulative temperature and durations. This study calculated the duration (days) and effective cumulative temperature (degree-days) required for each phenological stage of S. moorcroftiana from 2023 to 2024 using Eq. 1 (Table 2). The mean and standard deviation of the duration and effective temperature were calculated. The duration of the bud stage was 30 ± 1 d with an effective cumulative temperature of 48.95 ± 10.61 °C·d. The leaf development stage had the longest duration of 119 ± 6 d with an effective cumulative temperature of 1063.13 ± 112.32 °C·d. The inflorescence development stage lasted for 33 ± 1 d with an effective cumulative temperature of 147.98 ± 3.64 °C·d. The flowering stage lasted 21 ± 1 d with an effective temperature of 190.50 ± 23.69 °C·d. Fruit development was longer at 49 ± 2 d with an effective temperature of 534.00 ± 42.21 °C·d. And seed maturation lasted 38 ± 1 d with an effective temperature of 455.13 ± 10.92 °C·d.Table 2 Total duration and effective accumulated temperature of different growth stages (2023–2024).Full size tableDiscussionAs a cold-tolerant, drought-tolerant and barren-tolerant pioneer shrub, S. moorcroftiana thrives under harsh conditions in Tibet, with its development closely linked to the local ecological environment. Using the BBCH scale, we systematically examined the growth and developmental stages of S. moorcroftiana, filling the gap in previous research on its phenological characteristics and effective accumulated temperature requirements. This research contributes to a better understanding of the growth patterns of S. moorcroftiana and provides a solid foundation for ecological conservation, breeding of superior varieties, and resource utilization. We observed and documented the complete developmental loop of S. moorcroftiana, from bud dormancy to flowering, fruiting, and seed maturation, identifying eight major growth stages and 41 secondary growth stages. The BBCH scale offers a standardized tool for describing and comparing the growth stages of various plant species. The results indicated that, due to the low temperature in the growing environment, the buds of S. moorcroftiana remained dormant from November to early-March of the following year, with germination occurring approximately one month later than for other tree species11,21,23. S. moorcroftiana leaf buds develop at the same time, but their development time is longer, stopping at the early stage of fruit maturation, with new buds developing during fruit development12,25. At stage 34, the sterile branch ends of S. moorcroftiana developed into spikes. The systematic description of phenology provides comprehensive insights into S. moorcroftiana development, and the scale’s comparability enhances the utility of phenological observations globally. The use of the BBCH coding system in phenological studies is crucial.Plants are able to complete their growth and development under suitable temperature conditions, and cumulative temperature is an important parameter for measuring the rate of development and maturity of plants26. Based on the effective cumulative temperature data for 2023–2024, the duration and standard deviation of effective cumulative temperature (°C·d) for bud development, inflorescence formation, flowering, and fruit development stages are relatively small. In contrast, the leaf development stage requires the longest duration and the highest effective cumulative temperature. It is possible to provide better economic and ecological value by calculating the effective temperature and heat required by different species at different growth stages to select suitable places for introduction, such as grapes and olives17,25. Effective accumulated temperature aids in predicting developmental stages and optimizing management practices, such as fertilization, pruning, and pest control, to enhance growth efficiency and yield27. Additionally, S. moorcroftiana growth varies with altitude and climate, implying that local accumulated temperature conditions should be considered to increase survival rates and adaptability during introduction and cultivation28.As an important sand-fixing and wind-preventing plant, S. moorcroftiana has a positive impact on the ecological environment of the Tibetan region, and also shows great potential in the field of biomedicine. Studies have shown that its seeds and young leaves contain abundant active compounds, such as flavonoids, alkaloids, and terpenoids, with notable antiviral, antitumor, and immunomodulatory properties4,7. This study’s phenological observations provide a clear time frame for harvesting these medicinal components. Harvesting seeds at full maturity (BBCH stage 89) ensures maximum extraction efficiency of active ingredients29. The study also emphasised the challenges of biotic stress during seed maturity. Observations revealed severe pest infestation at stage 77, with more than 70% of the seeds being damaged by insects, such as Robinia pseudoacacia bee, which exhibits asynchronous emergence. These bees spend most of their developmental period inside the seeds, consuming the majority of the kernel while the seeds are still immature. Therefore, preventive measures are of primary importance. Chemical control should be applied during the seed ripening stage, typically from late August to early September of the preceding year. Additionally, fumigation should be carried out at stages 70 to 73 during the second half of June of the following year2. Accurate phenological monitoring enables early pest control, reducing damage and increasing seed yield. This phenology-based pest management strategy has proven effective in other crops, significantly reducing pesticide use and protecting the environment30.Although this study provided systematic data on the phenological characteristics of S. moorcroftiana, the observations lasted only two years and were not sufficient to fully reveal the effects of all growth environment variables on its development. Therefore, follow-up studies should extend the observation period to cover the precise phenological stage. Multi-point data from various distribution areas and altitudes of S. moorcroftiana should be incorporated to validate and supplement the phenological model and heat demand specified in this study. It is evident that observations pertaining to insect pests were not meticulously documented, thereby precluding the possibility of proposing specific control measures. In addition, the establishment of the phenological period of S. moorcroftiana should be combined with long-term meteorological data, especially with the comprehensive analysis of other local climatic characteristics, in order to improve the accuracy of its acclimatization study and phenological prediction in the plateau region. Future research should concentrate on investigating the dynamics of alkaloid content in S. moorcroftiana plants at differing developmental stages. Furthermore, the optimal harvesting period to enhance medicinal value should be ascertained. Additionally, the genetic diversity of S. moorcroftiana and its influence on phenology should be examined. Finally, molecular biology techniques should be utilised to study the relationship between its growth and environmental adaptability. Given intensifying global climate change, the applicability of these results to field conditions requires verification. Continuous monitoring of S. moorcroftiana and other plateau plants’ phenology will provide crucial insights for ecological conservation and resource management.ConclusionThis study systematically describes eight primary and 41 secondary phenological stages of S. moorcroftiana using the BBCH scale, filling a significant gap in phenological research for this species. Following two years of uninterrupted observation and temperature analysis, the effective cumulative temperatures required for each developmental stage were determined, and the critical role of temperature in growth was demonstrated. The results of the study demonstrate the applicability and comparability of the BBCH scale in the context of highland plant research. This standardized climatic description and effective cumulative temperature can predict the growth and development process of S. moorcroftiana, thus guiding the appropriate harvesting period, determining suitable areas for planting, and improving the accuracy of pest control through BBCH. Furthermore, it provides a scientific basis for field management, collection of medicinal components, and selection of superior varieties of S. moorcroftiana. Additionally, as an important plant for ecological restoration and biomedical resources, recognizing key phenological stages in S. moorcroftiana’s life cycle will enhance ecological and economic benefits. In the future, there is still a need to combine long-term observations and broader regional studies to further improve the phenology model of S. moorcroftiana and optimize resource use and conservation strategies.

Data availability

All data are available from the corresponding author upon reasonable request.
ReferencesWang, W. J., He, D. H. & Liu, J. H. Current status of research on Sophora moorcroftiana in Tibet and evaluation of its utilization value. Anhui Agric. Sci. 36(33), 14513–14515 (2008).
Google Scholar 
Zang, J. C., Xin, F. M. & Wang, Z. H. A study on the damage and control of Robinia pseudoacacia bee on Sophora moorcroftiana seeds. Anhui Agric. Sci. 36(32), 14179–14180 (2008).
Google Scholar 
Lin, S. M. Characterization of seed germination of Sophora moorcroftiana in Tibet. Grass Sci. 05, 30–32 (2002).
Google Scholar 
Yin, X. et al. Sophora moorcroftiana genome analysis suggests association between sucrose metabolism and drought adaptation. Plant physiol. 191(2), 844–848 (2022).Article 
PubMed Central 

Google Scholar 
Li, H., Yao, W. J., Fu, Y. R., Li, S. K. & Guo, Q. Q. De novo assembly and discovery of genes that are involved in drought tolerance in Tibetan Sophora moorcroftiana. PLoS ONE 10(1), e111054. https://doi.org/10.1371/journal.pone.0111054 (2015).Article 
CAS 
PubMed 
PubMed Central 
ADS 

Google Scholar 
Wang, Y. T., Pubu, C. R., Ma, W. J. & Danzeng, L. B. Seed alkaloid content and its correlation analysis of Sophora moorcroftiana populations from different growing sites in Tibet. Northwest J Botany 38(10), 1913–1917 (2018).
Google Scholar 
Li, Z. G. Progress of pharmacological and clinical studies on matrine. West China J. Pharm. 06, 435–437 (2003).
Google Scholar 
Xue, J. Z. et al. Brief introduction and research progress on feeding value of Caragana korshinskii, Haloxylon ammodendron and Sophora moorcroftiana. Anim. Feed Sci. 39(10), 40–43 (2018).
Google Scholar 
Meier, U. et al. The BBCH system to coding the phenological growth stages of plants-history and publications. J. für Kulturpflanzen 61(2), 41–52 (2009).
Google Scholar 
Hack, H. et al. Einheitliche codierung der phänologischen entwicklungsstadien mono-und dikotyler pflanzen–erweiterte BBCH-Skala, Allgemein. Nachrichtenblatt des Deutschen Pflanzenschutzdienstes 44(12), 265–270 (1992).
Google Scholar 
Jiang, Y. G. et al. Phenological growth stages of Xanthoceras sorbifolium Bunge: Codification and description according to the BBCH scale. Sci. Hortic. 329, 113011 (2024).Article 
CAS 

Google Scholar 
Delgado, P. H. et al. Phenological growth stages of mango (Mangifera indica L.) according to the BBCH scale. Sci. Hortic. 130(3), 536–540 (2011).Article 

Google Scholar 
Liu, K. D. et al. Identification of phenological growth stages of sugar apple (Annona squamosa L.) using the extended BBCH-scale. Sci. Hortic. 181, 76–80 (2015).Article 

Google Scholar 
Fadón, E., Herrero, M. & Rodrigo, J. Flower development in sweet cherry framed in the BBCH scale. Sci. Hortic. 192, 141–147 (2015).Article 

Google Scholar 
García-Carbonell, S. et al. Phenological growth stages of the persimmon tree (Diospyros kaki). Ann. Appl. Biol. 141(1), 73–76 (2002).Article 

Google Scholar 
Atli, H. S., Arpaci, S., Tekin, H. & Yaman, A. Determination of the most suitable total temperature and harvest time of some pistachio cultivars[J]. Acta Hort. 470, 502–506 (1998).Article 

Google Scholar 
Odabaşıoğlu, M. İ. Comparison of various effective heat summation requirement (growing degree-day) calculation methods on different grape cultivars. Appl. Ecol. Environ. Res. 21(6), 5141–5162 (2023).Article 

Google Scholar 
Payero, J. Growing degree-day calculator for the southeast USA. Coop. Ext. AC10, 4 (2017).
Google Scholar 
Özbek, S. General fruit growing. Cukurova University Faculty of Agriculture Publications Adana 111(6), 386 (1977).Li, S. Q., Zhu, Y. P., Liu, H. L., Li, S. J. & Liu, S. P. Model construction and 3D visualisation of leaf height in winter wheat after regrowth. China Agric. Sci. Technol. Bull. 19(11), 59–67 (2017).
Google Scholar 
Zhao, G. C. et al. The phenological growth stages of Sapindus mukorossi according to BBCH scale. Forests 10(6), 462 (2019).Article 

Google Scholar 
Santos, F. R. et al. Phenological Growth Stages Based on the BBCH Scale and Thermal Requirements of Spondias dulcis Parkinson[J]. Appl. Fruit Sci. 67(4), 207–207 (2025).Article 

Google Scholar 
Amenan, J. K. et al. Phenological growth stages of shea tree (Vitellaria paradoxa subsp. paradoxa) according to the BBCH scale[J]. Ann. Appl. Biol. 182(1), 131–139 (2022).
Google Scholar 
Küsmüş, S. et al. Determination of phenological characters and effective heat summation requirements of local grape cultivars grown in Malatya province. J. Kırsehir Ahi Evran Univ Fac Agric 2(1), 9–23 (2022).
Google Scholar 
Tunç Yazgan, T., Yaman, M. & Yilmaz-Kadir-Uğurtan, K. U. Determination of phenotypic diversity and effective temperature sum times in some olive (Olea europaea L.) varieties by using phenological stages with multivariate analysis. Appl. Fruit Sci. 66(3), 1151–1161 (2024).Article 

Google Scholar 
Zou, L. et al. A study on the correlation between flowering lateness and meteorological factors in Magnoliaceae. Temp. For. Res. 7(03), 1–7 (2024).
Google Scholar 
Zamani-Noor, N. & Rodemann, B. Reducing the build-up of Plasmodiophora brassicae inoculum by early management of oilseed rape volunteers[J]. Plant. Pathol. 67(2), 426–432 (2018).Article 
CAS 

Google Scholar 
Li, X. J., Mu, Y. T., Zhang, X. M. & Wang, J. J. Phenological growth stages of Astragalus membranaceus var. Mongholicus according to the Biologische Bundesanstalt Bundessortenamt and chemical industry (BBCH) scale. Ann. Appl. Biol. 183(3), 302–319 (2023).Article 
CAS 

Google Scholar 
Wang, J. L. et al. Comparison of Picrasidine and Oxidised Picrasidine Contents in Flowers, Stems, Leaves and Seeds of Sophora moorcroftiana in Different Regions of Tibet at Different Harvesting Times. J. Tibet Univ. (Nat. Sci. Ed.) 27(01), 28–31 (2012).
Google Scholar 
Hogmire, H. W. & Biggs, A. R. Reduced pesticide programme for peach based on tree phenology. Crop Prot. 13(4), 277–280 (1994).Article 
CAS 

Google Scholar 
Download referencesFundingThis work was supported by Tibet Autonomous Region Science and Technology Program (XZ202401YD0026) and Phase I of the Forestry Doctoral Program at Xizang Agricultural and Animal Husbandry University (533325001).Author informationAuthor notesYanling Wan and Fumei Xin contributed equally to this work.Authors and AffiliationsXizang Agricultural and Animal Husbandry University, Nyingchi, 860000, People’s Republic of ChinaYanling Wan, Fumei Xin, Jiba Bianba, Xinlu Guo, Huanhuan Xie & Chenlong ZhangState Key Laboratory of Efficient Production of Forest Resources, Beijing Forestry University, 35 E Qinghua Rd., Beijing, 100083, People’s Republic of ChinaJiming LiuNational Energy R&D Center for Non-Food Biomass, Beijing Forestry University, Beijing, 100083, People’s Republic of ChinaJiming LiuSchool of Biological Sciences, Nanyang Technological University, Singapore, SingaporeJiming LiuAuthorsYanling WanView author publicationsSearch author on:PubMed Google ScholarFumei XinView author publicationsSearch author on:PubMed Google ScholarJiba BianbaView author publicationsSearch author on:PubMed Google ScholarXinlu GuoView author publicationsSearch author on:PubMed Google ScholarHuanhuan XieView author publicationsSearch author on:PubMed Google ScholarChenlong ZhangView author publicationsSearch author on:PubMed Google ScholarJiming LiuView author publicationsSearch author on:PubMed Google ScholarContributionsYanling Wan: Investigation, Writing—original draft, photography. Fumei Xin: Conceptualization, Validation, Visualization. Bianba Jiba: Methodology. Xinlu Guo: Investigation. Huanhuan Xie: photography. Chenlong Zhang: data curation. Jiming Liu: Validation, Supervision.Corresponding authorsCorrespondence to
Fumei Xin or Jiming Liu.Ethics declarations

Competing interests
The authors declare no competing interests.

Additional informationPublisher’s noteSpringer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.Supplementary InformationBelow is the link to the electronic supplementary material.Supplementary Material 1Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.
Reprints and permissionsAbout this articleCite this articleWan, Y., Xin, F., Bianba, J. et al. Establishing the phenological development stages of Sophora moorcroftiana using the BBCH scale.
Sci Rep 15, 44091 (2025). https://doi.org/10.1038/s41598-025-27760-wDownload citationReceived: 16 December 2024Accepted: 05 November 2025Published: 18 December 2025Version of record: 18 December 2025DOI: https://doi.org/10.1038/s41598-025-27760-wShare this articleAnyone you share the following link with will be able to read this content:Get shareable linkSorry, a shareable link is not currently available for this article.Copy shareable link to clipboard
Provided by the Springer Nature SharedIt content-sharing initiative
Keywords
Sophora
BBCH scaleDevelopment stagePhenologyAccumulated temperature More

AbstractBiodiversity loss threatens ecosystems and human well-being, making accurate, large-scale monitoring crucial. Environmental DNA (eDNA) has enabled species detection from substrates such as water, without the need for direct observation. Lately, airborne eDNA has been showing promise for tracking organisms from insects to mammals in terrestrial ecosystems. Conventional biodiversity assessments are often labor-intensive and limited in scope, leaving gaps in our understanding of ecosystem response to environmental change. Here, we demonstrate that airborne eDNA can detect organisms across the tree of life, quantify changes in abundance congruent with traditional monitoring, and reveal land-use induced regional decline of diversity in a northern boreal ecosystem over more than three decades. By analyzing 34 years of archived aerosol filters, we reconstruct weekly temporal relative abundance data for more than 2700 genera using non-targeted methods. This study provides unified, ecosystem-scale biodiversity surveillance spanning multiple decades, with data collected at weekly intervals on both the individual species and community level. Previously, large scale analyses of ecosystem changes, targeting all types of organisms, has been prohibitively expensive and difficult to attempt. Here, we present a way of holistically doing this type of analysis in a single framework.

Similar content being viewed by others

First national survey of terrestrial biodiversity using airborne eDNA

Article
Open access
02 June 2025

Shotgun sequencing of airborne eDNA achieves rapid assessment of whole biomes, population genetics and genomic variation

Article
Open access
03 June 2025

Archived natural DNA samplers reveal four decades of biodiversity change across the tree of life

Article
Open access
01 August 2025

Data availability

The sequencing data generated in this study are deposited in the NCBI Sequence Read Archive (SRA) under accession code PRJNA808200. The processed relative abundance data are available in Supplementary Data 6. External datasets used are land cover data (Swedish National Land Cover Database, www.naturvardsverket.se/en/services-and-permits/maps-and-map-services/national-land-cover-database/), map vector data (Natural Earth, www.naturalearthdata.com/), weather data (Copernicus Climate Change Service, https://doi.org/10.24381/cds.e2161bac, National Centers for Environmental Prediction (NCEP) and National Center for Atmospheric Research (NCAR) Reanalysis project, psl.noaa.gov/data/gridded/reanalysis/, Swedish Meteorological and Hydrological Institute (SMHI), www.smhi.se/data/hitta-data-for-en-plats/ladda-ner-vaderobservationer, Climatology Lab, www.climatologylab.org/terraclimate.html, National Oceanic and Atmospheric Administration (NOAA) – Climate Prediction Center, www.cpc.ncep.noaa.gov, Expert Team on Climate Change Detection and Indices (ETCCDI), etccdi.pacificclimate.org/data.shtml), reference sequence data (National Center for Biotechnology Information (NCBI), www.ncbi.nlm.nih.gov/nucleotide/, accession numbers for all sequences used in the Kraken database are available at https://doi.org/10.5281/zenodo.17778887), species observational data (Swedish Species Observation System database, artportalen.se, Global Biodiversity Information Facility (GBIF), www.gbif.org, Swedish Bird Survey, www.fageltaxering.lu.se, Sámi Parliament of Sweden (Sámediggi), sametinget.se/renstatistik), and forestry data (The Swedish National Forest Inventory (NFI), www.slu.se/en/about-slu/organisation/departments/forest-resource-management/miljoanalys/nfi/, Swedish Forest Agency, www.skogsstyrelsen.se/laddanergeodata).
Code availability

StringMeUp, a computer program developed in-house and used in the classification of the sequence data, and the Kraken 2 fork are both available under DOIs https://doi.org/10.5281/zenodo.17569636 and https://doi.org/10.5281/zenodo.17570001, respectively.
ReferencesNewbold, T. et al. Global effects of land use on local terrestrial biodiversity. Nature 520, 45–50 (2015).
Google Scholar 
Ceballos, G. et al. Accelerated modern human-induced species losses: Entering the sixth mass extinction. Sci. Adv. 1, e1400253 (2015).
Google Scholar 
Cristescu, M. E. & Hebert, P. D. N. Uses and misuses of environmental DNA in biodiversity science and conservation. Annu. Rev. Ecol. Evol. Syst. 49, 209–230 (2018).
Google Scholar 
Bálint, M. et al. Environmental DNA time series in ecology. Trends Ecol. Evol. 33, 945–957 (2018).
Google Scholar 
Seeber, P. A. & Epp, L. S. Environmental DNA and metagenomics of terrestrial mammals as keystone taxa of recent and past ecosystems. Mamm. Rev. 52, 538–553 (2022).
Google Scholar 
Djurhuus, A. et al. Environmental DNA reveals seasonal shifts and potential interactions in a marine community. Nat. Commun. 11, 254 (2020).
Google Scholar 
van der Heyde, M., Bunce, M. & Nevill, P. Key factors to consider in the use of environmental DNA metabarcoding to monitor terrestrial ecological restoration. Sci. Total Environ. 848, 157617 (2022).
Google Scholar 
Clare, E. L. et al. Measuring biodiversity from DNA in the air. Curr. Biol. 32, 693–700 (2022).
Google Scholar 
Lynggaard, C. et al. Airborne environmental DNA for terrestrial vertebrate community monitoring. Curr. Biol. 32, 701–707 (2022).
Google Scholar 
Littlefair, J. E. et al. Air-quality networks collect environmental DNA with the potential to measure biodiversity at continental scales. Curr. Biol. 33, R426–R428 (2023).
Google Scholar 
Després, V. R. et al. Primary biological aerosol particles in the atmosphere: A review. Tellus B Chem. Phys. Meteorol. 64, 15598 (2012).
Google Scholar 
Šantl-Temkiv, T., Amato, P., Casamayor, E. O., Lee, P. K. H. & Pointing, S. B. Microbial ecology of the atmosphere. FEMS Microbiol. Rev. 46, fuac009 (2022).
Google Scholar 
Fröhlich-Nowoisky, J. et al. Bioaerosols in the Earth system: Climate, health, and ecosystem interactions. Atmos. Res. 182, 346–376 (2016).
Google Scholar 
Métris, K. L. & Métris, J. Aircraft surveys for air eDNA: probing biodiversity in the sky. PeerJ. 11, e15171 (2023).
Google Scholar 
Karlsson, E. et al. Airborne microbial biodiversity and seasonality in Northern and Southern Sweden. PeerJ. 8, e8424 (2020).
Google Scholar 
Bowers, R. M. et al. Seasonal variability in bacterial and fungal diversity of the near-surface atmosphere. Environ. Sci. Technol. 47, 12097–12106 (2013).
Google Scholar 
Bowers, R. M., McLetchie, S., Knight, R. & Fierer, N. Spatial variability in airborne bacterial communities across land-use types and their relationship to the bacterial communities of potential source environments. ISME J. 5, 601–612 (2011).
Google Scholar 
Johnson, M. D., Cox, R. D., Grisham, B. A., Lucia, D. & Barnes, M. A. Airborne eDNA reflects human activity and seasonal changes on a landscape scale. Front. Environ. Sci. 8, 563431 (2021).
Google Scholar 
Johnson, M. D., Barnes, M. A., Garrett, N. R. & Clare, E. L. Answers blowing in the wind: Detection of birds, mammals, and amphibians with airborne environmental DNA in a natural environment over a yearlong survey. Environ. DNA 5, 375–387 (2023).
Google Scholar 
Lynggaard, C., Frøslev, T. G., Johnson, M. S., Olsen, M. T. & Bohmann, K. Airborne environmental DNA captures terrestrial vertebrate diversity in nature. Mol. Ecol. Resour. 24, e13840 (2024).
Google Scholar 
Roger, F. et al. Airborne environmental DNA metabarcoding for the monitoring of terrestrial insects—A proof of concept from the field. Environ. DNA 4, 790–807 (2022).
Google Scholar 
Pumkaeo, P., Takahashi, J. & Iwahashi, H. Detection and monitoring of insect traces in bioaerosols. PeerJ. 9, https://doi.org/10.7717/peerj.10862 (2021).Polling, M., Buij, R., Laros, I. & de Groot, G. A. Continuous daily sampling of airborne eDNA detects all vertebrate species identified by camera traps. Environ. DNA 6, e591 (2024).
Google Scholar 
Helin, A. et al. Characterization of free amino acids, bacteria and fungi in size-segregated atmospheric aerosols in boreal forest: Seasonal patterns, abundances and size distributions. Atmos. Chem. Phys. 17, 13089–13101 (2017).
Google Scholar 
Mamanova, L. et al. Target-enrichment strategies for next-generation sequencing. Nat. Methods 7, 111–118 (2010).
Google Scholar 
Gonzalez, A. et al. Avoiding pandemic fears in the subway and conquering the platypus. mSystems 1, e00050–16 (2016).
Google Scholar 
Lu, J. et al. Metagenome analysis using the Kraken software suite. Nat. Protoc. 17, 2815–2839 (2022).
Google Scholar 
Chen, T. & Guestrin, C. XGBoost: A scalable tree boosting system. In Proceedings of the ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, 785–794 (2016).GBIF.org GBIF Occurrence Download https://doi.org/10.15468/dl.cjxesu (2020).Yates, M. C., Fraser, D. J. & Derry, A. M. Meta-analysis supports further refinement of eDNA for monitoring aquatic species-specific abundance in nature. Environ. DNA 1, 5–13 (2019).
Google Scholar 
Yates, M. C. et al. The relationship between eDNA particle concentration and organism abundance in nature is strengthened by allometric scaling. Mol. Ecol. 30, 3068–3082 (2021).
Google Scholar 
Fediajevaite, J., Priestley, V., Arnold, R. & Savolainen, V. Meta-analysis shows that environmental DNA outperforms traditional surveys, but warrants better reporting standards. Ecol. Evol. 11, 4803–4815 (2021).Harrison, J. B., Sunday, J. M. & Rogers, S. M. Predicting the fate of eDNA in the environment and implications for studying biodiversity. Proc. R. Soc. B Biol. Sci. 286, 20191409 (2019).
Google Scholar 
Valentin, R. E. et al. Moving eDNA surveys onto land: Strategies for active eDNA aggregation to detect invasive forest insects. Mol. Ecol. Resour. 20, 746–755 (2020).
Google Scholar 
Kirtane, A., Kleyer, H. & Deiner, K. Sorting states of environmental DNA: Effects of isolation method and water matrix on the recovery of membrane-bound, dissolved, and adsorbed states of eDNA. Environ. DNA 5, 582–596 (2023).
Google Scholar 
Manninen, H. E. et al. Patterns in airborne pollen and other primary biological aerosol particles (PBAP), and their contribution to aerosol mass and number in a boreal forest. Boreal Environ. Res. 19, 383–405 (2014).
Google Scholar 
li, S. & Georgopoulos, P. A mechanistic modeling system for estimating large-scale emissions and transport of pollen and co-allergens. Atmos. Environ. 45, 2260–2276 (2011).
Google Scholar 
Nordén, J., Penttilä, R., Siitonen, J., Tomppo, E. & Ovaskainen, O. Specialist species of wood-inhabiting fungi struggle while generalists thrive in fragmented boreal forests. J. Ecol. 101, 701–712 (2013).Woo, C., An, C., Xu, S., Yi, S. M. & Yamamoto, N. Taxonomic diversity of fungi deposited from the atmosphere. ISME J. 12, 2051–2060 (2018).
Google Scholar 
Clauß, M. Particle size distribution of airborne micro-organisms in the environment-A review. Landbauforschung Volkenrode 65, 77–100 (2015).
Google Scholar 
Ruiz-Jimenez, J. et al. Determination of free amino acids, saccharides, and selected microbes in biogenic atmospheric aerosols – Seasonal variations, particle size distribution, chemical and microbial relations. Atmos. Chem. Phys. 21, 8775–8790 (2021).
Google Scholar 
Brook, J. R., Johnson, D. & Mamedov, A. Determination of the source areas contributing to regionally high warm season PM2.5 in eastern north america. J. Air Waste Manage Assoc. 54, 1162–1169 (2004).
Google Scholar 
Zhou, L., Hopke, P. K. & Liu, W. Comparison of two trajectory based models for locating particle sources for two rural New York sites. Atmos. Environ. 38, 1955–1963 (2004).Hopke, P. K. Review of receptor modeling methods for source apportionment. J. Air Waste Manage Assoc. 66, 237–259 (2016).
Google Scholar 
Belis, C. et al. European Guide on Air Pollution Apportionment with Receptor Models. (2019).Lavsund, S., Nygrén, T. & Solberg, E. J. Status of moose populations and challenges to moose management in Fennoscandia. Alces 39, 109–130 (2003).
Google Scholar 
Singh, N. J., Börger, L., Dettki, H., Bunnefeld, N. & Ericsson, G. From migration to nomadism: Movement variability in a northern ungulate across its latitudinal range. Ecol. Appl. 22, 2007–2020 (2012).
Google Scholar 
Watson, J. G., Chen, L. W. A., Chow, J. C., Doraiswamy, P. & Lowenthal, D. H. Source apportionment: Findings from the U.S. supersites program. J Air Waste Manage Assoc. 58, 265–288 (2008).
Google Scholar 
Blackman, R. et al. Environmental DNA: The next chapter. Mol. Ecol. 33, e17355 (2024).
Google Scholar 
Gloor, G. B., Macklaim, J. M., Pawlowsky-Glahn, V. & Egozcue, J. J. Microbiome datasets are compositional: And this is not optional. Front. Microbiol. 8, 2224 (2017).
Google Scholar 
Roche, K. E. & Mukherjee, S. The accuracy of absolute differential abundance analysis from relative count data. PLoS Comput. Biol. 18, e1010284 (2022).
Google Scholar 
Quinn, T. P., Richardson, M. F., Lovell, D. & Crowley, T. M. Propr: An R-package for identifying proportionally abundant features using compositional data analysis. Sci. Rep. 7, 16252 (2017).
Google Scholar 
Haas, J. C. et al. Microbial community response to growing season and plant nutrient optimisation in a boreal Norway spruce forest. Soil Biol. Biochem. 125, 197–209 (2018).
Google Scholar 
Bowers, R. M. et al. Sources of bacteria in outdoor air across cities in the midwestern United States. Appl. Environ. Microbiol. 77, 6350–6356 (2011).
Google Scholar 
van der Merwe, M., Ericson, L., Walker, J., Thrall, P. H. & Burdon, J. J. Evolutionary relationships among species of Puccinia and Uromyces (Pucciniaceae, Uredinales) inferred from partial protein coding gene phylogenies. Mycol Res. 111, 163–175 (2007).
Google Scholar 
Terhonen, E., Blumenstein, K., Kovalchuk, A. & Asiegbu, F. O. Forest tree microbiomes and associated fungal endophytes: Functional roles and impact on forest health. Forests 10, 42 (2019).
Google Scholar 
Ren, F. et al. Tissue microbiome of Norway spruce affected by heterobasidion-induced wood decay. Microb. Ecol. 77, 640–650 (2019).
Google Scholar 
Ross, A. A., Müller, K. M., Scott Weese, J. & Neufeld, J. D. Comprehensive skin microbiome analysis reveals the uniqueness of human skin and evidence for phylosymbiosis within the class Mammalia. Proc. Natl. Acad. Sci. USA 115, E5786–E5795 (2018).
Google Scholar 
Wiśniewska, K., Lewandowska, A. U. & Śliwińska-Wilczewska, S. The importance of cyanobacteria and microalgae present in aerosols to human health and the environment – Review study. Environ. Int. 131, 104964 (2019).
Google Scholar 
Vázquez, D. P., Gianoli, E., Morris, W. F. & Bozinovic, F. Ecological and evolutionary impacts of changing climatic variability. Biol. Rev. 92, 22–42 (2017).
Google Scholar 
Reeve, R. et al. How to partition diversity. Preprint at https://doi.org/10.48550/arXiv.1404.6520 (2016).Leinster, T. Entropy and Diversity: The Axiomatic Approach. (Cambridge University Press, Cambridge, 2021).Hill, M. O. Diversity and evenness: A unifying notation and its consequences. Ecology 54, 427–432 (1973).
Google Scholar 
Sax, D. F. & Gaines, S. D. Species diversity: From global decreases to local increases. Trends Ecol. Evol. 18, 561–566 (2003).
Google Scholar 
Clavel, J., Julliard, R. & Devictor, V. Worldwide decline of specialist species: Toward a global functional homogenization?. Front. Ecol. Environ. 9, 222–228 (2011).
Google Scholar 
Ylisirniö, A. L. et al. Dead wood and polypore diversity in natural post-fire succession forests and managed stands – Lessons for biodiversity management in boreal forests. For. Ecol. Manage 286, 16–27 (2012).
Google Scholar 
Uboni, A., Blochel, A., Kodnik, D. & Moen, J. Modelling occurrence and status of mat-forming lichens in boreal forests to assess the past and current quality of reindeer winter pastures. Ecol. Indic. 96, 99–106 (2019).
Google Scholar 
Jonsson, B. G. et al. Rapid changes in ground vegetation of mature boreal forests—an analysis of Swedish national forest inventory data. Forests 12, 475 (2021).
Google Scholar 
SLU Artdatabanken. Rödlistade Arter i Sverige 2020. (SLU, Uppsala, 2020).Sandström, J. et al. Impacts of dead wood manipulation on the biodiversity of temperate and boreal forests. A systematic review. J. Appl. Ecol. 56, 1770–1781 (2019).
Google Scholar 
Bergstedt, J., Hagner, M. & Milberg, P. Effects on vegetation composition of a modified forest harvesting and propagation method compared with clear-cutting, scarification and planting. Appl. Veg. Sci. 11, 159–168 (2008).
Google Scholar 
Edman, M., Gustafsson, M., Stenlid, J., Jonsson, B. G. & Ericson, L. Spore deposition of wood-decaying fungi: Importance of landscape composition. Ecography 27, 103–111 (2004).
Google Scholar 
Siitonen, P., Lehtinen, A. & Siitonen, M. Effects of forest edges on the distribution, abundance, and regional persistence of wood-rotting fungi. Conserv. Biol. 19, 250–260 (2005).
Google Scholar 
Lewin, H. A. et al. Earth BioGenome project: sequencing life for the future of life. Proc. Natl. Acad. Sci. USA 115, 4325–4333 (2018).
Google Scholar 
Masson, O. et al. Airborne concentrations and chemical considerations of radioactive ruthenium from an undeclared major nuclear release in 2017. Proc. Natl. Acad. Sci. USA 116, 16750–16759 (2019).
Google Scholar 
The Preparatory Commission for the Comprehensive Nuclear-Test-Ban Treaty Organization (CTBTO). Annu. Rep. 2022. (2023).Söderström, C., Ban, S., Jansson, P., Lindh, K. & Tooloutalaie, N. Radionuclides in Ground Level Air in Sweden Year 2006. (2007).Dabney, J. et al. Complete mitochondrial genome sequence of a Middle Pleistocene cave bear reconstructed from ultrashort DNA fragments. Proc. Natl. Acad. Sci. USA 110, 15758–15763 (2013).
Google Scholar 
Slon, V. et al. Neandertal and Denisovan DNA from Pleistocene sediments. Science 356, 605–608 (2017).
Google Scholar 
Martin, M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet J. 17, 10–12 (2011).
Google Scholar 
Bushnell, B.BBMap Short Read Aligner. Joint Genome Institute, Department of Energy (2014).Wood, D. E., Lu, J. & Langmead, B. Improved metagenomic analysis with Kraken 2. Genome Biol 20, 1–13 (2019).
Google Scholar 
Martín-Fernández, J. A., Hron, K., Templ, M., Filzmoser, P. & Palarea-Albaladejo, J. Bayesian-multiplicative treatment of count zeros in compositional data sets. Stat. Modelling 15, 134–158 (2015).
Google Scholar 
Palarea-Albaladejo, J. & Martín-Fernández, J. A. zCompositions — R package for multivariate imputation of left-censored data under a compositional approach. Chemometr. Intell. Lab. Syst. 143, 85–96 (2015).
Google Scholar 
Seabold, S. & Perktold, J. Statsmodels: econometric and statistical modeling with Python. InProceedings of the 9th Python in Science Conference 92–96 (2010).van den Boogaart, K. G. & Tolosana-Delgado, R. ‘compositions’: A unified R package to analyze compositional data. Comput. Geosci. 34, 320–338 (2008).
Google Scholar 
Oksanen, J. et al. vegan: Community Ecology Package. R package at https://CRAN.R-project.org/package=vegan (2020).GBIF.org GBIF Occurrence Download https://doi.org/10.15468/dl.xnyctg. (2020).Danecek, P. et al. Twelve years of SAMtools and BCFtools. Gigascience 10, giab008 (2021).
Google Scholar 
Camacho, C. et al. BLAST+: architecture and applications. BMC Bioinform. 10, 421 (2009).
Google Scholar 
Lindqvist, J. En Stokastisk Partikelmodell i Ett Icke-Metriskt Koordinatsystem. FOI-R–99-01086-862-SE, Swedish Defence Research Agency (1999).Muñoz-Sabater, J. ERA5-Land hourly data from 1950 to present. Copernicus Climate Change Service (C3S) Climate Data Store (CDS). https://doi.org/10.24381/cds.e2161bac (Accessed September 2019) (2019).Canty, A. & Ripley, B. boot: Bootstrap Functions (Originally by Angelo Canty for S). R package at https://cran.r-project.org/package=boot (2022).Davison, A. C. & Hinkley, D. V. Bootstrap Methods and Their Application. Bootstrap Methods and their Application (Cambridge University Press, Cambridge, 1997).Stein, A. F. et al. NOAA’s HYSPLIT atmospheric transport and dispersion modeling system. Bull. Am. Meteorol. Soc. 96, 2059–2077 (2015).
Google Scholar 
Kalnay, E. et al. The NCEP/NCAR 40-year reanalysis project. Bull. Am. Meteorol. Soc. 77, 437–472 (1996).
Google Scholar 
Carslaw, D. C. & Ropkins, K. Openair – An r package for air quality data analysis. Environ. Model. Softw. 27, 52–61 (2012).
Google Scholar 
Scott, S. L. & Varian, H. R. Predicting the present with Bayesian structural time series. Int. J. Math. Model. Num. Optimis. 5, 4–23 (2014).
Google Scholar 
Scott, S. L. bsts: Bayesian Structural Time Series. R package at https://CRAN.R-project.org/package=bsts (2022).Vehtari, A., Gelman, A. & Gabry, J. Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC. Stat. Comput. 27, 1413–1432 (2017).
Google Scholar 
Bürkner, P. C., Gabry, J. & Vehtari, A. Approximate leave-future-out cross-validation for Bayesian time series models. J. Stat. Comput. Simul. 90, 2499–2523 (2020).
Google Scholar 
Durbin, J. & Koopman, S. J. Time Series Analysis by State Space Methods. Time Series Analysis by State Space Methods (Oxford University Press, Oxford, 2012).Commandeur, J. J. F. & Koopman, S. J. An Introduction to State Space Time Series Analysis. (Oxford University Press, Incorporated, 2007).Geweke, J. Evaluating the Accuracy of Sampling-Based Approaches to the Calculation of Posterior Moments. in Bayesian Statistics (eds. Bernardo, J. M., Berger, O., Dawid, A. P. & Smith, A. F. M.) vol. 4 169–193 (Clarendon Press, Oxford, 1992).Raftery, A. E. & Lewis, S. M. Comment: One long run with diagnostics: Implementation strategies for markov chain monte carlo. Stat. Sci. 7, 493–497 (1992).
Google Scholar 
Plummer, M., Best, N., Cowles, K. & Vines, K. CODA: Convergence Diagnosis and Output Analysis for MCMC. R News 6, 7–11 (2006).
Google Scholar 
GBIF.org. GBIF Occurrence Download https://doi.org/10.15468/dl.k76kgd (2021).Holmes, E. E., Ward, E. J. & Wills, K. MARSS: Multivariate autoregressive state-space models for analyzing time-series data. R J 4, 11–19 (2012).
Google Scholar 
Holmes, E. E., Scheuerell, M. D. & Ward, E. J. Detecting a signal from noisy sensors. in Applied Time Series Analysis for Fisheries and Environmental Data. (2021).Download referencesAcknowledgementsWe thank Catharina Söderström and Johan Kastlander (CBRN Defense and Security, Swedish Defense Research Agency) for providing access to the air filter archive, and Benedicte Albrectsen and Göran Englund for their feedback on previous versions of this manuscript. We also wish to thank five anonymous reviewers for constructive criticism. We acknowledge support from the Science for Life Laboratory and the National Genomics Infrastructure (NGI) for providing assistance in massive parallel sequencing. The computations were enabled by resources provided by the National Academic Infrastructure for Supercomputing in Sweden (NAISS) and the Swedish National Infrastructure for Computing (SNIC) at UPPMAX and HPC2N, partially funded by the Swedish Research Council through grant agreement nos. 2022-06725 and 2018-05973. Thomas Ågren provided the organism illustrations in Figs. 1–3. Modified Copernicus Climate Change Service information 2020 was used for the catchment area analysis. Neither the European Commission nor the European Center for Medium-Range Weather Forecasts (ECMWF) is responsible for any use that may be made of the Copernicus information or data it contains. This study was supported by Formas (grant agreement nos. 2016-01371: PS, MF; 2019-00579: P.S., T.B., and M.F.; 2021-02155: PS, MF; 2024-01990: P.S., T.B., M.F., and N.S.), together with grants from Vetenskapsrådet (2021-06283: P.S. and M.F.), SciLifeLab Biodiversity fund (NP00048: P.S., M.F., and T.B.), Kempe foundation (JCK-1919: P.S., M.F., and T.B.), Umeå University Industrial research school (P.S.) and Swedish Defense Research Agency (M.F.)FundingOpen access funding provided by Umea University.Author informationAuthor notesThese authors contributed equally: Alexis R. Sullivan, Edvin Karlsson.Authors and AffiliationsDepartment of Ecology and Environmental Sciences, Umeå University, Umeå, SwedenAlexis R. Sullivan, Edvin Karlsson, Daniel Svensson, Jose Antonio Villegas, Daniel Bellieny, Abu Bakar Siddique, Per-Anders Esseen & Per StenbergDepartment of Wildlife, Fish and Environmental Studies, Swedish University of Agricultural Sciences, Umeå, SwedenAlexis R. Sullivan, Anita Norman, Navinder J. Singh & Tomas BrodinCBRN Defence and Security, Swedish Defence Research Agency (FOI), Umeå, SwedenEdvin Karlsson, Björn Brindefalk, Håkan Grahn, David Sundell, Andreas Sjödin, Mats Forsman & Per StenbergDepartment of Molecular Biosciences, The Wenner-Gren Institute, Stockholm University, Stockholm, SwedenBjörn BrindefalkUmeå Plant Science Centre, Department of Plant Physiology, Umeå University, Umeå, SwedenAmanda MikkoDepartment of Plant Biology, Swedish University of Agricultural Sciences, Uppsala, SwedenAbu Bakar SiddiqueDepartment of Molecular Biology, Umeå University, Umeå, SwedenAnna-Mia JohanssonAuthorsAlexis R. SullivanView author publicationsSearch author on:PubMed Google ScholarEdvin KarlssonView author publicationsSearch author on:PubMed Google ScholarDaniel SvenssonView author publicationsSearch author on:PubMed Google ScholarBjörn BrindefalkView author publicationsSearch author on:PubMed Google ScholarJose Antonio VillegasView author publicationsSearch author on:PubMed Google ScholarAmanda MikkoView author publicationsSearch author on:PubMed Google ScholarDaniel BellienyView author publicationsSearch author on:PubMed Google ScholarAbu Bakar SiddiqueView author publicationsSearch author on:PubMed Google ScholarAnna-Mia JohanssonView author publicationsSearch author on:PubMed Google ScholarHåkan GrahnView author publicationsSearch author on:PubMed Google ScholarDavid SundellView author publicationsSearch author on:PubMed Google ScholarAnita NormanView author publicationsSearch author on:PubMed Google ScholarPer-Anders EsseenView author publicationsSearch author on:PubMed Google ScholarAndreas SjödinView author publicationsSearch author on:PubMed Google ScholarNavinder J. SinghView author publicationsSearch author on:PubMed Google ScholarTomas BrodinView author publicationsSearch author on:PubMed Google ScholarMats ForsmanView author publicationsSearch author on:PubMed Google ScholarPer StenbergView author publicationsSearch author on:PubMed Google ScholarContributionsP.S., M.F., T.B., and E.K. conceived and designed the study; E.K. and A.M.J. extracted DNA; DSv constructed the database and performed read classification; E.K., A.R.S., D.B., D.S.v. pre-processed the data; A.R.S. designed and implemented the machine learning approach; and H.G. constructed the particle models. A.R.S. and E.K. conducted most of the data analysis, with support from D.S.v., D.B., A.B.S., J.A.V., A.M., D.S.u., B.B., A.N., A.S., N.S., and P.A.E. E.K., A.R.S., D.S.v., B.B., P.S., and N.S. wrote the first draft of the manuscript. All authors contributed intellectual input and approved the final versionCorresponding authorCorrespondence to
Per Stenberg.Ethics declarations

Competing interests
The Authors declare no competing interests.

Peer review

Peer review information
Nature Communications thanks David Schmale III and the other anonymous reviewer(s) for their contribution to the peer review of this work. A peer review file is available.

Additional informationPublisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.Supplementary informationSupplementary informationDescriptions of Additional Supplementary FilesSupplementary Data 1Supplementary Data 2Supplementary Data 3Supplementary Data 4Supplementary Data 5Supplementary Data 6Supplementary Data 7Supplementary Data 8Supplementary Data 9Supplementary Data 10Supplementary Data 11Supplementary Data 12Supplementary Data 13Reporting SummaryTransparent Peer Review fileRights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
Reprints and permissionsAbout this articleCite this articleSullivan, A.R., Karlsson, E., Svensson, D. et al. Airborne eDNA captures three decades of ecosystem biodiversity.
Nat Commun (2025). https://doi.org/10.1038/s41467-025-67676-7Download citationReceived: 17 April 2025Accepted: 05 December 2025Published: 18 December 2025DOI: https://doi.org/10.1038/s41467-025-67676-7Share this articleAnyone you share the following link with will be able to read this content:Get shareable linkSorry, a shareable link is not currently available for this article.Copy shareable link to clipboard
Provided by the Springer Nature SharedIt content-sharing initiative More

AbstractVocalizations play crucial roles in dolphin biological activities. Analysis of dolphin vocalizations provides valuable insights into their behaviors and population status. In this data descriptor, we present a dataset of underwater vocalizations of Indo-Pacific humpback dolphins (Sousa chinensis) recorded in Xiamen Bay, China. The dataset comprises a diverse range of dolphin emissions, including 143 whistles (100 of which were classified as high-quality), and 897 pulse trains, categorized as echolocation clicks, burst pulses, and buzzes. A range of acoustic parameters was measured to characterize these acoustic signals. The presented dataset serves as an essential contribution to addressing existing data gaps regarding vocalizations of the Indo-Pacific humpback dolphin population in Xiamen Bay. It provides an important resource for studying vocalization patterns and temporal variability in the acoustic behaviors of the Indo-Pacific humpback dolphins, offering key insights to inform conservation strategies for this endangered population. Additionally, the dataset holds potential for population connectivity research, enabling acoustic comparisons between dolphin populations across different geographic regions to assess potential isolation or interaction.

Similar content being viewed by others

Vocal universals and geographic variations in the acoustic repertoire of the common bottlenose dolphin

Article
Open access
04 June 2021

A WAV file dataset of bottlenose dolphin whistles, clicks, and pulse sounds during trawling interactions

Article
Open access
22 September 2023

Comparing visual and acoustic detectability of two coastal cetacean species off Sindhudurg, India, to better inform integrated survey protocol

Article
Open access
24 August 2025

Background & SummaryTo adapt to the complex and dark underwater environment, dolphins have evolved the ability to use sound for sensing their surroundings, navigation, foraging, and communication1. Dolphin vocalizations primarily consist of two main types: pulsed signals and frequency-modulated whistles. Pulsed signals, including echolocation clicks, burst pulses, and buzzes, are characterized by high-frequency, broadband properties and are typically produced in trains2. Among these, echolocation clicks are the most commonly observed acoustic signals. Dolphins emit highly directional echolocation clicks while navigating and detecting prey or other targets of interest2. Burst pulses share similar characteristics with echolocation clicks but generally have lower frequencies and higher repetition rates. These signals are primarily used for intraspecific communication3. Buzzes are associated with close-range echolocation, particularly during the capture phase of hunting, as indicated by a rapid increase in pulse repetition rate. In addition to their role in predation, buzzes also serve as functional signals in social interactions, such as mating behaviors and mother-calf communication4. Whistles, which are narrow-band signals with modulated frequencies, are predominantly used for communication1. This type of emission plays crucial roles in social contexts, including reproductive gathering5, group cohesion6, individual identification7, coordinating group activities8.Dolphin vocalizations, both whistles and pulsed signals, are highly dynamic and flexible. The characteristics of dolphin sound production are known to vary depending on the environmental conditions9,10,11,12,13,14,15. For instance, Jensen et al.9 demonstrated that Irrawaddy dolphins (Orcaella brevirostris) and Ganges river dolphins (Platanista gangetica gangetica) produce echolocation clicks with significantly lower peak frequencies and sound amplitudes in shallow waters compared to deep waters, likely to mitigate reverberation. The minimum, maximum, and peak frequencies of whistles produced by common bottlenose dolphins (Tursiops truncatus) were found to increase with the ambient noise10. Changes in both click and whistle parameters in response to vessel presence have also been reported11,12,13,14. Additionally, alterations in dolphin sound production patterns have been linked to factors such as group size, group composition, and behavioral contexts12,15,16,17. For example, Akiyama and Ohta18 observed that common bottlenose dolphins increase whistle production during feeding, with a preference for upsweep-contour whistles. In dolphin groups containing calves, whistles tend to have lower end and minimum frequencies and longer durations12. Dolphins also exhibit adaptive modifications in their click characteristics during prey perception and target detection. They dynamically adjust their click rate, output amplitude, and acoustic directivity across different phases of target approach to achieve precise detection and recognition19,20,21,22. Recently, research on the adaptive acoustic control of dolphin biosonar during target detection has intensified, driven by the dual objectives of understanding biological systems and inspiring biomimetic applications.In summary, dolphin vocalizations are complex and dynamic yet critical for their survival. The collection and analysis of these acoustic signals can help to better understand their behaviors and population status. As one of the vulnerable dolphin species assessed by the IUCN23, the Indo-Pacific humpback dolphin (Sousa chinensis) in Chinese waters is sporadically distributed along the southeastern coast of China, including Xiamen Bay. The population inhabiting Xiamen Bay is under significant survival pressure due to increasing anthropogenic disturbances, with a notable reduction in population size over the past two decadesSousa chinensis) in Chinese Waters. Aquatic Mammals 30, 149–158 (2004).” href=”https://www.nature.com/articles/s41597-025-06253-5#ref-CR24″ id=”ref-link-section-d279053545e574″>24,25,26,27,28,29,30. Enhanced conservation efforts are urgently needed to mitigate ongoing threats to this endangered population. A comprehensive assessment of population status, including population size and structure, distribution patterns, surface behaviors, and vocalizations, is essential for developing effective conservation strategies and management initiatives. While studies have reported on vocalizations of Indo-Pacific humpback dolphins in Chinese waters31,32,33,34,35,36,37,38,39, limited research focused on the Xiamen Bay population, necessitating more targeted bioacoustic investigations.In this paper, we present a dataset of acoustic recordings of Indo-Pacific humpback dolphin vocalizations, collected during our regular surveys of the species’ population status. These data represent an important supplement to the Indo-Pacific humpback dolphin sound library in Chinese waters and, to a large extent, fill the data gaps regarding the dolphin population in Xiamen Bay. In addition to sharing the complete acoustic data, a quantitative analysis was conducted on dolphin vocalizations, both whistles and pulsed signals (classified as echolocation clicks, burst pulses, and buzzes), to determine characteristic parameters of each signal, and the results are incorporated into the dataset. It offers a valuable resource to study the vocalization patterns of the Indo-Pacific humpback dolphin population in Xiamen Bay, as well as to investigate population connectivity through acoustic comparisons across different geographic regions, helping to assess potential population isolation or interaction. Notably, this dataset spans a period of up to three years, allowing for investigations into potential temporal variability in dolphin acoustic behaviors over a long timeline. The Indo-Pacific humpback dolphins inhabit shallow coastal waters where they face significant echolocation challenges due to high reverberation. Dolphins present remarkable biosonar plasticity, adjusting their acoustic signals in response to varying environmental conditions9,10,11,12,13,14,15. This dataset contains echolocation click trains recorded across a wide range of water depths from 2.1 m to 24.3 m, providing valuable acoustic materials for exploring biosonar operational mechanisms— particularly how dolphins adapt their signals in shallow-water environments with strong reverberation. The comprehensive collection of dolphin vocalization signals also serves as an important resource for bionic applications40,41,42,43.MethodsData collectionWe conducted a total of 22 regular vessel-based surveys in Xiamen Bay, China, from June 2022 to September 2024, to investigate the abundance and distribution patterns of the Indo-Pacific humpback dolphin population, as well as their underwater vocalizations (Fig. 1). During the surveys, the vessel navigated at a speed of 5–7 knots, with three observers equipped with Navigator 7 × 50 binoculars (magnification: 7×, objective lens diameter: 50 mm, field of view: 419 ft at 1000 meters, Steiner company, Germany) positioned on the forward deck to search for dolphin presence. Upon sighting dolphins, the vessel approached at a reduced speed and stopped its engine once a distance was reached less than 50 m. During this period, a calibrated underwater acoustic recorder, SoundTrap 300HF (Ocean Instruments, Auckland, New Zealand), was deployed to record the underwater sounds produced by the dolphins. The recorder was equipped with a 16-bit analog-to-digital converter (ADC) and a pre-amplified hydrophone exhibiting a linear frequency response across the frequency range from 20 Hz to 150 kHz, with a sensitivity of −189 ± 3 dB (re 1 V/μPa) in low-gain mode. The recorder was housed in a steel holder and positioned 1.5–2 m underwater using a steel pipe. A sampling rate of 576 kS/s was used, and the recorded sound data were stored as WAV audio files in real time. During acoustic recording, the vessel engine was turned off to minimize low-frequency sound interference. Dolphins were photographed for individual identification, and detailed information about each acoustic recording session was documented, including geographic location, water depth, dolphin group size, and dolphin behavioral state. From our field surveys, a total of 1019 minutes of original sound data were collected. These recordings were subsequently processed to extract communication whistles and high-frequency pulsed signals produced by the Indo-Pacific humpback dolphin.Fig. 1(a) Map of the survey area, with red dots indicating locations where acoustic recordings of dolphin vocalizations were conducted. (b) Aerial photographs of the survey vessel and an Indo-Pacific humpback dolphin captured by an unmanned aerial vehicle. (c) The SoundTrap 300HF underwater acoustic logger used in the survey.Full size imageDetection of whistlesThe original WAV files were displayed in the time-frequency domain using Adobe Audition (Version 2021) and manually reviewed by an experienced observer to identify whistles within consecutive 3-second time windows. A continuous tonal contour without temporal breakpoints on the spectrogram was identified as a single whistle. Additionally, consecutive contours were also considered a single whistle if the gap between them was shorter than 200 ms and less than the duration of the contours, following established methodologies35,44. Each identified whistle signal was subsequently extracted and saved as an individual WAV file, and the number and position of whistles within the original acoustic file were documented. Whistles were visually divided into three grades based on signal-to-noise ratio (SNR), referencing previous studies12,45. Grade 1 includes whistles with weak contour, but is visible in the spectrogram; Grade 2 includes whistles with clear and unambiguous contours; and Grade 3 includes whistles with contour prominent in the spectrogram. Whistles of Grade 2 and 3 were considered to be of high quality and selected for detailed analysis. All whistle contours were visually categorized into six tonal types32,35,36,44 (Fig. 2): (a) constant, (b) upsweep, (c) downsweep, (d) concave, (e) convex, and (f) sinusoidal. For each high-quality whistle, 13 acoustic characteristics were manually measured from the spectrograms: (a) duration (ms), (b) start frequency (kHz), (c) end frequency (kHz), (d) frequency change (kHz), (e) absolute frequency gradient (kHz/s), (f) minimum frequency (kHz), (g) maximum frequency (kHz), (h) delta frequency (kHz), (i) number of extrema, (j) number of inflection points, (k) number of saddle points, (l) number of breaks, and (m) presence/absence of harmonics. Detailed definitions of these acoustic parameters were consistent with those described in Marley et al.45.Fig. 2Example frequency contours illustrating the six whistle types: (a) constant, (b) upsweep, (c) downsweep, (d) concave, (e) convex, (f) sinusoidal.Full size imageDetection of pulsed signalsThe raw acoustic data contained various pulsed sounds produced by the Indo-Pacific humpback dolphins, including echolocation clicks, burst pulses, and buzzes. These pulsed signals have been shown to exhibit a Gabor-type waveform structure46,47, characterized by a distinct Gaussian envelope in the output of the Teager–Kaiser energy operator (TKEO). Based on this feature, a mature processing method proposed by Madhusudhana et al.48 was applied to extract pulsed signals from the recordings. Initially, raw WAV files were split into 30-second segments to reduce the computational load. A Butterworth high-pass filter with a cut-off frequency of 5 kHz was used to remove most low-frequency noise from the acoustic data.Subsequently, the TKEO output was calculated for each filter data segment as follows49:$${Psi }_{d}[{x}_{n}]={x}_{n}^{2}-{x}_{n-1}{x}_{n+1}$$
(1)
where xn represents the sampled points of each 30-second signal segment. A Gaussian-weighted averaging filter (MAF1) was then applied to highlight short-duration energy surges in TKEO outputs. The Gaussian filter was defined as:$$MAF1(n)=frac{{T}_{s}}{{sigma }_{G}sqrt{2pi }}{e}^{-{(n{T}_{s})}^{2}/2{{sigma }_{G}}^{2}}$$
(2)
where n = −N, …, 0, …, N, is the index of the sampled point in the filter, Ts is the sampling interval, and ({sigma }_{G}) is the standard deviation of the Gaussian, given by:$${sigma }_{G}=frac{FWHM}{2sqrt{2,mathrm{ln}(2)}}$$
(3)
where FWHM is the width of the Gaussian at half its peak value, which is taken as 1 × 10−4 s, based on the length of a representative pulsed signal produced by the dolphin.The value of N is determined based on ({sigma }_{G}), given as the following formula:$$N=frac{5{sigma }_{G}}{{T}_{s}}$$
(4)
where Ts is the sampling interval.Convolution of TKEO output with MAF1 can be expressed as:$${h}_{MAF1}(n)=frac{{T}_{s}}{{sigma }_{G}sqrt{2pi }}mathop{sum }limits_{i=-N}^{N}{e}^{-{(i{T}_{s})}^{2}/2{{sigma }_{G}}^{2}}{x}_{n+i}$$
(5)
TKEO outputs were also filtered by a rectangular averaging filter (MAF2) with identical length and filter gain to MAF1:$$MAF2(n)=frac{mathop{sum }limits_{m=-N}^{N}MAF1(m)}{2N+1}$$
(6)
The filter difference ratio (FDR) was then computed to measure the difference in the responses of the two filters:$$FDR(n)=frac{{h}_{MAF1}(n)-{h}_{MAF2}(n)}{{h}_{MAF1}(n)}$$
(7)
The obtained FDR is expected to yield local maximums at locations of Gaussian-like spikes in the TKEO outputs, indicating the presence of pulsed signals in the original recordings (Fig. 3). Since the FDR values remains relatively constant irrespective of signal amplitude, a detection threshold was set at 85% of the peak FDR value (FDRpeak) following Madhusudhana et al.48 Pulsed signals were confirmed only when the FDR value exceeded this threshold and the signal-to-noise ratio (SNR) of the signal was above 10 dB. For all identified pulses, the inter-pulse intervals (IPIs) were measured, which is defined as the time interval between two consecutive pulses.Fig. 3(a) Example data segment of dolphin vocalizations containing clicks, buzzes, and burst pulses. (b) A representative single pulsed signal, along with its outputs from (c) the Teager–Kaiser energy operator (TEKO) and (d) the filter difference ratio (FDR). (e) Positions of dolphin pulsed signals within the data segment as identified by the automated processing procedure.Full size imagePulse trains were identified by using an adaptive IPI threshold50. Consecutive pulses with a gradual change in the IPI were recognized to belong to the same pulse train, and a pulse train was considered terminated when an abrupt IPI increase occurred. Each pulse train detected was visually confirmed, and those containing signals with high reverberation and from more than one animal were removed. Based on IPI thresholds established by Wang et al.37, the identified pulse trains were subsequently classified into echolocation clicks, burst pulses, and buzzes, according to their mean IPI values. Pulse trains with a mean IPI value less than 4.9 ms were recognized as buzzes, those with mean IPIs greater than 15.5 ms were considered echolocation clicks, and pulse trains with intermediate mean IPIs (4.9–15.5 ms) were categorized as burst pulses. Each identified pulse train, along with its constituent individual pulses, were saved as separate WAV files within the dataset. Pulse trains were sequentially numbered based on their chronological order of appearance in the original recordings, and the position and total number of pulse trains in each original recording were documented. For each identified pulse train, six acoustic parameters were measured for the pulsed signals according to previous studies2,51,52: (a) inter-pulse-interval (IPI), (b) sound pressure level (SPLpp), (c) duration, (d) peak frequency (Fpeak), (e) -3dB bandwidth (-3dBBW), and (f) -10dB bandwidth (-10dBBW). Additionally, the mean values of these acoustic parameters were calculated and recorded for each pulse train.Data RecordsThe acoustic dataset contained a total of 143 whistle signals, of which 100 were identified as high-quality based on visual inspection. Additionally, 897 pulse trains were included, comprising 832 echolocation click trains, 15 burst pulse trains, and 50 buzz trains. The complete dataset is publicly accessible via an unrestricted repository at figshare53, consisting of WAV audio files, TXT text document files, Excel data sheets, and Portable Network Graphic (PNG) image files.The acoustic recordings are provided in different folders based on the signal type:

1.

Original Audio File contains the complete original recordings collected from field surveys, including 35 WAV audio files sequentially named according to the recording time (e.g., Ori_Recording_01.wav, Ori_Recording_02.wav, Ori_Recording_03.wav).

2.

Whistles comprises 143 extracted whistle signals saved as separate WAV files (e.g., Whistle_001.wav, Whistle_002.wav, Whistle_003.wav), along with PNG image files illustrating spectrograms of each whistle (e.g., Whistle_001.png, Whistle_002. png, Whistle_003.png).

3.

Click Trains, Burst Pulse Trains, and Buzz Trains contain WAV audio files corresponding to the detected click trains, burst pulse trains, and buzz trains, respectively. Each pulse train, along with its individual pulsed signals, is saved in a separate subfolder. These subfolders are sequentially named according to the serial number of the pulse train (e.g., PulseTrain_001, PulseTrain_002, PulseTrain_003). Each subfolder additionally contains a tab-separated TXT file named “PulseParameters.txt”, documenting acoustic parameters for each pulsed signal.

Results.xlsx is an Excel file containing detailed information on the original acoustic recordings and quantitative data on whistles and pulsed signals detected within each original acoustic file:

1.

recording dates when the original acoustic data were collected;

2.

geographic coordinates (latitude and longitude) of the recording locations;

3.

start and end times of each original acoustic file;

4.

number of pulse trains detected within each original acoustic file;

5.

number of identified whistles, including counts of high-quality whistles (Grade 2 and 3) within each original acoustic file;

6.

additional information about the recordings, including season, tide, water depth, dolphin group size, and behavioral state of dolphins.

Whistles.xlsx is an Excel file that provides detailed descriptions of each whistle signal:

1.

the original acoustic file from which the whistle signal was extracted;

2.

start and end times of each whistle signal within the original acoustic file, relative to the file’s start, expressed in seconds;

3.

contour type classification for each whistle signal (constant, upsweep, downsweep, concave, convex, and sinusoidal);

4.

quality grading for each whistle (Grade 1, Grade 2, and Grade 3);

5.

acoustic parameters of each whistle signal: Duration, Start Frequency, End Frequency, Frequency Change, Absolute Frequency Gradient, Minimum Frequency, Maximum Frequency, Delta Frequency, Number of Extrema, Number of Inflection Points, Number of Saddle Points, Number of Breaks, and Presence/Absence of Harmonics.

ClickTrains.xlsx, BurstPulseTrains.xlsx, and BuzzTrains.xlsx are Excel files containing detailed information on each identified pulse train type:

1.

serial number of pulse trains denoting their order of appearance within original recordings;

2.

original acoustic file from which each pulse train was extracted;

3.

start and end times of each pulse train within the original acoustic file, relative to the file’s start, expressed in seconds;

4.

lengths of each pulse train (Length);

5.

number of pulsed signals within each pulse train (NumP);

6.

mean values of acoustic parameters of the pulsed signals within each pulse train: mean SNR (Mean_SNR), mean IPI (Mean_IPI), mean sound pressure level (Mean_SPLpp), mean duration (Mean_Duration), mean peak frequency (Mean_Fpeak), mean -3dB bandwidth (Mean_-3dBBW), and mean -10dB bandwidth (Mean_-10dB BW).

Technical ValidationTo ensure the reliability of the presented dataset, the signal processing procedures in this study strictly adhered to the methodologies reported in published studies32,35,36,44,45,48,54,55. The identification, quality grading, classification, and subsequent characteristic measurements of the whistles were manually completed and reviewed by trained and experienced operators using the specialized acoustic signal processing software (Adobe Audition, Version 2021). Statistical results for the whistle parameters are summarized in Table 1, including the range (minimum to maximum), mean, and standard deviation for the 12 measured characteristic parameters of the 100 high-quality whistles.Table 1 Summary statistics of the characteristic parameters for the 100 high-quality whistles.Full size tableLocal extrema, inflections, saddles, breaks, and harmonics were observed in 56%, 36%, 30%, 42%, and 65% of all high-quality whistles, respectively. The constant-shaped contour was the most dominant tonal type, accounting for 48.25% of all whistles identified in this dataset (Fig. 4), consistent with previous reported data for Indo-Pacific humpback dolphin populations in Zhanjiang, Sanniang Bay, and Hainan waters32,35,36. This consistency provides support for the reliability of the presented dataset.Fig. 4Proportions of the six contour types in all whistles in the dataset.Full size imageDetection and classification of pulsed signals were conducted following standard procedures and methodologies widely adopted in published research54,55,56,57,58,59,60. Clicks were the most frequent pulse type, comprising approximately 93% of the identified 897 pulse trains, with a total of 832 occurrences. The acoustic characteristics of the three pulse train types are statistically summarized in Table 2, and each parameter was averaged per train. Consistent with previous findings37, the pulsed signals in this dataset were characterized by high-frequency, broadband features, with click trains exhibiting longer train lengths and higher peak frequencies compared to burst pulses and buzzes. Among the three sound types, burst pulses had the lowest peak frequency and bandwidth. These results further support the reliability of the presented dataset.Table 2 Descriptive statistics for acoustic characteristics of pulse trains produced by Indo-Pacific humpback dolphins in Xiamen Bay, China.Full size table

Data availability

The acoustic dataset described in the current paper is publicly accessible via an unrestricted repository at figshare (https://doi.org/10.6084/m9.figshare.29143727).
Code availability

Data processing followed established published methodologies, and no custom code was used in this paper. For detailed procedures on pulsed signal detections, refer to Madhusudhana et al.44.
ReferencesAu, W. W. L. & Hastings, M. C. Principles of Marine Bioacoustics. (Springer, New York, 2008).Au, W. W. L. The Sonar of Dolphins. (Springer Science & Business Media, 1993).Yoshida, Y. M. et al. Sound variation and function in captive Commerson’s dolphins (Cephalorhynchus commersonii). Behavioural Processes 108, 11–19 (2014).Article 
PubMed 

Google Scholar 
Martin, M. J., Elwen, S. H., Kassanjee, R. & Gridley, T. To buzz or burst-pulse? The functional role of Heaviside’s dolphin, Cephalorhynchus heavisidii, rapidly pulsed signals. Animal Behaviour 150, 273–284 (2019).Article 

Google Scholar 
Belikov, R. A. & Bel’kovich, V. M. Whistles of beluga whales in the reproductive gathering off Solovetskii Island in the White Sea. Acoust. Phys. 53, 528–534 (2007).Article 
ADS 

Google Scholar 
King, S. L. & Janik, V. M. Bottlenose dolphins can use learned vocal labels to address each other. Proc. Natl. Acad. Sci. USA 110, 13216–13221 (2013).Article 
ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
Janik, V. M. & Sayigh, L. S. Communication in bottlenose dolphins: 50 years of signature whistle research. J Comp Physiol A 199, 479–489 (2013).Article 

Google Scholar 
King, S. L., Guarino, E., Donegan, K., McMullen, C. & Jaakkola, K. Evidence that bottlenose dolphins can communicate with vocal signals to solve a cooperative task. R. Soc. open sci. 8, rsos.202073, 202073 (2021).Jensen, F. H. et al. Clicking in Shallow Rivers: Short-Range Echolocation of Irrawaddy and Ganges River Dolphins in a Shallow, Acoustically Complex Habitat. PLoS ONE 8, e59284 (2013).Article 
ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
van Ginkel, C., Becker, D. M., Gowans, S. & Simard, P. Whistling in a noisy ocean: bottlenose dolphins adjust whistle frequencies in response to real-time ambient noise levels. Bioacoustics 27, 391–405 (2018).Article 

Google Scholar 
Luís, A. R., Couchinho, M. N. & Dos Santos, M. E. Changes in the acoustic behavior of resident bottlenose dolphins near operating vessels. Mar Mam Sci 30, 1417–1426 (2014).Article 

Google Scholar 
Heiler, J., Elwen, S. H., Kriesell, H. J. & Gridley, T. Changes in bottlenose dolphin whistle parameters related to vessel presence, surface behaviour and group composition. Animal Behaviour 117, 167–177 (2016).Article 

Google Scholar 
Rako Gospić, N. & Picciulin, M. Changes in whistle structure of resident bottlenose dolphins in relation to underwater noise and boat traffic. Marine Pollution Bulletin 105, 193–198 (2016).Article 

Google Scholar 
Antichi, S., Urbán, R. J., Martínez-Aguilar, S. & Viloria-Gómora, L. Changes in whistle parameters of two common bottlenose dolphin ecotypes as a result of the physical presence of the research vessel. PeerJ 10, e14074 (2022).Article 
PubMed 
PubMed Central 

Google Scholar 
Quick, N. J. & Janik, V. M. Whistle rates of wild bottlenose dolphins (Tursiops truncatus): Influences of group size and behavior. Journal of Comparative Psychology 122, 305–311 (2008).Article 
PubMed 

Google Scholar 
Hernandez, E. N., Solangi, M. & Kuczaj, S. A. II. Time and frequency parameters of bottlenose dolphin whistles as predictors of surface behavior in the Mississippi Sound. The Journal of the Acoustical Society of America 127, 3232–3238 (2010).Article 
ADS 
PubMed 

Google Scholar 
Díaz López, B. Whistle characteristics in free-ranging bottlenose dolphins (Tursiops truncatus) in the Mediterranean Sea: Influence of behaviour. Mamm Biol 76, 180–189 (2011).Article 

Google Scholar 
Akiyama, J. & Ohta, M. Increased Number of Whistles of Bottlenose Dolphins, Tursiops truncatus, Arising from Interaction with People. J. Vet. Med. Sci. 69, 165–170 (2007).Article 
PubMed 

Google Scholar 
Jensen, F. H., Bejder, L., Wahlberg, M. & Madsen, P. T. Biosonar adjustments to target range of echolocating bottlenose dolphins(Tursiops sp.) in the wild. Journal of Experimental Biology 212, 1078–1086 (2009).Article 
CAS 
PubMed 

Google Scholar 
Kloepper, L. N. et al. Cognitive Adaptation of Sonar Gain Control in the Bottlenose Dolphin. PLOS ONE 9, e105938 (2014).Article 
ADS 
PubMed 
PubMed Central 

Google Scholar 
Ladegaard, M., Jensen, F. H., Beedholm, K., Da Silva, V. M. F. & Madsen, P. T. Amazon river dolphins (Inia geoffrensis) modify biosonar output level and directivity during prey interception in the wild. Journal of Experimental Biology jeb.159913, https://doi.org/10.1242/jeb.159913 (2017).Ladegaard, M. et al. Dolphin echolocation behaviour during active long-range target approaches. Journal of Experimental Biology jeb.189217, https://doi.org/10.1242/jeb.189217 (2018).Jefferson, T. A. & Smith, B. D. Re-assessment of the Conservation Status of the Indo-Pacific Humpback Dolphin (Sousa chinensis) Using the IUCN Red List Criteria. in Advances in Marine Biology (eds. Jefferson, T. A. & Curry, B. E.) vol. 73 1–26 (Academic Press, 2016).Jefferson, T. A. & Hung, S. K. A Review of the Status of the Indo-Pacific Humpback Dolphin (<I>Sousa chinensis</I>) in Chinese Waters. Aquatic Mammals 30, 149–158 (2004).Article 

Google Scholar 
Chen, B. et al. Abundance, distribution and conservation of Chinese White Dolphins (Sousa chinensis) in Xiamen, China. Mammalian Biology 73, 156–164 (2008).Article 

Google Scholar 
Chen, B., Zheng, D., Yang, G., Xu, X. & Zhou, K. Distribution and conservation of the Indo–Pacific humpback dolphin in China. Integrative Zoology 4, 240–247 (2009).Article 
PubMed 

Google Scholar 
Chen, B. et al. Survival rate and population size of Indo-Pacific humpback dolphins (Sousa chinensis) in Xiamen Bay, China. Marine Mammal Science 34, 1018–1033 (2018).Article 
ADS 

Google Scholar 
Wu, F., Wang, X., Ding, X., Miao, X. & Zhu, Q. Distribution Pattern of Indo-Pacific Humpback Dolphins (Sousa chinensis) along Coastal Waters of Fujian Province, China. Aquatic Mammals 40, 341–349 (2014).Article 

Google Scholar 
Wang, X. et al. External Injuries of Indo-Pacific Humpback Dolphins (Sousa chinensis) in Xiamen, China, and Its Adjacent Waters as an Indicator of Potential Fishery Interactions. Aquatic Mammals 44, 285–292 (2018).Article 

Google Scholar 
Zeng, Q. et al. Modeling demographic parameters of an edge-of-range population of Indo-Pacific humpback dolphin in Xiamen Bay, China. Regional Studies in Marine Science 40, 101462 (2020).Article 

Google Scholar 
Xu, X., Zhang, L. & Wei, C. Whistles of Indo-Pacific humpback dolphins (Sousa chinensis). AIP Conference Proceedings 1495, 556–562 (2012).Article 
ADS 

Google Scholar 
Wang, Z., Fang, L., Shi, W., Wang, K. & Wang, D. Whistle characteristics of free-ranging Indo-Pacific humpback dolphins (Sousa chinensis) in Sanniang Bay, China. The Journal of the Acoustical Society of America 133, 2479–2489 (2013).Article 
ADS 
PubMed 

Google Scholar 
Wang, Z.-T. et al. Apparent source levels and active communication space of whistles of free-ranging Indo-Pacific humpback dolphins (Sousa chinensis) in the Pearl River Estuary and Beibu Gulf, China. PeerJ 4, e1695 (2016).Article 
PubMed 
PubMed Central 

Google Scholar 
Sims, P. Q., Vaughn, R., Hung, S. K. & Würsig, B. Sounds of Indo-Pacific humpback dolphins (Sousa chinensis) in West Hong Kong: A preliminary description. The Journal of the Acoustical Society of America 131, EL48–EL53 (2012).Article 
ADS 
PubMed 

Google Scholar 
Dong, L. et al. Whistles emitted by Indo-Pacific humpback dolphins (Sousa chinensis) in Zhanjiang waters, China. The Journal of the Acoustical Society of America 145, 3289–3298 (2019).Article 
ADS 
PubMed 

Google Scholar 
Dong, L. et al. Whistle characteristics of a newly recorded Indo-Pacific humpback dolphin () population in waters southwest of Hainan Island, China, differ from other humpback dolphin populations. Marine Mammal Science 37, 1341–1362 (2021).Article 

Google Scholar 
Wang, X. et al. Three types of pulsed signal trains emitted by Indo-Pacific humpback dolphins (Sousa chinensis) in Beibu Gulf, South China Sea. Front. Mar. Sci. 9, (2022).Fang, L. et al. Echolocation signals of free-ranging Indo-Pacific humpback dolphins (Sousa chinensis) in Sanniang Bay, China. The Journal of the Acoustical Society of America 138, 1346–1352 (2015).Article 
ADS 
PubMed 

Google Scholar 
Fang, L. et al. The echolocation transmission beam of free-ranging Indo-Pacific humpback dolphins (Sousa chinensis). The Journal of the Acoustical Society of America 142, 771–779 (2017).Article 
ADS 
PubMed 
PubMed Central 

Google Scholar 
Au, W. W. L., Branstetter, B. K., Benoit-Bird, K. J. & Kastelein, R. A. Acoustic basis for fish prey discrimination by echolocating dolphins and porpoises. The Journal of the Acoustical Society of America 126, 460–467 (2009).Article 
ADS 
PubMed 

Google Scholar 
Xiang, W. et al. Underwater target classification based on the combination of dolphin click trains and convolutional neural networks. The Journal of the Acoustical Society of America 157, 647–658 (2025).Article 
ADS 
PubMed 

Google Scholar 
Jiang, J. et al. Covert underwater communication based on combined encoding of diverse time-frequency characteristics of sperm whale clicks. Applied Acoustics 171, 107660 (2021).Article 

Google Scholar 
Ahn, J., Park, G.-H., Kim, W., Kim, H.-M. & Lee, D.-H. Analysis of Research Trends and Technological Maturity of Biomimetic Underwater Communication. IEEE J. Oceanic Eng. 50, 1676–1702 (2025).Article 
ADS 

Google Scholar 
Bazúa-Durán, C. & Au, W. W. L. The whistles of Hawaiian spinner dolphins. The Journal of the Acoustical Society of America 112, 3064–3072 (2002).Article 
ADS 
PubMed 

Google Scholar 
Marley, S. A., Erbe, C. & Kent, C. P. S. Underwater recordings of the whistles of bottlenose dolphins in Fremantle Inner Harbour, Western Australia. Sci Data 4, 170126 (2017).Article 
PubMed 
PubMed Central 

Google Scholar 
Kamminga, C. Investigations on cetacean sonar X: A comparative analysis of underwater echolocation clicks of Inia spp. and Sotalia spp. Aquatic Mammals 19, 31 (1993).
Google Scholar 
Kamminga, C. & Stuart, A. B. C. Wave shape estimation of delphinid sonar signals, a parametric model approach. Acoust. lett 19, 70–76 (1995).
Google Scholar 
Madhusudhana, S., Gavrilov, A. & Erbe, C. Automatic detection of echolocation clicks based on a Gabor model of their waveform. The Journal of the Acoustical Society of America 137, 3077–3086 (2015).Article 
ADS 
PubMed 

Google Scholar 
Kaiser, J. F. On a simple algorithm to calculate the ‘energy’ of a signal. in International Conference on Acoustics, Speech, and Signal Processing 381–384 vol.1, https://doi.org/10.1109/ICASSP.1990.115702 (1990).Di Nardo, F., De Marco, R., Lucchetti, A. & Scaradozzi, D. A WAV file dataset of bottlenose dolphin whistles, clicks, and pulse sounds during trawling interactions. Sci Data 10, 650 (2023).Article 
PubMed 
PubMed Central 

Google Scholar 
Madsen, P., Kerr, I. & Payne, R. Source parameter estimates of echolocation clicks from wild pygmy killer whales (Feresa attenuata) (L). Journal of The Acoustical Society of America – J ACOUST SOC AMER 116, 1909–1912 (2004).Article 
ADS 
CAS 

Google Scholar 
Madsen, P. T. & Wahlberg, M. Recording and quantification of ultrasonic echolocation clicks from free-ranging toothed whales. Deep Sea Research Part I: Oceanographic Research Papers 54, 1421–1444 (2007).Article 
ADS 

Google Scholar 
Fu, W. J. et al. Acoustic recordings of underwater vocalizations of Indo-Pacific humpback dolphins in Xiamen Bay, China. figshare https://doi.org/10.6084/m9.figshare.29143727 (2025).Oswald, J., Erbe, C., Gannon, W., Madhusudhana, S. & Thomas, J. Detection and Classification Methods for Animal Sounds. in 269–317, https://doi.org/10.1007/978-3-030-97540-1_8 (2022).Kandia, V. & Stylianou, Y. Detection of sperm whale clicks based on the Teager–Kaiser energy operator. Applied Acoustics 67, 1144–1163 (2006).Article 

Google Scholar 
Klinck, H. & Mellinger, D. K. The energy ratio mapping algorithm: A tool to improve the energy-based detection of odontocete echolocation clicks. The Journal of the Acoustical Society of America 129, 1807–1812 (2011).Article 
ADS 
PubMed 

Google Scholar 
Yang, W., Luo, W. & Zhang, Y. Automatic detection method for monitoring odontocete echolocation clicks. Electronics Letters 53, 367–368 (2017).Article 
ADS 

Google Scholar 
Yang, W., Luo, W. & Zhang, Y. Classification of odontocete echolocation clicks using convolutional neural network. The Journal of the Acoustical Society of America 147, 49–55 (2020).Article 
ADS 
PubMed 

Google Scholar 
Liu, J., Yang, X., Wang, C. & Tao, Y. A Convolution Neural Network for Dolphin Species Identification Using Echolocation Clicks Signal. 4, https://doi.org/10.1109/ICSPCC.2018.8567796 (2018).Luo, W., Yang, W., Song, Z. & Zhang, Y. Automatic Species Recognition Using Echolocation Clicks from Odontocetes. 5, https://doi.org/10.1109/ICSPCC.2017.8242503 (2017).Download referencesAcknowledgementsThis work was funded by the National Natural Science Foundation of China (Grant Nos: 62231011), the Scientific Research Foundation of Third Institute of Oceanography, Ministry of Natural Resources (No. 2020017), the Natural Science Foundation of Fujian Province of China (Nos: 2024J01019), the marine ecological early warning and monitoring Foundation of Ministry of Natural Resources (No. S-HR04-230701-24), Open Fund Project of Hanjiang National Laboratory (No.KF2024030), the National Key R&D Program of China (grant number: 2024YFD2401402), the fund (No. Pilab2409) from State key laboratory of precision measuring technology and instruments (Tianjin University), and the MEL-RLAB Joint Fund for Marine Science & Technology Innovation (Grant number: M&R202402).Author informationAuthors and AffiliationsKey Laboratory of Underwater Acoustic Communication and Marine Information Technology of the Ministry of Education, College of Ocean and Earth Sciences, Xiamen University, Xiamen, 361005, ChinaWeijie Fu, Xuming Peng, Fei Zhang, Chuang Zhang, Wenjie Xiang, Zhongchang Song & Yu ZhangThird Institute of Oceanography, Ministry of Natural Resources, Xiamen, 361005, ChinaFuxing WuObservation and Research Station of Coastal Wetland Ecosystem in Beibu Gulf, Ministry of Natural Resources, Beihai, 536015, ChinaFuxing WuAuthorsWeijie FuView author publicationsSearch author on:PubMed Google ScholarXuming PengView author publicationsSearch author on:PubMed Google ScholarFuxing WuView author publicationsSearch author on:PubMed Google ScholarFei ZhangView author publicationsSearch author on:PubMed Google ScholarChuang ZhangView author publicationsSearch author on:PubMed Google ScholarWenjie XiangView author publicationsSearch author on:PubMed Google ScholarZhongchang SongView author publicationsSearch author on:PubMed Google ScholarYu ZhangView author publicationsSearch author on:PubMed Google ScholarContributionsWeijie Fu: writing – original draft, writing – review & editing, methodology, field investigation, data curation, formal analysis, validation. Xuming Peng: field investigation, writing – review & editing, data curation, software. Fuxing Wu: funding acquisition, field investigation, supervision, writing – review & editing. Fei Zhang: field investigation, data curation. Chuang Zhang: field investigation, data curation. Wenjie Xiang: data curation, software. Zhongchang Song: funding acquisition, field investigation, supervision, writing – review & editing. Yu Zhang: funding acquisition, supervision, writing – review & editing.Corresponding authorsCorrespondence to
Fuxing Wu or Zhongchang Song.Ethics declarations

Competing interests
The authors declare no competing interests.

Additional informationPublisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
Reprints and permissionsAbout this articleCite this articleFu, W., Peng, X., Wu, F. et al. Acoustic recordings of underwater vocalizations of Indo-Pacific humpback dolphins in Xiamen Bay, China.
Sci Data 12, 1960 (2025). https://doi.org/10.1038/s41597-025-06253-5Download citationReceived: 09 June 2025Accepted: 03 November 2025Published: 18 December 2025Version of record: 18 December 2025DOI: https://doi.org/10.1038/s41597-025-06253-5Share this articleAnyone you share the following link with will be able to read this content:Get shareable linkSorry, a shareable link is not currently available for this article.Copy shareable link to clipboard
Provided by the Springer Nature SharedIt content-sharing initiative More