Ecology

1.Vareda, J. P., Valente, A. J. M. & Durães, L. Assessment of heavy metal pollution from anthropogenic activities and remediation strategies: A review. J. Environ. Manage. 246, 101–118 (2019).CAS 
PubMed 

Google Scholar 
2.Chanamé, F., Custodio, M., Poma-Chávez, C. & Huamán, A. Nutrient concentrations and trophic state of three Andean lakes from Junín, Perú. Rev. Ambient Agua 15, 1–9 (2020).
Google Scholar 
3.Bhardwaj, R., Gupta, A. & Garg, J. K. Evaluation of heavy metal contamination using environmetrics and indexing approach for River Yamuna, Delhi stretch, India. Water Sci. 31, 52–66 (2017).
Google Scholar 
4.Custodio, M. et al. Human risk from exposure to heavy metals and arsenic in water from rivers with mining influence in the Central Andes of Peru. Water (Switzerland) 12, 1–20 (2020).
Google Scholar 
5.Arisekar, U., Jeya, R., Shalini, R. & Jeyasekaran, G. Human health risk assessment of heavy metals in aquatic sediments and freshwater fish caught from Thamirabarani River, the Western Ghats of South Tamil Nadu. Mar. Pollut. Bull. 159, 111496 (2020).CAS 
PubMed 

Google Scholar 
6.Chabukdhara, M. & Nema, A. K. Assessment of heavy metal contamination in Hindon River sediments: A chemometric and geochemical approach. Chemosphere 87, 945–953 (2012).CAS 
PubMed 
ADS 

Google Scholar 
7.Chai, L. et al. Heavy metals and metalloids in the surface sediments of the Xiangjiang River, Hunan, China: Distribution, contamination, and ecological risk assessment. Environ. Sci. Pollut. Res. 24, 874–885 (2017).CAS 

Google Scholar 
8.Liu, T. T. & Yang, H. Comparative analysis of the total and active bacterial communities in the surface sediment of Lake Taihu. FEMS Microbiol. Ecol. 96, 1–11 (2020).CAS 
ADS 

Google Scholar 
9.Custodio, M. et al. Evaluation of surface sediment quality in rivers with fish farming potential (Peru) using indicators of contamination, accumulation and ecological risk of heavy metals and arsenic. J. Ecol. Eng. 22, 78–87 (2021).
Google Scholar 
10.Zhang, Z. et al. Assessment of heavy metal contamination, distribution and source identification in the sediments from the Zijiang River, China. Sci. Total Environ. 645, 235–243 (2018).CAS 
PubMed 
ADS 

Google Scholar 
11.Sojka, M., Jaskula, J. & Siepak, M. Heavy metals in bottom sediments of reservoirs in the lowland area of western Poland: Concentrations, distribution, sources and ecological risk. Water (Switzerland) 11, 1–20 (2018).
Google Scholar 
12.Xu, Z., Te, S. H., Xu, C., He, Y. & Gin, K. Y. H. Variations of bacterial community composition and functions in an estuary reservoir during spring and summer alternation. Toxins (Basel) 10, 1–22 (2018).CAS 

Google Scholar 
13.Xiao, F. et al. The impact of anthropogenic disturbance on bacterioplankton communities during the construction of Donghu Tunnel (Wuhan, China). Microb. Ecol. 77, 277–287 (2019).CAS 
PubMed 

Google Scholar 
14.Wang, B. et al. Bacterial community responses to tourism development in the Xixi National Wetland Park, China. Sci. Total Environ. 720, 137570 (2020).CAS 
PubMed 
ADS 

Google Scholar 
15.Deng, W. et al. Heavy metals, antibiotics and nutrients affect the bacterial community and resistance genes in chicken manure composting and fertilized soil. J. Environ. Manage. 257, 109980 (2020).CAS 
PubMed 

Google Scholar 
16.Gubelit, Y. et al. Nutrient and metal pollution of the eastern Gulf of Finland coastline: Sediments, macroalgae, microbiota. Sci. Total Environ. 550, 806–819 (2016).CAS 
PubMed 
ADS 

Google Scholar 
17.Wang, J. et al. Contribution of heavy metal in driving microbial distribution in a eutrophic river. Sci. Total Environ. 712, 136295 (2020).CAS 
PubMed 
ADS 

Google Scholar 
18.Liao, H. et al. Profiling microbial communities in a watershed undergoing intensive anthropogenic activities. Sci. Total Environ. 647, 1137–1147 (2019).CAS 
PubMed 
ADS 

Google Scholar 
19.Liu, J. et al. Spatiotemporal dynamics of the archaeal community in coastal sediments: Assembly process and co-occurrence relationship. ISME J. 14, 1463–1478 (2020).PubMed 
PubMed Central 

Google Scholar 
20.Liao, H., Yen, J. Y., Guan, Y., Ke, D. & Liu, C. Differential responses of stream water and bed sediment microbial communities to watershed degradation. Environ. Int. 134, 105198 (2020).CAS 
PubMed 

Google Scholar 
21.Song, H., Li, Z., Du, B., Wang, G. & Ding, Y. Bacterial communities in sediments of the shallow Lake Dongping in China. J. Appl. Microbiol. 112, 79–89 (2012).CAS 
PubMed 

Google Scholar 
22.Ligi, T. et al. Characterization of bacterial communities in soil and sediment of a created riverine wetland complex using high-throughput 16S rRNA amplicon sequencing. Ecol. Eng. 72, 56–66 (2014).
Google Scholar 
23.Wilmes, P. et al. Natural acidophilic biofilm communities reflect distinct organismal and functional organization. ISME J. 3, 266–270 (2009).CAS 
PubMed 

Google Scholar 
24.Mavromatis, K. et al. Use of simulated data sets to evaluate the fidelity of metagenomic processing methods. Nat. Methods. 4, 495–500 (2007).CAS 
PubMed 

Google Scholar 
25.Yuan, X., Zhang, L., Li, J., Wang, C. & Ji, J. Sediment properties and heavy metal pollution assessment in the river, estuary and lake environments of a fluvial plain, China. CATENA 119, 52–60 (2014).CAS 

Google Scholar 
26.Lin, Q., Liu, E., Zhang, E., Li, K. & Shen, J. Spatial distribution, contamination and ecological risk assessment of heavy metals in surface sediments of Erhai Lake, a large eutrophic plateau lake in southwest China. CATENA 145, 193–203 (2016).CAS 

Google Scholar 
27.Guo, T. et al. Distribution of arsenic and its biotransformation genes in sediments from the East China Sea. Environ. Pollut. 253, 949–958 (2019).CAS 
PubMed 

Google Scholar 
28.Taylor, S. R. & Mclennan, S. M. The geochemical the continental evolution crust. Rev. Miner. Geochem. 33, 241–265 (1995).
Google Scholar 
29.Lastauskienė, E. et al. The impact of intensive fish farming on pond sediment microbiome and antibiotic resistance gene composition. Front. Vet. Sci. 8, 1–12 (2021).
Google Scholar 
30.Ragab, S., Sikaily, A. E., Nemr, A. E. & Sea, R. Concentrations and sources of pesticides and PCBs in surficial sediments of the Red Sea coast, Egypt. Egypt. J. Aquat. Res. 42, 365–374 (2016).
Google Scholar 
31.Kavita, V. & Pandey, J. Heavy metal accumulation in surface sediments of the Ganga River (India): Speciation, fractionation, toxicity, and risk assessment. Environ. Monit. Assess. 191, 20 (2019).
Google Scholar 
32.Haghnazar, H. et al. Chemosphere Potentially toxic elements contamination in surface sediment and indigenous aquatic macrophytes of the Bahmanshir River, Iran: Appraisal of phytoremediation capability. 285, (2021).33.Perera, P. C. T., Sundarabarathy, T. V., Sivananthawerl, T., Kodithuwakku, S. P. & Edirisinghe, U. Arsenic and cadmium contamination in water, sediments and fish is a consequence of paddy cultivation: Evidence of river pollution in Sri Lanka. Achiev. Life Sci. 10, 144–160 (2016).
Google Scholar 
34.Kalantzi, I., Rico, A., Mylona, K., Pergantis, S. A. & Tsapakis, M. Fish farming, metals and antibiotics in the eastern Mediterranean Sea: Is there a threat to sediment wildlife?. Sci. Total Environ. 764, 142843 (2021).CAS 
PubMed 
ADS 

Google Scholar 
35.Monroy, M., Maceda-Veiga, A. & de Sostoa, A. Metal concentration in water, sediment and four fish species from Lake Titicaca reveals a large-scale environmental concern. Sci. Total Environ. 487, 233–244 (2014).CAS 
PubMed 
ADS 

Google Scholar 
36.Rodbell, D. T., Delman, E., Abbott, M., Besonen, M. & Tapia, P. The heavy metal contamination of Lake Junín National Reserve, Peru: An unintended consequence of the juxtaposition of hydroelectricity and mining. GSA Today 24, 4–10 (2014).
Google Scholar 
37.Ni, C. et al. High concentrations of bioavailable heavy metals impact freshwater sediment microbial communities. Ann. Microbiol. 66, 1003–1012 (2016).CAS 

Google Scholar 
38.Huang, W. et al. Comparison among the microbial communities in the lake, lake wetland, and estuary sediments of a plain river network. Microbiologyopen https://doi.org/10.1002/mbo3.644 (2018).Article 
PubMed 
PubMed Central 

Google Scholar 
39.Abia, A. L. K., Alisoltani, A., Keshri, J. & Ubomba-Jaswa, E. Metagenomic analysis of the bacterial communities and their functional profiles in water and sediments of the Apies River, South Africa, as a function of land use. Sci. Total Environ. 616–617, 326–334 (2018).PubMed 
ADS 

Google Scholar 
40.Guo, X. et al. Characteristics of microbial community indicate anthropogenic impact on the sediments along the Yangtze Estuary and its coastal area, China. Sci. Total Environ. 648, 306–314 (2019).CAS 
PubMed 
ADS 

Google Scholar 
41.Betiku, O. C. et al. Evaluation of microbial diversity of three recreational water bodies using 16S rRNA metagenomic approach. Sci. Total Environ. 771, 144773 (2021).CAS 
PubMed 
ADS 

Google Scholar 
42.Zhang, T. et al. Suspended particles phoD alkaline phosphatase gene diversity in large shallow eutrophic Lake Taihu. Sci. Total Environ. 728, 138615 (2020).CAS 
PubMed 
ADS 

Google Scholar 
43.Shen, M. et al. Trophic status is associated with community structure and metabolic potential of planktonic microbiota in Plateau Lakes. Front. Microbiol. 10, 1–15 (2019).
Google Scholar 
44.Quero, G. M., Cassin, D., Botter, M., Perini, L. & Luna, G. M. Patterns of benthic bacterial diversity in coastal areas contaminated by heavy metals, polycyclic aromatic hydrocarbons (PAHs) and polychlorinated biphenyls (PCBs). Front. Microbiol. 6, 1–15 (2015).
Google Scholar 
45.Wang, Y. et al. Comparison of the levels of bacterial diversity in freshwater, intertidal wetland, and marine sediments by using millions of illumina tags. Appl. Environ. Microbiol. 78, 8264–8271 (2012).CAS 
PubMed 
PubMed Central 
ADS 

Google Scholar 
46.Long, Y. et al. The response of microbial community structure and sediment properties to anthropogenic activities in Caohai wetland sediments. Ecotoxicol. Environ. Saf. 211, 111936 (2021).CAS 
PubMed 

Google Scholar 
47.Yao, X., Zhang, J., Tian, L. & Guo, J. The effect of heavy metal contamination on the bacterial community structure at Jiaozhou Bay, China. Braz. J. Microbiol. 48, 71–78 (2017).CAS 
PubMed 

Google Scholar 
48.Hur, M. & Park, S. J. Identification of microbial profiles in heavy-metal-contaminated soil from full-length 16s rRNA reads sequenced by a pacbio system. Microorganisms 7, 25 (2019).
Google Scholar 
49.Zhuang, M., Sanganyado, E., Li, P. & Liu, W. Distribution of microbial communities in metal-contaminated nearshore sediment from Eastern Guangdong, China. Environ. Pollut. 250, 482–492 (2019).CAS 
PubMed 

Google Scholar 
50.Gu, Y. et al. Degradation shaped bacterial and archaeal communities with predictable taxa and their association patterns in Zoige wetland at Tibet plateau. Sci. Rep. 8, 1–11 (2018).ADS 

Google Scholar 
51.Newton, R. J., Jones, S. E., Eiler, A., McMahon, K. D. & Bertilsson, S. A guide to the natural history of freshwater lake bacteria. Microbiol. Mol. Biol. Rev. 75, 25 (2011).
Google Scholar 
52.Hu, A. et al. Strong impact of anthropogenic contamination on the co-occurrence patterns of a riverine microbial community. Environ. Microbiol. 19, 4993–5009 (2017).CAS 
PubMed 

Google Scholar 
53.Ren, Z. et al. Taxonomic and functional differences between microbial communities in Qinghai Lake and its input streams. Front. Microbiol. 8, 1–14 (2017).
Google Scholar 
54.Yin, X. et al. Cadmium isotope constraints on heavy metal sources in a riverine system impacted by multiple anthropogenic activities. Sci. Total Environ. 750, 141233 (2021).CAS 
PubMed 
ADS 

Google Scholar 
55.Yan, C. et al. Integrating high-throughput sequencing and metagenome analysis to reveal the characteristic and resistance mechanism of microbial community in metal contaminated sediments. Sci. Total Environ. 707, 136116 (2020).CAS 
PubMed 
ADS 

Google Scholar 
56.Coclet, C. et al. Trace metal contamination impacts predicted functions more than structure of marine prokaryotic biofilm communities in an anthropized coastal area. Front. Microbiol. 12, 1–16 (2021).
Google Scholar 
57.Esri Inc. ArcMap 10.8. Esri Inc. (2020). https://www.esri.com/en-us/arcgis/products/arcgis-pro/overview.58.Avalos, G. et al. Climate Change in the Mantaro River Basin (MINEN, 2013).
Google Scholar 
59.APHA. Standard methods for the examination of water and wastewater. Stand. Methods 541, 25 (2012).
Google Scholar 
60.Singh, H., Pandey, R., Singh, S. K. & Shukla, D. N. Assessment of heavy metal contamination in the sediment of the River Ghaghara, a major tributary of the River Ganga in Northern India. Appl. Water Sci. 7, 4133–4149 (2017).CAS 
ADS 

Google Scholar 
61.El-Amier, Y. A., Elnaggar, A. A. & El-Alfy, M. Evaluation and mapping spatial distribution of bottom sediment heavy metal contamination in Burullus Lake, Egypt. Egypt. J. Basic Appl. Sci. https://doi.org/10.1016/j.ejbas.2016.09.005 (2016).Article 

Google Scholar 
62.Miller, D. N., Bryant, J. E., Madsen, E. L. & Ghiorse, W. C. Evaluation and optimization of DNA extraction and purification procedures for soil and sediment samples. Appl. Environ. Microbiol. 65, 4715–4724 (1999).CAS 
PubMed 
PubMed Central 
ADS 

Google Scholar 
63.Custodio, M. et al. Metagenomic data on the composition of bacterial communities in lake environment sediments for fish farming by next generation Illumina sequencing. Data Br. 32, 106228 (2020).
Google Scholar 
64.Bolger, A. M., Lohse, M. & Usadel, B. Trimmomatic: A flexible trimmer for Illumina sequence data. Bioinformatics 30, 2114–2120 (2014).CAS 
PubMed 
PubMed Central 

Google Scholar 
65.Wood, D. E. & Salzberg, S. L. Kraken: Ultrafast metagenomic sequence classification using exact alignments. Genome Biol. 15, 25 (2014).
Google Scholar 
66.Edgar, R. C. UPARSE: Highly accurate OTU sequences from microbial amplicon reads. Nat. Methods 10, 996–998 (2013).CAS 
PubMed 

Google Scholar 
67.Gan, Y. et al. Multiple factors impact the contents of heavy metals in vegetables in high natural background area of China. Chemosphere 184, 1388–1395 (2017).CAS 
PubMed 
ADS 

Google Scholar 
68.Diallo, M. D. et al. Polymerase chain reaction denaturing gradient gel electrophoresis analysis of the N2-fixing bacterial diversity in soil under Acacia tortilis ssp. raddiana and Balanites aegyptiaca in the dryland part of Senegal. Environ. Microbiol. 6, 400–415 (2004).CAS 

Google Scholar 
69.R Core Team. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna (2020). https://www.R-project.org/.70.Li, C. et al. Effects of heavy metals on microbial communities in sediments and establishment of bioindicators based on microbial taxa and function for environmental monitoring and management. Sci. Total Environ. 749, 141555 (2020).CAS 
PubMed 
ADS 

Google Scholar 
71.Murtaza, N. et al. Analysis of the effects of dietary pattern on the oral microbiome of elite endurance athletes. Nutrients 11, 1–12 (2019).MathSciNet 

Google Scholar  More

Species dataSpecies CPUE data were obtained from the National Oceanographic and Atmospheric Administration (NOAA) Northeast Fishery Science Center (NEFSC) U.S. NES bottom trawl survey, which, for almost 50 years, has collected abundance data for over 250 fish species in the spring and fall. The survey employs a stratified random design, with stations allocated proportionally to the stratum area. A 12 mm mesh coded liner is used to retain small-bodied and juvenile fish. All fish caught are weighed and counted18. We downloaded the data from OceanAdapt.com, which calibrates the CPUE for each species depending on survey ship. We cleaned the data for the years from 1998 to 2020, excluding years prior to 1997 due to many missing values for chlorophyll (Chla). We only included strata that were consistently sampled in the spring and fall. To account for the seasonal migrations of many of the studied species, we modeled spring and fall seasons separately. We present the results for the fall CPUE, with the spring results and presence/absence results in the supplemental materials. We selected species that were present in at least 400 tows and with a biomass of at least 0.5 kg/tow (CPUE) in more than 100 tows. Finally, we removed roughly 400 tows per season with missing environmental covariates (see below). In the fall, we selected 30 species with 5217 observations, and in the spring, we selected 24 species with 5935 observations (see Supplemental Tables S1, S2).Environmental dataThe study region includes Southern New England and The Gulf of Maine. We selected environmental covariates known to influence marine fish distributions and abundances. Depth, temperature (bottom and surface) and salinity (bottom and surface) were measured in situ during trawl surveys. Missing values were augmented with the data-assimilative HYbrid Coordinate Ocean Model (HYCOM) daily and then monthly data. HYCOM is an oceanographic model that produces 32 vertical layers including ocean temperature, salinity, sea surface height, and wind stress as well as other 3- and 4-dimensional variables. The system uses the Navy Coupled Ocean Data Assimilation (NCODA) system19 for data assimilation. NCODA uses the model forecast as a first guess in a multivariate optimal interpolation (MVOI) scheme and assimilates available satellite altimeter observations (along track obtained via the Naval Oceanographic Office Altimeter Data Fusion Center satellite) and in situ sea surface temperature as well as available in situ vertical temperature and salinity profiles from expendable bathythermographs, Argo floats, and moored buoys20. Seven HYCOM models (HYCOM + NCODA Global 1/12° Reanalysis GLBu0.08 Expts 19.0, 19.1, 90.9, 91.0, 91.1, 91.2) were temporally concatenated to create a continuous dataset of BT and salinity, ranging from 1992 to 2017. These model runs differed slightly in their configurations (time steps, advection scheme, mixing, vertical structure, slight change in NCODA, and MVOI transition to 3-dimensional analysis in 2013), but the differences are not expected to influence the applicability of the output21. The numbers of filled in missing values were 787 (7.0%) surface salinity (SSAL), 735 (6.5%) surface temperature (SST), 809 (7.2%) bottom temperature (BT), and 850 (7.6%) bottom salinity (BSAL). Chla was obtained from the MODIS satellite (monthly rasters from 2003 to 2019) on a monthly time step22, with missing values filled using the SeaWIFS satellite23 (1998 to 2009). Temperature, salinity and Chla data that were not collected in situ were downloaded using Google Earth Engine (HYCOM and MODIS)24. Benthic substrate (grain size in mm, referred to as SEDSIZE), subregion (Gulf of Maine or Southern New England), benthic land position (high, low, mid), and seabed form data (depression, high flat, high slope, low slope, mid flat, side slope, steep) were obtained from the Nature Conservancy’s Northwest Atlantic Marine Ecoregional Assessment25 (Supplemental Fig. S1).GJAMTo study the influence of the environmental covariates on the joint distribution of marine fish and invertebrate species we use the generalized joint attribute model (GJAM)12 and the corresponding R package (Version 2.5)26. Briefly, this multivariate Bayesian model allows us to jointly model the marine fish community and accounts for responses to the environment that can include combinations of continuous and discrete responses (e.g., CPUE and zeros) and the dependence between species. GJAM returns all parameters on the observation scale, in this case, CPUE and presence-absence. Products of model fitting include a species‐by‐species covariance matrix (Σ), species responses to predictor variables (B), and predicted responses. The species‐by‐species covariance matrix Σ captures residual codependence between species after removing the main structure explained by the model (also referred to as the residual correlation matrix). As a result, Σ allows for conditional prediction of some species under different scenarios for the abundances of others27.CPUE is termed continuous abundance (CA) data in GJAM, meaning that observations are continuous with discrete zeros. Let yis be the CPUE for species s at location i. For CA data GJAM expands the tobit model for (univariate) regression to the multivariate setting, where a latent variable wis is equal to yis when yis is positive and zero otherwise,$$y_{i,s}^{0} = left{ {begin{array}{*{20}l} {w_{is} ,} hfill & {w_{is} > 0quad {text{continuous}}} hfill \ {0,} hfill & {w_{is} le 0quad {text{discrete zero}}} hfill \ end{array} } right.$$
(1)
The length-S vector of all species responses wi is continuous on the real line, and thus can be modeled with a multivariate normal. The model for wi is$$begin{aligned} left. {{mathbf{w}}_{i} } right|{mathbf{x}}_{i, } {mathbf{y}}_{i} & sim ,MVNleft( {{varvec{mu}}_{i} ,{Sigma }} right) times mathop prod limits_{s = 1}^{S} {mathcal{I}}_{is} \ u_{{varvec{i}}} & = {mathbf{B}}^{prime } {mathbf{x}}_{{varvec{i}}} \ {mathcal{I}}_{is} & = mathop prod limits_{k in C} I_{is,k}^{{Ileft( {y_{is} = k} right)}} left( {1 – I_{is,k} } right)^{{Ileft( {y_{is} ne k} right)}} \ end{aligned}$$
(2)
$$begin{aligned} {mathcal{I}}_{is} & = I(w_{is} le 0)^{{Ileft( {y_{is} = 0} right)}} left[ {1 – Ileft( {w_{is} le 0} right)} right]^{{Ileft( {y_{is} > 0} right)}} \ & quad I(w_{is} > 0)^{{Ileft( {y_{is} > 0} right)}} left[ {(1 – I(w_{is} > 0)} right]^{{Ileft( {y_{is} = 0} right)}} \ end{aligned}$$where the indicator function (I(cdot )) is equal to 1 when its argument is true and 0 otherwise. For presence-absence data, ({mathbf{p}}_{{varvec{i}}{varvec{s}}}boldsymbol{ }=boldsymbol{ }left(-boldsymbol{infty },boldsymbol{ }0,boldsymbol{ }boldsymbol{infty }right).) This is equivalent to Chib and Greenberg’s28 probit model which can be written as ({mathcal{I}}_{is}=I({w}_{is} >{0)}^{Ileft({y}_{is} >0right)}I({w}_{is}le {0)}^{1-{y}_{is}}).The mean vector ({varvec{mu}}_{i} = {mathbf{B}}^{prime } {mathbf{x}}_{{varvec{i}}}) contains the Q × S matrix of coefficients B and the length-Q design vector xi. Σ is a S × S covariance matrix. There is a correlation matrix associated with Σ,$${mathbf{R}}_{{S,S^{prime } }} = frac{{{{varvec{Sigma}}}_{{S,S^{prime } }} }}{{sqrt {{{varvec{Sigma}}}_{S,S} {{varvec{Sigma}}}_{{S^{prime } ,S^{prime } }} } }}$$
(3)
The predictive distribution is obtained as$$left[tilde{Y }left| tilde{X }right.right]=int left[ tilde{Y }left| tilde{X }right.,widehat{theta }right]left[widehat{theta } left|X, Yright.right]$$
(4)
The integrand contains the likelihood (Eq. (2)) followed by the posterior distribution for parameters, (widehat{theta }= left{widehat{mathbf{B}},boldsymbol{ }widehat{{varvec{Sigma}}}right}). Input (tilde{X }) can equal X (in-sample prediction) or not (out-of-sample prediction). We fitted both CPUE (continuous abundance) and presence-absence versions of the model. As a Bayesian method, GJAM provides probabilistic estimates of parameters with full dependence in data, including jointly distributed species. Model fitting is performed using Gibbs sampling, which is a Markov chain Monte Carlo (MCMC) technique.The sensitivity of an individual response variable s to an individual predictor q is given by the coefficient βqs (individual coefficients from the B matrix). The sensitivity that applies to the full response matrix is given by$${mathbf{f}} = diagleft( {{mathbf{B}}{Sigma }^{ – 1} {mathbf{B}}^{prime } } right)$$
(5)
The Q × S matrix B contains relationships of each species to the environment, the “signal”, but not to each another. Matrix E summarizes species similarities in terms of their response to an environment (stackrel{sim }{mathbf{x}}) and is given by$${mathbf{E}=mathbf{B}}^{boldsymbol{^{prime}}}mathbf{V}mathbf{B}$$
(6)
where V is a covariance matrix for (stackrel{sim }{mathbf{x}})(a vector of predictors) and contributes the environmental component of variation in (stackrel{sim }{mathbf{y}}). Similar species in E have similar columns in B. Those similarities and differences are amplified for predictors (stackrel{sim }{mathbf{x}}) with large variance. Conversely, species differences in B do not matter for variables in X that do not vary. The covariance in predictors could come from observed data, i.e., the variance of X (see12 for more details).Prior distributions for this study are non-informative. This is particularly helpful for the covariance, lending stability to Gibbs sampling and avoiding dominance by a prior. In cases this particular case, the direction of the prior effect may be known, but the magnitude is not.Variable selectionUnlike the familiar univariate setting, variable selection has to consider which species are included in the model. In a univariate model, there is one response and perhaps a number of potential predictor variables from which to choose. As in a univariate model, variable selection focuses on predictors held in the n by p design matrix X. Rather than a response vector, the multivariate model includes the n by S response matrix Y. Unlike the univariate model, the overall fit and predictive capacity depends not only on what is in X, but also on the species that are included in Y, each of which would be best explained by a different combination of variables. Rare species having no signal will not provide cross-correlations and thus can offer little learning from an analysis. For this reason, there may be no reason to include them in model fitting. Given that many species may be rare, and rare types will not be explained by the model, there will be decisions about what variables to include on both sides of the likelihood (i.e., predictors and responses).These considerations mean that simple rules for variable selection, such as the combination yielding the lowest DIC, may not be sensible. The combination of variables that yields the lowest DIC could miss variables that are important for subsets of species. In principle, one poorly-fitted species could dominate variable selection. The best model for responses ranging from rare to abundant will depend on precisely which species are included, both rare and abundant. Thus, in order to select variables, we utilize inverse prediction—predicting the environment from species – and the overall community sensitivity12.Inverse prediction provides a comprehensive estimate of the environmental importance for the entire community, because it determines the capacity of the community to predict (through the fitted model) the environment; it inverts the model12. A variable predicted by the community explains important variation in one to many species. A variable that is not predicted by the community does not explain important variation in any of them. To look at the importance of environmental variables for the entire community, we started with the saturated model that included the predictors BT, SST, depth, BSAL, SSAL, Chla, SEDSIZE, subregion, benthic position and an interaction between depth and BT, BSAL, SST and SSAL (Fig. 1a). Sensitivity was highest for the interaction between BT and depth and lowest for Chla and sediment size (see right subpanel on Fig. 1a for sensitivity). Inverse prediction confirmed that sediment size and Chla contribute little to community biomass, because the community cannot “predict” them (see left and middle subpanels on Fig. 1a for sensitivity). Inverse prediction results from a second model (Fig. 1b) showed that SSAL and the third model for benthic position also (Fig. 1c) contribute little to the community response. Using the combination of sensitivity and inverse prediction we obtained the final model that includes BT, depth, BSAL, SST, subregion and an interaction between depth and BT, BSAL and SST (Fig. 1d). Inverse prediction indicates that the CPUE predicts the environment well. In the final model, sensitivity is highest for depth. Subregion remains as a two-level factor and there is strong inverse prediction for that variable as well (Fig. 1d). In the variable-selection stage, each model was run on the entire fall dataset for 5000 iterations and a burn-in of 800. Inverse prediction results from the spring model indicated similar patterns; thus, the same variables were used for the spring and fall.Figure 1Inverse prediction and sensitivity for combinations of environmental parameters in GJAM. Starting with the most complicated model (a), sensitivity was highest for the interaction between BT and depth and lowest for Chla and sediment size (a). Inverse prediction confirms that sediment size and Chla contribute little to community biomass (a) and those are removed in the second model (b). SSAL contributes little to community response and are removed in the third model (c), The final model (d) includes terms that have strong inverse prediction and overall sensitivity. Inverse prediction for continuous and factor variables is on the left and center of each box, and overall sensitivity is on the right.Full size imageWe compare the model selected above using inverse prediction to a model selected using the more traditional method of out-of-sample prediction. For out-of-sample prediction, we fitted all combinations of 11 environmental variables (BT, BSAL, SST, SSAL, Chla, depth, sediment size, subregion, position, seabed form) plus interaction terms between depth and SEDSIZE, BT, BSAL, SST, SSAL and chlorophyll. These models were run with 1000 iterations and a burn-in of 400. All models included BT, BSAL, SST, SSAL, chlorophyll A and depth, as these variables have been shown to be important for these species. In total, 1,024 possible models were evaluated by training each potential model on 70% of the data (n = 3652 in the fall, n = 4155 in the spring), evaluating in-sample performance with DIC, and then testing out-of-sample performance on the remaining 30% (n = 1565 in the fall, n = 1780 in the spring). The 10 models with the lowest DIC in-sample were selected, and the final model was selected out of those 10 with the lowest out-of-sample R2. The selected model for fall CPUE had the following terms: ~ BT + depth + BSAL + SST + SSAL + chla + depth*BT + depth*SEDSIZE + depth*SSAL + depth*chla + SEDSIZE + Benthic position. Recall that inverse prediction selected a simpler model including the following terms: BT + depth + BSAL + SST + Subregion + depth*BT + depth*BSAL + depth*SST. The inclusion of SEDSIZE and benthic position in the model selected via out-of-sample prediction is probably a result of these predictor variables being important for a subset of species (i.e. benthic species29), but not the community as a whole. When we have a large number of response variables, as in this study, we need to consider the variables that are more important on a community level, rather than just for a few species. Thus, we use the model selected via inverse prediction for the remainder of the study.We fitted the selected model with 70% of the data for 20,000 iterations with a burn-in of 8,000 iterations (n = 3652 in the fall, n = 4155 in the spring). Out-of-sample prediction was performed on the remaining 30% (n = 1565 in the fall, n = 1780 in the spring) of the dataset and predicted versus observed values were evaluated (Supplemental Figs. S2 and S3) as well as residual versus fitted values (Supplemental Figs. S4 and S5). As has been shown in other research30,31, aggregating noisy predictions based on similar environmental preferences can improve performance, especially for larger datasets. Thus, we generated an aggregated data set that uses a k-means clustering of predictors (Supplemental Figs. S8 and S9). We performed the same analysis for the spring and the fall as well as with the presence absence data and recorded AUC as well as precision for each species (Supplemental Figs. S6 and S7). Precision is defined as the arithmetic mean of precision (proportion of predicted presences actually observed as presences) across all threshold values (at an interval of 0.01).Final modelWe ran the final model on 100% of the data with 20,000 iterations and a burn-in of 8000 iterations for the spring and fall for CPUE as well as presence absence for a total of 4 models. From the final model we obtained coefficients for the species-environment responses, β, covariance between species in how they respond to the environment E, and the residual correlation from the fitted model, R. We subtracted the absolute values from the presence/absence residual correlation matrix from the absolute values of the CPUE residual correlation matrix to observe where these results diverged. For MCMC chains and convergence of the final model as well as example models from both methods of variable selection see Supplemental Figs. S10–S12).Comparison to SSDMsWe built single species distribution models for each species in the form of GAMs using the mgcv package in R32. GAMs are a semiparametric extension of the generalized linear model and are a common modeling technique for species distribution modeling in this ecosystem33. For each species, we ran one GAM with CPUE as the response variable with a log-linked tweedie distribution that had penalized regression splines, a REML smoothing parameter with an outer Newton optimizer, 10 knots, and omitted NAs. We also ran GAMs for each species with a binary response variable indicating species presence with a binomial error distribution and a logit link function, penalized regression splines, a REML smoothing parameter with an outer Newton optimizer, 10 knots, and omitted NAs. We compared the out of sample observed versus predicted values for GAMs versus GJAM using RMSPE, R2, AUC, and precision. Root Mean Squared Prediction Error (RMSPE) is a measure of the average squared difference between the observed and predicted values, measured in the same units as the input data (kg/tow). R2 is a measure of the average squared difference between the observed and predicted values and is unitless. R2 is calculated as (1 − sum((predicted − observed)2)/sum((observed − mean(observed))2)) The ROC curve is a measure of model performance which plots true positive rate versus false positive rate, and the area under the ROC curve (AUC) provides a single measure of accuracy. A pairwise Wilcoxon test was used to compare means. We also compare the significance of predictors in both the GJAM model and GAM models. In this example, significance is defined for GJAM as a credible interval of the beta estimation that does not cross zero, and for the GAM as a p-value less than 0.0534.Spatial and temporal autocorrelationExamining the spatial and temporal autocorrelation of the modeled residuals can help specify missing endogenous (habitat selection or density dependence) and exogenous (covariate) effects that may be missing from the model. Thus, for each species modeled, we plot the spatial autocorrelation of residuals using a semi-variogram for the year 2015 and the temporal autocorrelation of the residuals using a partial autocorrelation function (PACF). We present the results for each species in the fall in the Supplemental materials (Supplemental Figs. S27–S57).All analysis and figure creation was performed in R version 3.6.235. Figures were created using the following R packages: ggplot236, ggpubr37, corrplot38, gridExtra39, cowplot40, lessR41, and ggcorrplot42. More

Study area and study populationData come from a long-term study of a pied flycatcher population breeding in nestboxes in central Spain (ca. 41°N, 3°W, 1200–1300 m.a.s.l.). The longitudinal data cover the period 1990–2016 (no data for 2003) and include records for 1436 males (yearly mean and SD: 107.4 and 34.2) and 1641 females (yearly mean and SD: 119.7 and 28.6). The study area consists of two plots in two different montane habitats separated by 1.1 km, including 237 nestboxes with an average occupancy rate around 54% (SD = 0.11). One habitat is an old deciduous oak (Quercus pyrenaica) forest, and the other one is a managed mixed coniferous (mainly Pinus sylvestris) forest. The nestboxes have remained in the same position since 1988 (pinewood) and 1995 (oakwood) (for details, see42,43).Fieldwork and data collectionNestboxes were regularly (every 3rd–4th day) checked during the breeding season (from mid-April to the beginning of July) to determine the date of the first egg laid, clutch size, hatching date, and the number of fledglings. Parents were captured with a nestbox trap while incubating (females) or feeding 8-day-old nestlings (both sexes; for details, see43 and marked with a numbered metal ring (both sexes). We used a unique combination of colour rings (males only) for individual identification before capture. Many breeding birds (53%) hatched in the nestboxes, and, therefore, their exact age was known44. Unringed breeders were aged as first-year or older based on plumage traits following ageing criteria described in44,]45. All nestlings were ringed at 13 days of age.Polygamous males were detected when captured and/or individually identified while repeatedly feeding young in two nests (see24 for details on capture protocol and mating status classification). We distinguished three classes of females according to their male mating status: (i) monogamous female, i.e. mated with a monogamous male; (ii) primary female, the first mated female of a polygynous male; and (iii) secondary female, the second mated female of a polygynous male. However, in some nests, it was not possible to know with certainty the mating status of the female (14.3% of times) or the male (3.7% of times, see below for how we dealt with this source of uncertainty).Ethics declarationThe study was reviewed by the ethical committees at the Doñana Biological Station and the Consejo Superior de Investigaciones Científicas headquarters (Spain) and adhered to Spain standards. All methods were carried out in accordance with relevant guidelines and regulations. Birds were caught and ringed with permission from the Spanish Ministry of Agriculture, Food, Fisheries, and Environment’s Ringing Office. The study complied with (Animal Research: Reporting of In Vivo Experiments) guidelines46.Multi-event capture-recapture modelsWe used multi-event capture-recapture (MECR hereafter) models47 to test, separately for females and males, how the mating status affected the probability of surviving (and not leaving the area permanently) and the probability of changing, or not, from one mating status to another. The MECR models accommodate uncertainty in state assignment by distinguishing between what is observed (the event) and what is inferred (the state). This approach allows estimating the effects of mating status on the parameters (e.g. probabilities of local survival and change in mating status) while accounting for the uncertainty, as outlined above, due to the unknown mating status of some captured individuals.MECR models are defined by three types of parameters: Initial State probabilities, Transition probabilities and Event probabilities (details in Appendices S5). As these parameter types may be broken into steps, we considered two Transition steps, Local survival and Mating Status Change, and two Event steps, Recapture and Mating Status Assignment. Accordingly, we considered the following parameters of the MECR model: (i) Initial State, the probability of being in a specific mating status at the first encounter (in our case the first known breeding event of an individual); (ii) Local survival, the probability of surviving and not emigrating permanently from the study area between year t and year t + 1; (iii) Mating Status Change, the probability that a live bird changes state between year t and t + 1; (iv) Recapture: the probability of recapture of a live and not permanently emigrated individual; (v) Mating Status Assignment: the probability that the mating status of a captured individual is ascertained in the field (assuming no state misclassification). In this study, we will use the term “parameter” to denote any of the probabilities (see i-v above) estimated in the MECR model. Also, note that, as is often the case, we cannot distinguish the probability of site fidelity from that of surviving. For simplicity, we will often use the term “survival” to refer to “local survival”.We used the encounter histories of all identified birds breeding in the study area at least once between 1990 and 2016. We ran separate analyses for each sex, considering four biological states for females: live monogamous breeder (MBF), live primary breeder (PBF), live secondary breeder (SBF) and dead or permanently emigrated (†); and five events, numbered as they appear in the encounter histories: (0) non-captured, (1) captured as a monogamous breeder, (2) captured as a primary breeder, (3) captured as a secondary breeder and (4) captured in an unknown mating status. Females of unknown mating status were those for which we did not know the mate’s identity after repeated identification attempts at the nestbox (see details in24). These females could be of any mating status, and the mate being absent (e.g. dead after pairing) or very sporadically visiting the nest. For males, however, we considered three biological states: live monogamous breeder (MBM), live polygynous breeder (PBM) and dead or permanently emigrated (†), mediated by four events: (0) non-captured, (1) captured as a monogamous breeder, (2) captured as a polygynous breeder, (3) captured in an unknown mating status. Males of unknown mating status were identified by reading their colour-rings combinations near a nestbox and not captured or seen again during the breeding season. For both sexes, we established two age classes: 1-year-old individuals (1-yo hereafter: 41.74% females; 26.46% males) and individuals older than 1 year ( > 1-yo hereafter: 58.26% females; 73.54% males) that we included as a control variable in our capture-recapture models. This classification allowed the inclusion of non-local breeders (immigrants) in our analyses.Models were built and fitted to the data using E-SURGE 2.2.048. As our data were annually collected and we had no data for 2003, we selected the “Unequal Time Intervals” option to account for the 2002–2004 interval. Details on the probabilistic framework and the limitations of the modelling approach are given in Appendix S4.Goodness of fitBefore running the capture-recapture analysis, we preliminary assessed the goodness of fit (GOF) of a general model to the data. Since GOF tests are not available for multi-event models, we tested the GOF of the Cormack-Jolly-Seber (CJS), a model accounting for just two states, alive and dead, and for temporal variation in survival (Transition) and recapture (Event) probabilities, using U-CARE 2.3.249. This approach is conservative because the CJS is coarser than the MECR model. Thus, if the former fits the data well, the latter will fit them. All the GOF tests were run for males and females separately. The global tests were not significant for both males [c2 = 72.57, df = 103, p = 0.99; N(0,1) statistic for transient ( > 0) =  − 0.49, p = 0.69; N(0,1) signed statistic for trap-dependence = − 0.84, p = 0.99] and females [c2 = 76.13, df = 122, p = 0.99; N(0,1) statistic for transient ( > 0) = − 2.51, p = 0.69; N(0,1) signed statistic for trap-dependence = − 1.22, p = 0.22], indicating acceptable fits of the Cormack-Jolly-Seber models to the data. For the complete results of 3.SR (transience) and 2.CT (trap-dependence) tests, see Appendix S5.Model selectionModel selection was based on Akaike Information Criterion corrected for small sample sizes (AICc)50. For each sex, in a preliminary analysis, we built a global model checking that there were no parameter identifiability issues48. The structure of the global model was: Initial State (mating status × time), Local survival (age + (mating status × time)), Mating Status Change (age × mating status), Recapture (mating status × time), Mating status Assignment (mating status × time).Our modelling approach consisted of two steps. In step one, starting from the global model, we followed a backwards model selection procedure to test various combinations of variables potentially influencing each parameter of the MECR model while simplifying the model’s structure. According to the classic approach for which the recapture part of the model is modelled before that of survival51,52, we followed the following order of model selection: Initial State, Mating Status Assignment, Recapture, Mating Status Change, and Local survival. After testing the model structure (set of effects) for a parameter, we set the best structure (lower AICc) for that parameter, and we then tested the models for the following parameter. Thus, at the end of step one, we examined the effect of mating status on the biologically relevant parameters, that is, on Local survival and Mating Status Change. In step two, we used the simplified model resulting from step one (final model 1) to test whether the frequency of the FSP differentially affected the biologically relevant parameters according to the mating status. First, we tested the effects of FSP on Mating Status Change and then on Local survival (by keeping in MSC the same structure of final model 1). In the Results section, we reported parameter’ estimates from a model that combined the best final structure (lowest AICc) found on all the parameters, when not stated otherwise.Linear regression analysis of FSP and fledging success of hatchlingsWe used a GLM model to test whether the FSP depends on the yearly average proportion of hatchlings that fledged. We used the simulateResiduals function of the DHARMa53 package in R54 to confirm the absence of over-dispersion and the good fit of the model. More

Surface active agents (surfactants) are a group of chemicals that have a polar hydrophilic headgroup and a non-polar lipophilic hydrocarbon tail group1. This structure of surfactants allows them to be used in households and industries to increase the solubility of non-water-soluble substances, such as cleaning agents and emulsifiers. Global production of synthetic surfactants was 7.2 million tons in 20002; since 2006, this value has risen to 12.5 million tons3 and these numbers will grow with the growth of the detergent and cosmetics industry. After use, the residual surfactants are discharged into the sewage system or directly into surface water, resulting in an increase in the level of surfactants in the environment and a significant impact on the ecosystem1.The toxicity of surfactants to organisms is well known4 and depends on the physico-chemical properties of the surfactants themselves. They are generally classified into anionic, cationic, amphoteric and nonionic, depending on the charge of their headgroup. Among the groups listed above, the anionic surfactants are the most common in everyday and industrial uses and are toxic to both humans and the environment. In particular, anionic surfactants can bind to peptides, enzymes and DNA and alter their spatial layout (folding) and surface charge5. Such interactions can change the biological functions of biomolecules. Sodium dodecyl sulfate (SDS) is one of the most commonly used anionic surfactants, producing more than 3.8 million tons globally for industrial applications in cosmetics, clothing, food, fuel, and medicine6. Such mass production and use of SDS results in releases to the environment, with a semi-lethal concentration of not more than 45 μg/ml7 for algae, fish and crustaceans. In addition, it is known that surfactants can accumulate in the human body and cause autoimmune diseases, brain, liver, kidney and lung damage8,9. Besides the permissible limits for surfactants is 1 mg/l in water and at 0.5 mg/l for potable water10. In order to prevent negative environmental and human impacts of anionic surfactants (in particular SDS) in a timely manner, methods are needed to detect this surfactant in both wastewater and surface waters and in the soil11, food9, dust12,13, etc.Spectrophotometric and potentiometric methods are the most common means of determining anionic surfactant, and chromatography is often used to concentrate and separate complex surfactant mixtures14. Most often, the ionometric determination of the surfactant is carried out using ionic electrodes, which makes it possible to determine the concentration of the substance under investigation in a short time (up to 30 min). However, this method has low sensitivity (280–600 μg/ml)15 and low selectivity, which does not allow the determination of surfactants in relatively complex samples. Spectrophotometric methods are also labour-free (10–30 min) and have a high sensitivity of 0.001 μg/ml16,17. The main disadvantage is the low specificity and dilution of the sample to the measuring limit of 0.01 μg/ml, which complicates the measurement procedure. These defects are corrected by chromatography, which allows separating the studied mixture and increasing the concentration of the required substance, but this procedure requires a minimum of several hours15.There is therefore a need to develop a fast, low-cost method for determining anionic surfactant with high selectivity (specificity) that allows measurements to be made at environmentally toxic concentrations (10–50 μg/ml). Therefore, a quick semi-quantitative or qualitative determination of the substance by means of various rapid tests, such as paper tests, is sufficient for a number of practical tasks to determine the surfactants before applying a more precise and labour-intensive method systems, tracer powders, fabrics, polymer films, tablets18,19,20,21. In particular, Dmitrienko’s work with co-authors presents a method based on adsorption of a red-colored polyurethane foam (PUF) complex of an anionic surfactant with cation 1,10-fenantrolinate iron complex(II)22. This method allows the determination of anionic surfactants between 1 and 30 μg/ml. But all these systems have a common disadvantage—the need to use toxic reagents.Thus, we propose a non-toxic diagnostic system based on polyelectrolyte microcapsules for quick, cheap and highly selective qualitative and semi-quantitative determination of SDS in the medium. More

1.Askew, A. P., Corker, D., Hodkinson, J. & Thompson, K. A new apparatus to measure the rate of fall of seeds. Funct. Ecol. 11, 121–125. https://doi.org/10.1046/j.1365-2435.1997.00049.x (1997).Article 

Google Scholar 
2.Bohrer, G., Katul, G. G., Nathan, R., Walko, R. L. & Avissar, R. Effects of canopy heterogeneity, seed abscission and inertia on wind-driven dispersal kernels of tree seeds. J. Ecol. 96, 569–580. https://doi.org/10.1111/j.1365-2745.2008.01368.x (2008).Article 

Google Scholar 
3.Nathan, R. et al. Mechanistic models of seed dispersal by wind. Thyroid Res. 4, 113–132. https://doi.org/10.1007/s12080-011-0115-3 (2011).ADS 
Article 

Google Scholar 
4.Abedi-Lartey, M., Dechmann, D. K. N., Wikelski, M., Scharf, A. K. & Fahr, J. Long-distance seed dispersal by straw-coloured fruit bats varies by season and landscape. Global Ecol. Conserv. 7, 12–24. https://doi.org/10.1016/j.gecco.2016.03.005 (2016).Article 

Google Scholar 
5.Garcıa, C., Klein, E. K. & Jordano, P. Dispersal processes driving plant movement: Challenges for understanding and predicting range shifts in a changing world. J. Ecol. 105, 1–5. https://doi.org/10.1111/1365-2745.12705 (2017).Article 

Google Scholar 
6.Sutherland, W. J. et al. Identification of 100 fundamental ecological questions. J. Ecol. 101, 58–67. https://doi.org/10.1111/1365-2745.12025 (2013).Article 

Google Scholar 
7.Beckman, N. G., Aslan, C. E. & Rogers, H. S. The role of seed dispersal in plant populations: Perspectives and advances in a changing world. AoB Plants 12, plaa010. https://doi.org/10.1093/aobpla/plaa010 (2020).Article 
PubMed 
PubMed Central 

Google Scholar 
8.Willson, M. F., Rice, B. L. & Westoby, M. Seed dispersal spectra: a comparison of temperate plant-communities. J. Veg. Sci. 1, 547–562 (1990).Article 

Google Scholar 
9.van Rheede, K., Oudtshoorn, V. & van Rooyen, M. W. Dispersal Biology of Desert Plants (Springer, 1999).
Google Scholar 
10.Clark, J. S., Silman, M., Kern, R. & Hillerislambers, M. J. Seed dispersal near and far: patterns across temperate and tropical forests. Ecology 80, 1475–1494. https://doi.org/10.2307/176541 (1999).Article 

Google Scholar 
11.Maler, A., Emig, W. & Leins, P. Dispersal patterns of some Phyteuma species (Campanulaceae). Plant Biol. 1, 408–417. https://doi.org/10.1111/j.1438-8677.1999.tb00723.x (1999).Article 

Google Scholar 
12.Baker, D. V. & Beck, G. The weed tunnel: building an experimental wind tunnel. Weed Technol. 22, 549–552. https://doi.org/10.1614/wt-07-162.1 (2008).Article 

Google Scholar 
13.Planchuelo, G., Catalán, P. & Delgado, J. A. Gone with the wind and the stream: Dispersal in the invasive species Ailanthus altissima. Acta Oecol. 73, 31–37. https://doi.org/10.1016/j.actao.2016.02.006 (2016).ADS 
Article 

Google Scholar 
14.Liang, W. et al. How do diaspore traits, wind speed and sand surface configuration interact to determine seed burial during wind dispersal?. Plant Soil 440, 357–368. https://doi.org/10.1007/s11104-019-04071-4 (2019).CAS 
Article 

Google Scholar 
15.Zhou, Q. et al. Responses of secondary wind dispersal to environmental characteristics and diaspore morphology of seven Calligonum species. Land Degrad. Dev. 31, 842–850. https://doi.org/10.1002/ldr.3489 (2020).Article 

Google Scholar 
16.Andersen, M. C. An analysis of variability in seed settling velocities of several wind-dispersed Asteraceae. Am. J. Bot. https://doi.org/10.1002/j.1537-2197.1992.tb13702.x (1992).Article 
PubMed 

Google Scholar 
17.Augspurger, C. K., Franson, S. E., Cushman, K. C. & Muller-Landau, H. C. Intraspecific variation in seed dispersal of a Neotropical tree and its relationship to fruit and tree traits. Ecol. Evol. 6, 1128–1142. https://doi.org/10.1002/ece3.1905 (2016).Article 
PubMed 
PubMed Central 

Google Scholar 
18.Greene, D. F. & Quesada, M. The differential effect of updrafts, downdrafts and horizontal winds on the seed abscission of Tragopogon dubius. Funct. Ecol. 25, 468–472. https://doi.org/10.1111/j.1365-2435.2010.01788.x (2011).Article 

Google Scholar 
19.Paice, M., Day, W., Rew, J. & Howard, A. A simulation model for evaluating the concept of patch spraying. Weed Res. 38, 373–388 (1998).Article 

Google Scholar 
20.Jongejans, E. & Telenius, A. Field experiments on seed dispersal by wind in ten umbelliferous species (Apiaceae). Plant Ecol. 152, 67–78. https://doi.org/10.1023/a:1011467604469 (2001).Article 

Google Scholar 
21.Augspurger, C. K. Morphology and dispersal potential of wind-dispersed diaspore of neotropical trees. Am. J. Bot. 73, 353–363 (1986).Article 

Google Scholar 
22.Zhu, J., Liu, M., Xin, Z., Liu, Z. & Schurr, F. M. A trade-off between primary and secondary seed dispersal by wind. Plant Ecol. 220, 541–552. https://doi.org/10.1007/s11258-019-00934-z (2019).Article 

Google Scholar 
23.Phartyal, S. S., Rosbakh, S., Ritz, C., Poschlod, P. & Bruun, H. H. Ready for change: Seed traits contribute to the high adaptability of mudflat species to their unpredictable habitat. J. Veg. Sci. 31, 331–342. https://doi.org/10.1111/jvs.12841 (2020).Article 

Google Scholar 
24.Saatkamp, A. et al. A research agenda for seed-trait functional ecology. New Phytol. 221, 1764–1775. https://doi.org/10.1111/nph.15502 (2019).Article 
PubMed 

Google Scholar 
25.Green, D. S. The terminal velocity and dispersal of spinning samaras. Am. J. Bot. 67, 1218–1224 (1980).Article 

Google Scholar 
26.Jongejans, E., Pedatella, N. M., Shea, K., Skarpaas, O. & Auhl, R. Seed release by invasive thistles: the impact of plant and environmental factors. Proc R Soc B Biol Sci 274, 2457–2464. https://doi.org/10.1098/rspb.2007.0190 (2007).Article 

Google Scholar 
27.Andersen, M. Mechanistic models for the seed shadows of wind-dispersed plants. Am. Nat. 137, 476–497. https://doi.org/10.1086/285178 (1991).Article 

Google Scholar 
28.Cousens, R. D. & Rawlinson, A. A. When will plant morphology affect the shape of a seed dispersal “kernel”?. J Thenr Biol. 211, 229–238. https://doi.org/10.1006/jtbi.2001.2341 (2001).ADS 
CAS 
Article 

Google Scholar 
29.Qin, X. et al. Shrub canopy interception of diaspores dispersed by wind. Seed Sci. Res. 30, 310–318. https://doi.org/10.1017/S0960258520000410 (2020).Article 

Google Scholar 
30.Kelly, N., Cousens, R. D., Taghizadeh, M. S. & Hanan, J. S. Plants as populations of release sites for seed dispersal: A structural-statistical analysis of the effects of competition on Raphanus raphanistrum. J. Ecol. 101, 878–888. https://doi.org/10.1111/1365-2745.12097 (2013).Article 

Google Scholar 
31.Donohue, K. Maternal determinants of seed dispersal in cakile edentula: Fruit, plant, and site traits. Ecology 79, 2771–2788 (1998).Article 

Google Scholar 
32.Bullock, J. M. & Moy, I. L. Plants as seed traps: inter-specific interference with dispersal. Acta Oecol. 25, 35–41. https://doi.org/10.1016/j.actao.2003.10.005 (2004).ADS 
Article 

Google Scholar 
33.Cousens, R., Dytham, C. & Law, R. Dispersal in plants. A population perspective. Ann. Bot. https://doi.org/10.1111/j.1442-9993.2010.02216.x (2008).Article 

Google Scholar 
34.Phillips, M. L. Dispersal in plants: A population perspective. Austral Ecol. 36, e27 (2011).
Google Scholar 
35.Pounden, E., Greene, D. F., Quesada, M. & Contreras Sánchez, J. M. The effect of collisions with vegetation elements on the dispersal of winged and plumed seed. J. Ecol. 96, 591–598 (2008).Article 

Google Scholar 
36.Poggi, D., Porporato, A., Ridolfi, L., Albertson, J. D. & Katul, G. G. The effect of vegetation density on canopy sub-layer turbulence. BoLMe 111, 565–587 (2004).ADS 

Google Scholar 
37.Greene, D. F. The role of abscission in long-distance seed dispersal by the wind. Ecology 86, 3105–3110 (2005).Article 

Google Scholar 
38.Matlack, G. R. Diaspore size, shape, and fall behavior in wind-dispersed plant species. Am. J. Bot. https://doi.org/10.1002/j.1537-2197.1987.tb08729.x (1987).Article 

Google Scholar 
39.Zhou, Q. et al. Relationship between seed morphological traits and wind dispersal trajectory. Funct. Plant Biol. 46, 1063–1071. https://doi.org/10.1071/FP19087 (2019).Article 
PubMed 

Google Scholar 
40.Heydel, F., Cunze, S., Bernhardt-Römermann, M. & Tackenberg, O. Long-distance seed dispersal by wind: Disentangling the effects of species traits, vegetation types, vertical turbulence and wind speed. Ecol. Res. 29, 641–651. https://doi.org/10.1007/s11284-014-1142-5 (2014).Article 

Google Scholar 
41.Greene, D. F. & Johnson, E. A. Long-distance wind dispersal of tree seeds. Can J Bot. 73, 1036–1045 (1995).Article 

Google Scholar 
42.Skarpaas, O., Auhl, R. & Shea, K. Environmental variability and the initiation of dispersal: turbulence strongly increases seed release. Proc. R. Soc. B Biol. Sci. 273, 751–756. https://doi.org/10.1098/rspb.2005.3366 (2006).Article 

Google Scholar 
43.Lipoma, M. L., Cuchietti, A., Gorne, L. D., Díaz, S. M. & Zobel, M. Not gone with the wind: Vegetation complexity increases seed retention during windy periods in the Argentine Semiarid Chaco. J. Veg. Sci. 30, 542–552. https://doi.org/10.1111/jvs.12747 (2019).Article 

Google Scholar 
44.Soons, M. B., Heil, G. W., Ran, N. & Katul, G. G. Determinants of long-distance seed dispersal by wind in grasslands. Ecology 85, 3056–3068. https://doi.org/10.1890/03-0522 (2004).Article 

Google Scholar 
45.Soons, M. B. & Bullock, J. M. Non-random seed abscission, long-distance wind dispersal and plant migration rates. J. Ecol. 96, 581–590. https://doi.org/10.1111/j.1365-2745.2008.01370.x (2008).Article 

Google Scholar 
46.Rotundo, J. L. & Aguiar, M. R. Litter effects on plant regeneration in arid lands: A complex balance between seed retention, seed longevity and soil–seed contact. J. Ecol. 93, 829–838. https://doi.org/10.1111/j.1365-2745.2005.01022.x (2005).Article 

Google Scholar 
47.Tackenberg, O., Poschlod, P. & Bonn, S. Assessment of wind dispersal potential in plant species. Ecol. Monogr. 73, 191–205 (2003).Article 

Google Scholar 
48.Schurr, F. M., Bond, W. J., Midgley, G. F. & Higgins, S. I. A mechanistic model for secondary seed dispersal by wind and its experimental validation. J. Ecol. 93, 1017–1028. https://doi.org/10.1111/j.1365-2745.2005.01018.x (2005).Article 

Google Scholar 
49.McGinley, M. & Brigham, E. Fruit morphology and terminal velocity in Tragopogon dubious (L). Funct. Ecol. 4, 489–496 (1989).Article 

Google Scholar  More

Current and past scenariosThe Arbutus unedo distribution presence database and eight environmental variables were used to predict the species potential distribution under past, current, and future climatic conditions. The current fitted model predicted the A. unedo suitable habitat with good performance, when evaluated using AUC (AUC = 0.92)31.The results suggest that the species distribution was mainly determined by the following attributes: (1) mean temperature of the driest quarter (BIO9), (2) annual mean temperature (BIO1), (3) slope and (4) temperature annual range (BIO7), summing up to 78.6% of the model’s total contributions (Table 1). In a recent study for the same species in Portugal, the authors considered that the variables precipitation seasonality (BIO15) and slope had a net significant influence on the species’ habitat suitability16. The variable slope was expected to influence the A. unedo distribution, since this tree is ecologically adapted to rocky slopes, which are edaphically unsuitable for woodland, where the vegetation is limited to shrub communities. The Arbutus unedo is particularly intolerant to shade, and this high light requirement is particularly needed for fruit production, restricting this species to settlement with open habitats21. This aspect was confirmed in a study about a woodland disturbance, where the A. unedo widespread initially, and, afterward, the species declined as canopy woodland re-developed35.The Arbutus unedo was adapted to a wide range of climatic conditions considering the response curves analysis (Fig. 1). According to Santiso24, this species developed a conservative strategy in the use of nutrients and water, when they are scarce. The Arbutus unedo plasticity and evolvability explained the current species persistence throughout its distribution range, which will be decisive in future response to climate change. The species distribution prediction for the present is in agreement with the species distribution according to Caudullo et al.38 (Fig. 6). Indeed, the A. unedo is widespread in Portugal, Galicia, and southwest of Spain, occupying the coastal belt from Tunisia to Morocco along the North of Africa, and from the Iberian Peninsula to Turkey across southern Europe, yet occurring in the northern Iberian Peninsula, western France, and south-western Ireland. The results showed that MaxEnt fitted model suitability prediction was high in those regions, which was further supported by the high observed AUC. Furthermore, it is important to stress the current potential distribution meaningful accuracy, and, therefore, the ability to extrapolate past and future prediction scenarios for A. unedo distribution in the Mediterranean Basin. Not surprisingly, the highest predicted suitability for the species occurs mostly in its native range, in regions where the mean January temperature is above 4ºC, a limit required for the species’ survival23. Unfortunately, this variable is only available in the WorldClim database for the present, but not for the future or past periods. Nevertheless, the species’ current spatial distribution was overlapped to the areas with mean January temperature higher than 4ºC, and a high overlap is reckoned, thus confirming this factor ecological importance17,23. The analysis of the species’ predicted suitability show that high-suitability areas [0.8–1.0] are more frequent in current conditions, compared to other scenarios, and thus confirming the suitability of climate current conditions for the species’ development (Fig. 2d).Figure 6Arbutus unedo L. distribution area and occurrences. Black dots represent the species’ occurrences (Supplementary Table S1), and the green area represents the species native distribution range for A. unedo downloaded from the data sets in38. Equal-area projection EPSG:3035. The map was made using the sf 1.0–3 and the terra 1.4–14 (https://rspatial.org/terra/) packages in R 4.1.032,33.Full size imageThe predicted habitat suitability for the LGM and MH scenarios differed, as expected, from the obtained for the current climatic conditions. The climatic conditions in the MH, were more suitable than during the LGM, due to climate warming, thus allowing the species expansion. According to the LGM projections, the species occurred in the Mediterranean islands (Sicily, Sardinia, and Balearic), North of Italy, south and eastern coast of Spain, Catalonia, North Africa, and spots in Portugal, Greece, and Turkey (Fig. 2a). These were suitable putative cryptic refugia areas that remained during the LGM. Keppel et al.39 defined refugia as locations to which species retreated during periods of adverse climate, and could potentially expand from, when environmental conditions turn out to be more suitable for the species. These locations/habitats were responsible for the species’ survival under changing environmental conditions for millennia. Identifying and characterizing climate refugia provides an important context for understanding the modern species distribution development, traits, and local adaptation40. The Iberian, Italian and Balkan peninsulas, which remained relatively ice-free during the ice ages were identified as the main glacial tree refugia areas in Europe41,42.The Betic mountains, in the Iberian Peninsula (IP), were sought to be a species’ refuge from the last ice age, nearby the sea (Valencia region, Spain), where fossil pollen pieces of evidence were found in the Canal de Navarrés peat deposit4, with 34 ka, and in the Siles lake (ca. 19 ka) located in the Segura mountains of southern Spain43 (Fig. 5, numbers 1 and 4, respectively). These mountains were ice-free during the Late Pleistocene, and the persistence of a mild climate could explain that other thermophilic species have prevailed in this particular region during the Full Glacial, proved by high genetic diversity levels44, and with the pollen spectra suggesting the region as a glacial refugium for temperate and Mediterranean trees9,13,43.The Iberian Peninsula border could have been a refuge for the species, due to the importance of the sea as a temperature regulator. Indeed, pollen evidence was found in the Basque mountains (northern IP, Fig. 5, number 3: 21 ka)45. Despite the information given by the MaxEnt modelling, the A. unedo reconstruction for the LGM was observed in other areas, since paleo-evidences were found in unexpected northern sites (Donatella Magri, personal communication). Certainly, fossil charcoal pieces of evidence were found in the central region of Portugal (Serra do Sicó), nearby the sea, supporting the presence of thermophilous taxa, including A. unedo, during the Full Glacial (24 ka, Fig. 5, number 2), and the presence of other Mediterranean taxa, such as olive tree and Pistacia lentiscus L., suggested that this may have been refuge zones for thermophilous plants46. Corroboration of the species expansion later during the warm and humid period of the Holocene (see Fig. 3b; green area) was confirmed by pollen remains (Fig. 5), particularly inland and northwardly, in the IP, and, also, in a molecular population genetics study made in Portugal, with northward haplotype migration47.Santiso et al.17 in a study about the A. unedo phylogeography, concluded that this species had two clades, separated during the late Pleistocene, before the LGM, suggesting that it may have coincided with the hardest glaciations recorded in the Quaternary. One clade occupied the Atlantic Iberia and, possibly, North Africa, while the other occurred in the western Mediterranean Basin in Spain. Besides, the A. unedo possibly persisted in the late Quaternary in the western Mediterranean, based on chloroplast DNA observations, and the results from the current study supported this interpretation, and by the fossil evidence (Supplementary Table 3). The same strong genetic differentiation between the western and eastern Mediterranean Basin was also found in olive lineages, and this pattern is congruent in other Mediterranean shrubs and tree species9,13,15. In another species (maritime pine), Bucci et al.44 confirmed the existence of a genetic divide between eastern and western lineages, also previously described by Burban and Petit48 based on mitochondrial DNA. Additionally, maximum haplotypic diversity for this species was found in south-eastern and central Spain, which, therefore, may be considered a biodiversity hotspot and a strong signal for refuge, as long-term populations tend to harbour more diversity than recently expanding ones41. In their study with A. unedo in Portugal using cpSSR, Ribeiro et al.47 found signals of two putative refugia in southern and central littoral in the country, also supported by macrofossil and pollen remains. Furthermore, according to Médail and Diadema41, several regions in the Western Mediterranean (large Mediterranean islands, North Africa and Catalonia), could have played a role in the case of A. unedo, and sought as refugia locations.The Arbutus unedo suitability maps obtained for past conditions support the claim that MH climatic conditions were more suitable for the species expansion than the LGM ones (Fig. 2a,b, and d). Consequently, spatial distribution changes analysis between past and current scenarios (Fig. 3) showed that A. unedo potential distribution expanded extensively from the LGM-current and the MH-current periods, indicating that more suitable areas were available for the species at present. Pollen records during the Early-Middle Holocene supported the Mediterranean presence of this species, in particular, nearby the sea (Fig. 5), which is under the MaxEnt predicted eastward dispersion (Fig. 3a, b). In another species bird-dispersed, the Myrtus communis L. (myrtle), a spread from west to east was verified, following genetic differentiation during the Pleistocene15, since the western region of the Mediterranean Basin had a milder climate, compared to the eastern one, during the LGM.The Arbutus unedo had, probably, a considerable ability to disperse, migrating over thousands of kilometres and even crossing sea stretches, allowing the species expansion during the MH, when the climatic conditions became more favourable17. In the future, A. unedo migration possibility will depend on climatic change pace and, also, on seed dispersal, particularly long-distance dispersal events49, since long-distance migration fitness will be useful in future change scenarios, allowing species to progress to newly suitable areas, as happened in the past. Nevertheless, from past evidence in Ireland, Britain, and across Europe, trees migrated faster than would be expected, during the warming period at the beginning of the Holocene, and that biotic and/or abiotic vectors must have been involved in dispersal50. This is the case with A. unedo, the seeds are dispersed by different type of birds that eat the soft fruit and, some of them are migratory species, as the European robin (Erithacus rubecula L.)51.Future scenariosThe impact of climate change on the A. unedo potential distribution was assessed under two representative concentration pathways (RCP 4.5 and RCP8.5) for the years 2050 and 2070. The results showed that climate change, under both moderate and high emission scenarios, will affect the species distribution range. A decrease in the predicted suitable areas was, generally, observed, since the climate becomes less favourable for the species in the future.According to the results, under the RCP 4.5 scenario, the potential distribution area will increase up to 2050, and decrease afterward according to the prediction for 2070 (Fig. 2d, e, f). Considering the RCP 8.5 scenario, the suitable area will exhibit a net decrease from the present up to 2050 and continue this trend until 2070. It is also expected that he medium–high suitable areas will gradually decrease from the present to the future, especially for scenario RCP 8.5. These results are in agreement with those obtained by other authors, confirming that Mediterranean species (e.g. Quercus sp.) distribution areas will be negatively affected by future climate change52,53.Several studies concluded that global warming will influence species distributions by causing expansions, contractions, or shifts in the species ranges2. As expected, the results from the current study showed that suitable areas will contract under future climate scenarios when compared to the current conditions, though suitable areas will emerge. These effects’ impact will depend on the climate change scenario severity. Despite contraction and expansion effects, the presence of the species will gradually decrease from the 2050s to the 2070s. Moreover, a species’ shift toward the North will be verified, because of suitable areas emergence, observed mainly in the RCP 8.5 scenario. Those areas will mostly emerge in the North of France, South of the United Kingdom, and Ireland, implying species’ latitudinal migration. The species’ northward displacement is consistent with climate change studies’ results obtained by several authors, including A. unedo16,20,54,55. Nevertheless, according to Gerassis et al.20, under climate change, the expected habitat disruption and fragmentation could lead to very adverse conditions for A. unedo survival in the future, which could undoubtedly conduct to a possible species’ presence decline in most of the current distribution area. Moreover, distribution models that predict climate‐induced range shifts do not account for spatial dispersal variation56, but adaptive dispersal evolution always reduced neutral genetic diversity across the species’ range. This means that the species’ genetic pool might be erased, depending on the climate change velocity, amongst other conditions, like landscape fragmentation and competition with other species/crops16,57. Additionally, the species ability to migrate mainly through seeds, dispersed by migratory frugivorous birds, bird abundance, and the velocity of climate change, are key issues for future species survival23,58.These results suggested that future changes in environmental conditions may lead to suitable habitat loss in areas where the species had persisted and with a possible range shift towards the North. These findings also revealed that with continuous future climate warmth, the current potential distribution A. unedo areas will become unsuitable or contract, leading to significant changes in the species’ current distribution pattern and putative presence loss. The possibility of species’ migration will ultimately depend on its capability to keep pace with the changing conditions and the velocity to adapt to environmental changes, such as those presented by habitat and climate modifications56.The Mediterranean Basin is one of the most vulnerable climate change hotspots in the world, thus understanding how future climate changes will disturb Mediterranean plant species distribution will be key for tree management planning and conservation design. Nevertheless, further investigation is needed for species well adapted in this region to assess the impacts of climate change in their current and future potential distributions. Those studies including past climate impact on species distribution should be complemented with phylogeographic methods and paleoclimate reconstructions to locate refuges. Other species distribution models, besides MaxEnt, could be tested, although the lack of absence records data limits considerably the modelling approach. Additionally, the human influence magnitude on the predicted ecological niche should be further studied in detail, including urbanization, industrialization and putative tourism pressure, particularly in coastal areas. More

1.Kelly, J. F. Stable isotopes of carbon and nitrogen in the study of avian and mammalian trophic ecology. Can. J. Zool. 78(1), 1–27 (2000).Article 

Google Scholar 
2.Barnes, C., Sweeting, C. J., Jennings, S., Barry, J. & TandPolunin, N. V. Effect of temperature and ration size on carbon and nitrogen stable isotope trophic fractionation. Funct. Ecol. 21(2), 356–362. https://doi.org/10.1111/j.1365-2435.2006.01224.x (2007).Article 

Google Scholar 
3.Chisholm, B. S., Nelson, D. E. & Schwarcz, H. P. Stable-carbon isotope ratios as a measure of marine versus terrestrial protein in ancient diets. Science 216(4550), 1131–1132. https://doi.org/10.1126/science.216.4550.1131 (1982).CAS 
Article 
PubMed 
ADS 

Google Scholar 
4.O’Leary, M. H. Carbon isotopes in photosynthesis. Bioscience 38(5), 328–336. https://doi.org/10.2307/1310735 (1988).CAS 
Article 

Google Scholar 
5.DeNiro, M. J. & Epstein, S. Influence of diet on the distribution of carbon isotopes in animals. Geochim. Cosmochim. Acta 42(5), 495–506. https://doi.org/10.1016/0016-7037(78)90199-0 (1978).CAS 
Article 
ADS 

Google Scholar 
6.DeNiro, M. J. & Epstein, S. Influence of diet on the distribution of nitrogen isotopes in animals. Geochim. Cosmochim. Acta 45(3), 341–351. https://doi.org/10.1016/0016-7037(81)90244-1 (1981).CAS 
Article 
ADS 

Google Scholar 
7.Post, D. M. Using stable isotopes to estimate trophic position: Models, methods, and assumptions. Ecology 83(3), 703–718. https://doi.org/10.1890/0012-9658(2002)083[0703:USITET]2.0.CO;2 (2002).Article 

Google Scholar 
8.McCutchan, J. H. Jr., Lewis, W. M. Jr., Kendall, C. & McGrath, C. C. Variation in trophic shift for stable isotope ratios of carbon, nitrogen, and sulfur. Oikos 102(2), 378–390. https://doi.org/10.1034/j.1600-0706.2003.12098.x (2003).CAS 
Article 

Google Scholar 
9.Tieszen, L. L., Boutton, T. W., Tesdahl, K. G. & Slade, N. A. Fractionation and turnover of stable carbon isotopes in animal tissues: Implications for δ 13 C analysis of diet. Oecologia 57(1–2), 32–37. https://doi.org/10.1007/BF00379558 (1983).CAS 
Article 
PubMed 
ADS 

Google Scholar 
10.Casey, M. M. & Post, D. M. The problem of isotopic baseline: Reconstructing the diet and trophic position of fossil animals. Earth Sci. Rev. 106(1–2), 131–148. https://doi.org/10.1016/j.earscirev.2011.02.001 (2011).CAS 
Article 
ADS 

Google Scholar 
11.West, J. B. et al. (eds) Isoscapes: Understanding Movement, Pattern, and Process on Earth Through Isotope Mapping (Springer, 2009).
Google Scholar 
12.Cheeseman, A. W. & Cernusak, L. A. Isoscapes: A new dimension in community ecology. Tree Physiol. 36(12), 1456–1459. https://doi.org/10.1093/treephys/tpw099 (2016).Article 

Google Scholar 
13.Hellmann, C., Rascher, K. G., Oldeland, J. & Werner, C. Isoscapes resolve species-specific spatial patterns in plant–plant interactions in an invaded Mediterranean dune ecosystem. Tree Physiol. 36(12), 1460–1470. https://doi.org/10.1093/treephys/tpw075 (2016).CAS 
Article 
PubMed 

Google Scholar 
14.Chiocchini, F., Portarena, S., Ciolfi, M., Brugnoli, E. & Lauteri, M. Isoscapes of carbon and oxygen stable isotope compositions in tracing authenticity and geographical origin of Italian extra-virgin olive oils. Food Chem. 202, 291–301. https://doi.org/10.1016/j.foodchem.2016.01.146 (2016).CAS 
Article 
PubMed 

Google Scholar 
15.Newton, J. An insect isoscape of UK and Ireland. Rapid Commu. Mass Spectrom. 1, e9126 (2021).
Google Scholar 
16.Veen, T. et al. Identifying the African wintering grounds of hybrid flycatchers using a multi–isotope (δ 2 H, δ 13 C, δ 15 N) assignment approach. PLoS ONE 9(5), e98075 (2014).Article 
ADS 

Google Scholar 
17.Schneider, K. et al. Trophic niche differentiation in soil microarthropods (Oribatida, Acari): Evidence from stable isotope ratios (15N/14N). Soil Biol. Biochem. 36(11), 1769–1774. https://doi.org/10.1016/j.soilbio.2004.04.033 (2004).CAS 
Article 

Google Scholar 
18.Menichetti, L. et al. Increase in soil stable carbon isotope ratio relates to loss of organic carbon: Results from five long-term bare fallow experiments. Oecologia 177(3), 811–821. https://doi.org/10.1007/s00442-014-3114-4 (2015).Article 
PubMed 
ADS 

Google Scholar 
19.Amundson, R. et al. Global patterns of the isotopic composition of soil and plant nitrogen. Glob. Biogeochem. Cycles https://doi.org/10.1029/2002GB001903 (2003).Article 

Google Scholar 
20.Craine, J. M. et al. Ecological interpretations of nitrogen isotope ratios of terrestrial plants and soils. Plant Soil 396(1–2), 1–26. https://doi.org/10.1007/s11104-015-2542-1 (2015).CAS 
Article 

Google Scholar 
21.Ben-David, M. & Flaherty, E. A. Stable isotopes in mammalian research: A beginner’s guide. J. Mammal. 93(2), 312–328. https://doi.org/10.1644/11-MAMM-S-166.1 (2012).Article 

Google Scholar 
22.del Rio, C. M. & Carleton, S. A. How fast and how faithful: The dynamics of isotopic incorporation into animal tissues. J. Mammal. 93(2), 353–359. https://doi.org/10.1644/11-MAMM-S-165.1 (2012).Article 

Google Scholar 
23.Clementz, M. T. New insight from old bones: Stable isotope analysis of fossil mammals. J. Mammal. 93(2), 368–380. https://doi.org/10.1644/11-MAMM-S-179.1 (2012).Article 

Google Scholar 
24.Inger, R. et al. Temporal and intrapopulation variation in prey choice of wintering geese determined by stable isotope analysis. J. Anim. Ecol. 75(5), 1190–1200. https://doi.org/10.1111/j.1365-2656.2006.01142.x (2006).Article 
PubMed 

Google Scholar 
25.Jackson, A. L., Inger, R., Parnell, A. C. & Bearhop, S. Comparing isotopic niche widths among and within communities: SIBER–Stable Isotope Bayesian Ellipses in R. J. Anim. Ecol. 80(3), 595–602. https://doi.org/10.1111/j.1365-2656.2011.01806.x (2011).Article 
PubMed 

Google Scholar 
26.Jackson, M. C. et al. Population-level metrics of trophic structure based on stable isotopes and their application to invasion ecology. PLoS ONE https://doi.org/10.1371/journal.pone.0031757 (2012).Article 
PubMed 
PubMed Central 

Google Scholar 
27.Semmens, B. X. et al. Statistical basis and outputs of stable isotope mixing models: Comment on Fry (2013). Mar. Ecol. Prog. Ser. 490, 285–289. https://doi.org/10.3354/meps10535 (2013).Article 
ADS 

Google Scholar 
28.Parnell, A. C., Inger, R., Bearhop, S. & Jackson, A. L. Source partitioning using stable isotopes: Coping with too much variation. PLoS ONE https://doi.org/10.2307/1310735 (1988).Article 

Google Scholar 
29.Phillips, D. L. et al. Best practices for use of stable isotope mixing models in food-web studies. Can. J. Zool. 92(10), 823–835. https://doi.org/10.1139/cjz-2014-0127 (2010).Article 

Google Scholar 
30.Judge, J., Wilson, G. J., Macarthur, R., McDonald, R. A. & Delahay, R. J. Abundance of badgers (Meles meles) in England and Wales. Sci. Rep. 7(1), 1–8. https://doi.org/10.1038/s41598-017-00378-3 (2017).CAS 
Article 

Google Scholar 
31.Allen, A. et al. Genetic evidence further elucidates the history and extent of badger introductions from Great Britain into Ireland. R. Soc. Open Sci. 7(4), 200–288. https://doi.org/10.1098/rsos.200288 (2020).CAS 
Article 

Google Scholar 
32.Davies, J. M., Lachno, D. R. & Roper, T. J. The anal gland secretion of the European badger (Meles meles) and its role in social communication. J. Zool. 216(3), 455–463. https://doi.org/10.1111/j.1469-7998.1988.tb02441.x (1988).CAS 
Article 

Google Scholar 
33.Lüps, P., Roper, T. J. & Stocker, G. Stomach contents of badgers (Meles meles L.) in central Switzerland. Mammalia 51(4), 559–570. https://doi.org/10.1515/mamm.1987.51.4.559 (1987).Article 

Google Scholar 
34.Roper, T. J. The European badger Meles meles: Food specialist or generalist?. J. Zool. 234(3), 437–452. https://doi.org/10.1111/j.1469-7998.1994.tb04858.x (1994).Article 

Google Scholar 
35.Roper, T. J. Badger Meles meles setts–architecture, internal environment and function. Mamm. Rev. 22(1), 43–53. https://doi.org/10.1111/j.1365-2907.1992.tb00118.x (1992).Article 

Google Scholar 
36.Feore, S. & Montgomery, W. I. Habitat effects on the spatial ecology of the European badger (Meles meles). J. Zool. 247(4), 537–549. https://doi.org/10.1111/j.1469-7998.1999.tb01015.x (1999).Article 

Google Scholar 
37.Robertson, A., McDonald, R. A., Delahay, R. J., Kelly, S. D. & Bearhop, S. Individual foraging specialisation in a social mammal: The European badger (Meles meles). Oecologia 176(2), 409–421. https://doi.org/10.1007/s00442-014-3019-2 (2014).Article 
PubMed 
ADS 

Google Scholar 
38.Haussmann, N. S. Soil movement by burrowing mammals: A review comparing excavation size and rate to body mass of excavators. Prog. Phys. Geogr. 41(1), 29–45. https://doi.org/10.1177/0309133316662569 (2017).Article 

Google Scholar 
39.Cabana, G. & Rasmussen, J. B. Comparison of aquatic food chains using nitrogen isotopes. Proc. Acad. Natl. Sci. 93(20), 10844–10847. https://doi.org/10.1073/pnas.93.20.10844 (1996).CAS 
Article 
ADS 

Google Scholar 
40.Phillips, D. L. & Gregg, J. W. Source partitioning using stable isotopes: Coping with too many sources. Oecologia 136(2), 261–269. https://doi.org/10.1007/s00442-003-1218-3 (2003).Article 
PubMed 
ADS 

Google Scholar 
41.Wright, D. M. et al. Herd-level bovine tuberculosis risk factors: Assessing the role of low-level badger population disturbance. Sci. Rep. 5(1), 1–11. https://doi.org/10.1038/srep13062 (2015).CAS 
Article 

Google Scholar 
42.Britain, G. The strategy for achieving officially bovine tuberculosis free status for England. Department for Environment, Food & Rural Affairs. https://www.gov.uk/government/publications/a-strategy-for-achieving-officially-bovine-tuberculosis-free-status-for-england. (2014).43.Ireland, G. Spending Review 2019 Animal Health: TB Eradication. Economics and Planning Division, Department of Agriculture, Food and the Marine. http://budget.gov.ie/Budgets/2020/Documents/Budget/Animal%20Health%20-%20TB%20Eradication.pdf. (2019).44.Kruuk, H. Spatial organization and territorial behaviour of the European badger Meles meles. J. Zool. 184(1), 1–19. https://doi.org/10.1111/j.1469-7998.1978.tb03262.x (1978).Article 

Google Scholar 
45.Macdonald, D. W., Newman, C. & Buesching, C. D. Badgers in the rural landscape—Conservation paragon or farmland pariah? Lessons from the Wytham Badger Project. Wildl. Conserv. Farmland 2, 65–95 (2015).
Google Scholar 
46.McDonald, J. L., Robertson, A. & Silk, M. J. Wildlife disease ecology from the individual to the population: Insights from a long-term study of a naturally infected European badger population. J. Anim. Ecol. 87(1), 101–112. https://doi.org/10.1111/1365-2656.12743 (2018).Article 
PubMed 

Google Scholar 
47.Rogers, L. M., Cheeseman, C. L., Mallinson, P. J. & Clifton-Hadley, R. The demography of a high-density badger (Meles meles) population in the west of England. J. Zool. 242(4), 705–728. https://doi.org/10.1111/j.1469-7998.1997.tb05821.x (1997).Article 

Google Scholar 
48.Desktop, E. A. Release 10 437–438 (Environmental Systems Research Institute, 2011).
Google Scholar 
49.Kostka, B.I., Landscape ecology, diet composition and energetics of the Eurasian badger (Meles meles). Unpublished PhD thesis, Queen’s University Belfast. https://ethos.bl.uk/OrderDetails.do?uin=uk.bl.ethos.579755. (2012).50.Scheppers, T. L. et al. Estimating social group size of Eurasian badgers Meles meles by genotyping remotely plucked single hairs. Wildl. Biol. 13(2), 195–207. https://doi.org/10.2981/0909-6396(2007)13[195:ESGSOE]2.0.CO;2 (2007).Article 

Google Scholar 
51.Geological Survey Ireland. Tellus Geochemical Survey: Shallow Topsoil Data from the Border and West of Ireland. Department of the Environment, Climate and Communications. https://secure.dccae.gov.ie/GSI_DOWNLOAD/Geochemistry/Reports/Tellus_A_geochemistry_data_report_2020_v1.0.pdf. Accessed 7Jun 2021.52.Smyth, D. Methods used in the Tellus Geochemical Mapping of Northern Ireland. http://nora.nerc.ac.uk/id/eprint/14008. (2007).53.Murray, R, McCann, T. P. & Cooper, A. A land classification and landscape ecological survey of Northern Ireland. Report, University of Ulster, Coleraine (1992).54.Stewart, P. D. & Macdonald, D. W. Age, sex, and condition as predictors of moult and the efficacy of a novel fur-clip technique for individual marking of the European badger (Meles meles). J. Zool. 241(3), 543–550. https://doi.org/10.1111/j.1469-7998.1997.tb04846.x (1997).Article 

Google Scholar 
55.R Core Team. R: A Language and Environment for Statistical Computing. (R Foundation for Statistical Computing, 2021)56.Met Office. UK Daily Temperature Data, Part of the Met Office Integrated Data Archive System (MIDAS). NCAS British Atmospheric Data Centre, (2006). Accessed 2 Sep 2019.57.Mardia, K. V., Kent, J. T. & Bibby, J. M. Multivariate Analysis (Academic Press Inc, 1979).MATH 

Google Scholar 
58.Venables, W. N. & Ripley, B. D. Modern Applied Statistics with S (Springer, 2002).Book 

Google Scholar 
59.Kendon, M., McCarthy, M., Jevrejeva, S., Matthews, A. & Legg, T. State of the UK climate 2017. Int. J. Climatol. 38, 1–35. https://doi.org/10.1139/z99-165 (2018).Article 

Google Scholar 
60.Kassambara, A. & Mundt, F., Package ‘factoextra’. Extract and Visualize the Results of Multivariate Data Analyses, 76. https://cran.microsoft.com/snapshot/2016-11-30/web/packages/factoextra/factoextra.pdf. (2017).61.Funck, J., Bataille, C., Rasic, J. & Wooller, M. A bio-available strontium isoscape for eastern Beringia: A tool for tracking landscape use of Pleistocene megafauna. J. Quat. Sci. 36(1), 76–90. https://doi.org/10.1002/jqs.3262 (2021).Article 

Google Scholar 
62.Reddin, C. J., Bothwell, J. H., O’Connor, N. E. & Harrod, C. The effects of spatial scale and isoscape on consumer isotopic niche width. Funct. Ecol. 32(4), 904–915. https://doi.org/10.1111/1365-2435.13026 (2018).Article 

Google Scholar 
63.Layman, C. A. et al. Applying stable isotopes to examine food-web structure: An overview of analytical tools. Biol. Rev. 87(3), 545–562. https://doi.org/10.1111/j.1469-185X.2011.00208.x (2012).Article 
PubMed 

Google Scholar 
64.Fabrizio, M. et al. Habitat suitability vs landscape connectivity determining roadkill risk at a regional scale: A case study on European badger (Meles meles). Eur. J. Wildl. Res. 65(1), 7. https://doi.org/10.1007/s10344-018-1241-7 (2019).Article 

Google Scholar 
65.Rosalino, L. M. et al. Climate and landscape changes as driving forces for future range shift in southern populations of the European badger. Sci. Rep. 9(1), 1–15. https://doi.org/10.1038/s41598-019-39713-1 (2019).CAS 
Article 

Google Scholar 
66.Potts, J. R., Fagan, W. F. & Mourão, G. Deciding when to intrude on a neighbour: Quantifying behavioural mechanisms for temporary territory expansion. Thyroid Res. 12(3), 307–318. https://doi.org/10.1007/s12080-018-0396-x (2019).Article 

Google Scholar 
67.Noonan, M. J. et al. Knowing me, knowing you: Anal gland secretion of European Badgers (Meles meles) codes for individuality, sex and social group membership. J. Chem. Ecol. 45(10), 823–837. https://doi.org/10.1007/s10886-019-01113-0 (2019).CAS 
Article 
PubMed 

Google Scholar 
68.Kurek, P. Topsoil mixing or fertilization? Forest flora changes in the vicinity of badgers’ (Meles meles L.) setts and latrines. Plant Soil 437(1–2), 327–340. https://doi.org/10.1007/s11104-019-03984-4 (2019).CAS 
Article 

Google Scholar 
69.Abduriyim, S. et al. Variation in pancreatic amylase gene copy number among Eurasian badgers (Carnivora, Mustelidae, Meles) and its relationship to diet. J. Zool. 308(1), 28–36. https://doi.org/10.1111/jzo.12649 (2019).Article 

Google Scholar 
70.Balestrieri, A., Remonti, L., Saino, N. & Raubenheimer, D. The ‘omnivorous badger dilemma’: Towards an integration of nutrition with the dietary niche in wild mammals. Mamm. Rev. 49(4), 324–339. https://doi.org/10.1111/mam.12164 (2019).Article 

Google Scholar 
71.Noonan, M. J. et al. Climate and the individual: Inter-annual variation in the autumnal activity of the European badger (Meles meles). PLoS ONE https://doi.org/10.1371/journal.pone.0083156 (2014).Article 
PubMed 
PubMed Central 

Google Scholar 
72.Tsunoda, M., Newman, C., Buesching, C. D., Macdonald, D. W. & Kaneko, Y. Badger setts provide thermal refugia, buffering changeable surface weather conditions. J. Therm. Biol. 74, 226–233. https://doi.org/10.1016/j.jtherbio.2018.04.005 (2018).Article 
PubMed 

Google Scholar  More

Forest fires are disasters that cause extensive damage to the entire world in economic, ecological, and environmental ways. These fires can be caused by natural reasons, such as high temperatures that can create spontaneous combustion of dry fuel such as sawdust, leaves, lightning, etc., or by human activities, such as unextinguished campfires, arson, inappropriately burned debris, etc1. According to research, 90% of the world’s forest fire incidents have occurred as a result of the abovementioned human carelessness1. The increase in carbon dioxide levels in the atmosphere due to forest fires contributes to the greenhouse effect and climate change. Additionally, ash destroys much of the nutrients in the soil and can cause erosion, which may result in floods and landslides.
At earlier times, forest fires were detected using watchtowers, which were not efficient because they were based on human observations. In recent history and even the present day, several forest fire detection methods have been implemented, such as watchtowers, satellite image processing methods, optical sensors, and digital camera-based methods2, although there are many drawbacks, such as inefficiency, power consumption, latency, accuracy and implementation costs. To address these drawbacks, a forest fire detection system using wireless sensor networks is proposed in this paper.Wireless sensor networks (WSNs) are self-configured and infrastructure-free wireless networks that help monitor physical or environmental conditions and pass these data through the network to a designated location or sink where the data can be observed and analyzed3. Efficiency and low power consumption are the major advantages of a WSN. In the proposed detection system, wireless sensor nodes are deployed according to cellular architecture to cover the entire area with sensors to monitor temperature, relative humidity, light intensity level, and carbon monoxide (CO) level using a microcontroller, transceiver module, and power components. The power supply to the sensor node is provided using batteries as the primary power supply, and solar panels are used as the secondary power supply. These sensor nodes are specially designed with a spherical shape to withstand damage caused by environmental conditions as well as animals.The sensor readings for each parameter are checked with a preset threshold ratio and a ratio that is calculated continuously in the node in real time, and only the ratios that exceed the preset ratio are sent from the sensor node to the base station for further analytical processing. The network utilized for this transmission is in the architecture of tree topology considering facts such as low power consumption, reduced latency, less complexity, etc. Cluster heads are used in this network to gather data from several sensor nodes and pass them on to the base station or the gateway node. The gateway node is an interface that connects the network with the secondary analysis process.For the analysis process, a machine learning regression model was used along with threshold ratio analysis to enhance detection accuracy. For the training and testing process of the model, data were collected during the fire and no fire situations in different areas and under different climatic conditions. During the data collection process, 7000 data samples were collected, where a data sample included temperature, relative humidity, light intensity level, and CO level at a particular time. Eighty percent of the collected data were randomly used as training data for the model, and the remaining 20% were used as test data.If the outcome of the machine learning model indicates a fire in a specific area, a text message will be sent to the mobile phone numbers of the authorized officers in responsible units. As this process is designed with a minimum delay, the fire can be detected within the initial stage, and the responsible parties can take necessary actions in a shorter period, which will minimize the damage.Related workForest fire detection has been a focus of many researchers for the last decade because of increased forest fire case reports from all over the world due to severe damage to society and the environment. Many methods have been proposed to detect forest fires, such as camera-based systems, WSN-based systems, and machine learning application-based systems, with both positive and negative aspects and performance figures of detection. Due to the higher probability of accurate and early detection due to the use of multiple sensor sources and deployment of sensor nodes in areas not visible to satellites, wireless sensor networks have a more positive outlook, and they have become the more applicable technology in many fields4.Many researchers have focused on environmental parameters, such as air temperature, relative humidity, barometric pressure, sound, light intensity, soil moisture, and wind speed and direction, along with gases, such as CO, CO2, methane, H2, and hydrocarbons apart from smoke, to detect forest fire conditions by considering the variations in these parameters during a fire5,6, and sensors have been selected according to the range, sensitivity, power consumption, and cost7,8.As supplying power to a sensor node is a challenging task in forested areas, utilizing only battery options is difficult because they do not last long, and distributing power using a wire would require a higher cost to deploy in a large forest. Therefore, many researchers have proposed solar-powered systems as secondary power sources along with rechargeable batteries as the main power source4,6, while some researchers have proposed solar batteries because they last longer9. To reduce the power consumption of sensor nodes, techniques such as keeping selected components active while others are deactivated have been proposed10,11,12.Most WSN-based detection systems are centered around a base station due to the memory and processing limitations of the nodes. Important and partially processed data are transferred to the base station through wireless media for processing and enabling relevant actions, while the base station also acts as the gateway between the sensor nodes and the system user4,9.When constructing a WSN, communicating data among the relevant entities is the main objective, and star topology and mesh topology-based networks have been proposed in many papers because of the different attributes in their systems. A mesh topology was chosen over a star topology because of its ability to self-organize, self-configure, and automatically establish among nodes in a network13. As a smaller number of nodes involved for transmission results in minimum energy consumption, concepts based on cluster heads have been used14. To minimize the loss of energy and data packages during transfer, a cluster-tree network topology structure was proposed15. Considering the sensing range of a node, fault tolerance, and energy consumption, a paper has proposed applying the on-demand k-coverage technique that provides event detection using static nodes with variable sensing ranges. This technique utilizes the maximum detection performance with the minimum power consumption for an event16. A survey on rare event detection has mentioned many event detection strategies that deliver maximized detection capability, minimized detection delay and low energy consumption, such as duty cycle, component deactivation, overpopulation/node redundancy, collaboration, and energy harvesting17.To reduce deployment cost and power consumption, a paper proposes a novel localization scheme that divides the whole forest area into different grids and allocates them to respective zones with another 8 neighboring grids. One centroid node from those grids, which is called the initiator node, predicts whether the zone is highly active (HA), medium active (MA), and low active (LA). Here, HA zones send data continuously to the base node through the interior node, MA zones send data periodically, and LA zones do not transfer data in the status that manages power consumption effectively18. Another author proposed obtaining data from the sensors every 2 min if there is the potential of a forest fire or obtaining data every 15 min otherwise to reduce the energy wastage19.To place sensor nodes in the most effective configuration to detect fire conditions, a sensor node was proposed at three different heights to perform different parameter measurements, while some authors have suggested covering sensor nodes to avoid direct sunlight exposure and minimize the false alarm rate4,10. As the network connectivity of service providers in forest areas is not robust, communication techniques that use dedicated network paths such as LoRa, ZigBee, and XBee have been used as the communication infrastructure. When considering attributes such as transmission range, high security, low power consumption (LPWAN protocol), and other relevant configurations, most papers have suggested using the LoRa module for transmission13,19.Most papers have suggested having threshold value-based fire detection on a sensor node, and if the exceeded threshold has remained the same, then a sink determines the location and will send an alarm to the fire department11. Because of the environmental parameter variations according to the place and time, threshold values are configured by the user considering geographic situations, climatic changes, seasonal changes, etc. after sensors obtain the data from the surrounding10.A fusion information process was proposed, where information from multiple sources is considered in making the final decision, which is better than using those sources individually, and two algorithms based on the threshold ratio method and Dempster-Shafer theory were used4. To enhance the detection accuracy, machine learning applications have been proposed in many papers based on different machine learning approaches, such as support vector machine (SVM) classification14 and regression techniques, such as logistic regression. However, applying machine learning techniques to fire detection systems has many limitations, such as the limited amount of energy, the energy required for data processing, the short range of communication and limited computations, the complexity of ML algorithms when executing on sensor nodes, and the difficulty of being distributed on every sensor node20,21. More