Ecology

1.Zehr, J. P. & Capone, D. G. Changing perspectives in marine nitrogen fixation. Science 9514, 729 (2020).
Google Scholar 
2.Karl, D. et al. Dinitrogen fixation in the world’s oceans. Biogeochemistry 57–58, 47–98 (2002).
Google Scholar 
3.Dugdale, R. & Wilkerson, F. in Primary Productivity and Biogeochemical Cycles in the Sea (eds Falkowski, P. G. et al.) 107–122 (Springer, 1992).4.Carpenter, E. J. & Capone, D. G. in Nitrogen in the Marine Environment 2nd edn (eds Capone, D. G., Bronk, D. A., Mulholland, M. R. & Carpenter, E. J.) Ch. 4 (Elsevier, 2008).5.Gruber, N. & Sarmiento, J. L. Global patterns of marine nitrogen fixation and denitrification. Global Biogeochem. Cycles 11, 23–266 (1997).
Google Scholar 
6.Buchanan, P. J., Chase, Z., Matear, R. J., Phipps, S. J. & Bindoff, N. L. Marine nitrogen fixers mediate a low latitude pathway for atmospheric CO2 drawdown. Nat. Commun. https://doi.org/10.1038/s41467-019-12549-z (2019).7.Monteiro, F. M., Follows, M. J. & Dutkiewicz, S. Distribution of diverse nitrogen fixers in the global ocean. Global Biogeochem. Cycles 24, 1–16 (2010).
Google Scholar 
8.Church, M. J., Björkman, K. M., Karl, D. M., Saito, M. A. & Zehr, J. P. Regional distributions of nitrogen-fixing bacteria in the Pacific Ocean. Limnol. Oceanogr. 53, 63–77 (2008).CAS 

Google Scholar 
9.Monteiro, F. M., Dutkiewicz, S. & Follows, M. J. Biogeographical controls on the marine nitrogen fixers. Global Biogeochem. Cycles 25, 1–8 (2011).
Google Scholar 
10.Dutkiewicz, S., Ward, B. A., Monteiro, F. & Follows, M. J. Interconnection of nitrogen fixers and iron in the Pacific Ocean: theory and numerical simulations. Global Biogeochem. Cycles 26, 1–16 (2012).
Google Scholar 
11.Walworth, N. G. et al. Nutrient-colimited Trichodesmium as a nitrogen source or sink in a future ocean. Appl. Environ. Microbiol. 84, 1–14 (2018).CAS 

Google Scholar 
12.McGillicuddy, D. J. Jr. Do Trichodesmium spp. populations in the North Atlantic export most of the nitrogen they fix? Global Biogeochem. Cycles 28, 103–114 (2014).CAS 

Google Scholar 
13.Carpenter, E. J. & Romans, K. Major role of the cyanobacterium Trichodesmium in nutrient cycling in the North Atlantic Ocean. Science 254, 1989–1992 (1991).
Google Scholar 
14.Bergman, B., Sandh, G., Lin, S., Larsson, J. & Carpenter, E. J. Trichodesmium – a widespread marine cyanobacterium with unusual nitrogen fixation properties. FEMS Microbiol. Rev. 37, 286–302 (2013).CAS 
PubMed 

Google Scholar 
15.Capone, D. G. Trichodesmium, a globally significant marine cyanobacterium. Science 276, 1221–1229 (1997).CAS 

Google Scholar 
16.Gallon, J. R. The oxygen sensitivity of nitrogenase: a problem for biochemists and micro-organisms. Trends Biochem. Sci. 6, 19–23 (1981).CAS 

Google Scholar 
17.Saito, M. A. et al. Iron conservation by reduction of metalloenzyme inventories in the marine diazotroph Crocosphaera watsonii. Proc. Natl Acad. Sci. USA 108, 2184–2189 (2011).CAS 
PubMed 
PubMed Central 

Google Scholar 
18.Dron, A. et al. Light-dark (12:12) cycle of carbon and nitrogen metabolism in Crocosphaera watsonii WH8501: relation to the cell cycle. Environ. Microbiol. 14, 967–981 (2012).CAS 
PubMed 

Google Scholar 
19.Mohr, W., Intermaggio, M. P. & LaRoche, J. Diel rhythm of nitrogen and carbon metabolism in the unicellular, diazotrophic cyanobacterium Crocosphaera watsonii WH8501. Environ. Microbiol. 12, 412–421 (2010).CAS 
PubMed 

Google Scholar 
20.Flores, E. & Herrero, A. Compartmentalized function through cell differentiation in filamentous cyanobacteria. Nat. Rev. Microbiol. 8, 39–50 (2010).CAS 
PubMed 

Google Scholar 
21.Burnat, M., Herrero, A. & Flores, E. Compartmentalized cyanophycin metabolism in the diazotrophic filaments of a heterocyst-forming cyanobacterium. Proc. Natl Acad. Sci. USA 111, 3823–3828 (2014).CAS 
PubMed 
PubMed Central 

Google Scholar 
22.Sherman, D. M., Tucker, D. & Sherman, L. A. Heterocyst development and localization of cyanophycin in N2-fixing cultures of Anabaena sp. PCC 7120 (Cyanobacteria). J. Phycol. 941, 932–941 (2000).
Google Scholar 
23.Lamont, H. C., Silvester, W. B. & Torrey, J. G. Nile red fluorescence demonstrates lipid in the envelope of vesicles from N2-fixing cultures of Frankia. Can. J. Microbiol. 34, 656–660 (1988).CAS 

Google Scholar 
24.Saino, T. Diel variation in nitrogen fixation by a marine blue-green alga, Trichodesmium thiebautii. Deep Sea Res. 25, 1259–1263 (1978).
Google Scholar 
25.Saino, T. & Hattori, A. Aerobic nitrogen fixation by the marine non-heterocystous cyanobacterium Trichodesmium (Oscillatoria) spp.: its protective mechanism against oxygen. Mar. Biol. 70, 251–254 (1982).
Google Scholar 
26.Berman-Frank, I. et al. Segregation of nitrogen fixation and oxygenic photosynthesis in the marine cyanobacterium Trichodesmium. Science 294, 1534–1537 (2001).CAS 
PubMed 

Google Scholar 
27.Ohki, K. & Taniuchi, Y. Detection of nitrogenase in individual cells of a natural population of Trichodesmium using immunocytochemical methods for fluorescent cells. J. Oceanogr. 65, 427–432 (2009).CAS 

Google Scholar 
28.Eichner, M. et al. N2 fixation in free-floating filaments of Trichodesmium is higher than in transiently suboxic colony microenvironments. New Phytol. 222, 852–863 (2019).CAS 
PubMed 

Google Scholar 
29.Ohki, K. Intercellular localization of nitrogenase in a non-heterocystous cyanobacterium (cyanophyte), Trichodesmium sp. NIBB1067. J. Oceanogr. 64, 211–216 (2008).CAS 

Google Scholar 
30.Ohki, K., Zehr, F. & Fujita, Y. Regulation of nitrogenase activity in relation to the light-dark regime in the filamentous non-heterocystous cyanobacterium Trichodesmium sp. NIBB 1067. J. Gen. Microbiol. 138, 2679–2685 (1992).CAS 

Google Scholar 
31.Finzi-Hart, J. A. et al. Fixation and fate of C and N in the cyanobacterium Trichodesmium using nanometer-scale secondary ion mass spectrometry. Proc. Natl Acad. Sci. USA 106, 9931 (2009).CAS 

Google Scholar 
32.Sandh, G., El-Shehawy, R., Díez, B. & Bergman, B. Temporal separation of cell division and diazotrophy in the marine diazotrophic cyanobacterium Trichodesmium erythraeum IMS101. FEMS Microbiol. Lett. 295, 281–288 (2009).CAS 
PubMed 

Google Scholar 
33.Küpper, H. et al. Traffic lights in Trichodesmium. Regulation of photosynthesis for nitrogen fixation studied by chlorophyll fluorescence kinetic microscopy. Plant Physiol. 135, 2120–2133 (2019).
Google Scholar 
34.Ohki, K. & Fujita, Y. Aerobic nitrogenase activity measured as acetylene reduction in the marine non-heterocystous cyanobacterium Trichodesmium spp. grown under artificial conditions. Mar. Biol. 98, 111–114 (1988).CAS 

Google Scholar 
35.Waterbury, J. B. & Willey, J. M. Isolation and growth of marine planktonic Cyanobacteria. Methods Enzymol. 167, 100–105 (1988).CAS 

Google Scholar 
36.Chen, Y. B., Zehr, J. P. & Mellon, M. Growth and nitrogen fixation of the diazotrophic filamentous nonheterocystous cyanobacterium Trichodesmium sp. IMS 101 in defined media: evidence for a circadian rhythm. J. Phycol. 32, 916–923 (1996).
Google Scholar 
37.Berman-Frank, I., Bidle, K. D., Haramaty, L. & Falkowski, P. G. The demise of the marine cyanobacterium, Trichodesmium spp., via an autocatalyzed cell death pathway. Limnol. Oceanogr. 49, 997–1005 (2004).
Google Scholar 
38.Bell, P. R. F. et al. Laboratory culture studies of Trichodesmium isolated from the Great Barrier Reef lagoon, Australia. Hydrobiologia 532, 9–21 (2005).
Google Scholar 
39.Tzubari, Y., Magnezi, L., Be’Er, A. & Berman-Frank, I. Iron and phosphorus deprivation induce sociality in the marine bloom-forming cyanobacterium Trichodesmium. ISME J. 12, 1682–1693 (2018).CAS 
PubMed 
PubMed Central 

Google Scholar 
40.Held, N. A., McIlvin, M. R., Moran, D. M., Laub, M. T. & Saito, M. A. Unique patterns and biogeochemical relevance of two-component sensing in marine bacteria. mSystems 4, 1–16 (2019).
Google Scholar 
41.Aryal, U. K. & Sherman, L. A. in Cyanobacteria Omics Manipulation (ed. Los, D. A.) Ch. 6 (Caister Academic Press, 2017).42.Held, N. A. et al. Co-occurrence of Fe and P stress in natural populations of the marine diazotroph Trichodesmium. Biogeosciences 17, 2537–2551 (2020).
Google Scholar 
43.Klugkist, J., Haaker, H., Wassink, H. & Veeger, C. The catalytic activity of nitrogenase in intact Azotobacter vinelandii cells. Eur. J. Biochem. 146, 509–515 (1985).CAS 
PubMed 

Google Scholar 
44.Zehr, J. P., Wyman, M., Miller, V., Capone, D. G. & Duguay, L. Modification of the Fe protein of nitrogenase in natural populations of Trichodesmium thiebautii. Appl. Environ. Microbiol. 59, 669–676 (1993).CAS 
PubMed 
PubMed Central 

Google Scholar 
45.Rodriguez, I. B. & Ho, T.-Y. Diel nitrogen fixation pattern of Trichodesmium: the interactive control of light and Ni. Sci. Rep. 4, 4445 (2014).PubMed 
PubMed Central 

Google Scholar 
46.Eichner, M., Kranz, S. A. & Rost, B. Combined effects of different CO2 levels and N sources on the diazotrophic cyanobacterium Trichodesmium. Physiol. Plant. 152, 316–330 (2014).CAS 
PubMed 
PubMed Central 

Google Scholar 
47.Hutchins, D. A. et al. Irreversibly increased nitrogen fixation in Trichodesmium experimentally adapted to elevated carbon dioxide. Nat. Commun. 6, 1–7 (2015).
Google Scholar 
48.Levitan, O. et al. Combined effects of CO2 and light on the N2-fixing cyanobacterium Trichodesmium IMS101: a mechanistic view. Plant Physiol. 154, 346–356 (2010).CAS 
PubMed 
PubMed Central 

Google Scholar 
49.Villareal, T. A. & Carpenter, E. J. Buoyancy regulation and the potential for vertical migration in the oceanic cyanobacterium Trichodesmium. Microb. Ecol. 45, 1–10 (2003).CAS 
PubMed 

Google Scholar 
50.Rabouille, S., Staal, M., Stal, L. J. & Soetaert, K. Modeling the dynamic regulation of nitrogen fixation in the cyanobacterium Trichodesmium sp. Appl. Environ. Microbiol. 72, 3217–3227 (2006).CAS 
PubMed 
PubMed Central 

Google Scholar 
51.Breitbarth, E., Wohlers, J., Kläs, J., LaRoche, J. & Peeken, I. Nitrogen fixation and growth rates of Trichodesmium IMS-101 as a function of light intensity. Mar. Ecol. Prog. Ser. 359, 25–36 (2008).CAS 

Google Scholar 
52.Chen, Y. B. et al. Circadian rhythm of nitrogenase gene expression in the diazotrophic filamentous nonheterocystous cyanobacterium Trichodesmium sp. strain IMS101. J. Bacteriol. 180, 3598–3605 (1998).CAS 
PubMed 
PubMed Central 

Google Scholar 
53.Rabouille, S., Staal, M., Stal, L. J. & Soetaert, K. Modeling the dynamic regulation of nitrogen fixation in the Cyanobacterium Trichodesmium sp. Appl. Environ. Microbiol. 72, 3217–3227 (2006).CAS 
PubMed 
PubMed Central 

Google Scholar 
54.Capone, D. G., O’Neill, J. M., Zehr, J. & Carpenter, E. J. Basis for diel variation in nitrogenase activity in the marine planktonic cyanobacterium Trichodesmium thiebautti. Appl. Environ. Microbiol. 56, 3532–3536 (1990).CAS 
PubMed 
PubMed Central 

Google Scholar 
55.Gründel, M., Scheunemann, R., Lockau, W. & Zilliges, Y. Impaired glycogen synthesis causes metabolic overflow reactions and affects stress responses in the cyanobacterium Synechocystis sp. PCC 6803. Microbiology 158, 3032–3043 (2012).PubMed 

Google Scholar 
56.Jackson, S. A., Eaton-Rye, J. J., Bryant, D. A., Posewitz, M. C. & Davies, F. K. Dynamics of photosynthesis in a glycogen-deficient glgC mutant of Synechococcus sp. strain PCC 7002. Appl. Environ. Microbiol. 81, 6210–6222 (2015).CAS 
PubMed 
PubMed Central 

Google Scholar 
57.Boatman, T. G., Davey, P. A., Lawson, T. & Geider, R. J. The physiological cost of diazotrophy for Trichodesmium erythraeum IMS101. PLoS ONE 13, 1–24 (2018).
Google Scholar 
58.Chappell, P. D., Moffett, J. W., Hynes, A. M. & Webb, E. A. Molecular evidence of iron limitation and availability in the global diazotroph Trichodesmium. ISME J. 6, 1728–1739 (2012).CAS 
PubMed 
PubMed Central 

Google Scholar 
59.Chappell, P. D. & Webb, E. A. A molecular assessment of the iron stress response in the two phylogenetic clades of Trichodesmium. Environ. Microbiol. 12, 13–27 (2010).CAS 
PubMed 

Google Scholar 
60.Walsby, A. E. The properties and buoyancy-providing role of gas vacuoles in Trichodesmium Ehrenberg. Br. Phycol. J. 13, 103–116 (1978).
Google Scholar 
61.Villareal, T. A. & Carpenter, E. J. Diel buoyancy regulation in the marine diazotrophic cyanobacterium Trichodesmium thiebautii. Limnol. Oceanogr. 35, 1832–1837 (1990).
Google Scholar 
62.Romans, K. M., Carpenter, E. J. & Bergman, B. Buoyancy regulation in the colonial diazotrophic cyanobacterium Trichodesmium tenue: ultrastructure and storage of carbohydrate, polyphosphate, and nitrogen. J. Phycol. 30, 935–942 (1994).
Google Scholar 
63.Wang, L. et al. Molecular structure of glycogen in Escherichia coli. Biomacromolecules 20, 2821–2829 (2019).CAS 
PubMed 

Google Scholar 
64.Berman-Frank, I., Cullen, J. T., Shaked, Y., Sherrell, R. M. & Falkowski, P. G. Iron availability, cellular iron quotas, and nitrogen fixation in Trichodesmium. Limnol. Oceanogr. 46, 1249–1260 (2001).CAS 

Google Scholar 
65.Kustka, A. B. et al. Iron requirements for dinitrogen- and ammonium-supported growth in cultures of Trichodesmium (IMS 101): comparison with nitrogen fixation rates and iron:carbon ratios of field populations. Limnol. Oceanogr. 49, 1224 (2004).CAS 

Google Scholar 
66.Paerl, H. W., Prufert-Bebout, I. L. E., Guo, C. & Carolina, N. Iron-stimulated N2 fixation and growth in natural and cultured populations of the planktonic marine cyanobacteria Trichodesmium spp. Appl. Environ. Microbiol. 60, 1044–1047 (1994).CAS 
PubMed 
PubMed Central 

Google Scholar 
67.Rubin, M., Berman-Frank, I. & Shaked, Y. Dust- and mineral-iron utilization by the marine dinitrogen-fixer Trichodesmium. Nat. Geosci. 4, 529–534 (2011).CAS 

Google Scholar 
68.Polyviou, D. et al. Desert dust as a source of iron to the globally important diazotroph Trichodesmium. Front. Microbiol. 8, 1–12 (2018).
Google Scholar 
69.Basu, S. & Shaked, Y. Mineral iron utilization by natural and cultured Trichodesmium and associated bacteria. Limnol. Oceanogr. 63, 2307–2320 (2018).CAS 

Google Scholar 
70.Held, N. A. et al. Mechanisms and heterogeneity of in situ mineral processing by the marine nitrogen fixer Trichodesmium revealed by single-colony metaproteomics. ISME Commun. 1, 35 (2021).
Google Scholar 
71.Basu, S., Gledhill, M., de Beer, D., Prabhu Matondkar, S. G. & Shaked, Y. Colonies of marine cyanobacteria Trichodesmium interact with associated bacteria to acquire iron from dust. Commun. Biol. 2, 1–8 (2019).CAS 

Google Scholar 
72.Tyrrell, T. et al. Large-scale latitudinal distribution of Trichodesmium spp. in the Atlantic Ocean. J. Plankton Res. 25, 405–416 (2003).CAS 

Google Scholar 
73.Robson, R. L. & Postgate, J. R. Oxygen and hydrogen in biological nitrogen fixation. Annu. Rev. Microbiol. 34, 183–207 (1980).74.Zehr, J. P. Nitrogen fixation by marine cyanobacteria. Trends Microbiol. 19, 162–173 (2011).CAS 
PubMed 

Google Scholar 
75.Bergman, B. & Carpenter, E. J. Nitrogenase confined to randomly distributed trichomes in the marine cyanobacterium Trichodesmium thiebautii. J. Phycol. 27, 158–165 (1991).CAS 

Google Scholar 
76.Inomura, K., Wilson, S. T. & Deutsch, C. Mechanistic model for the coexistence of nitrogen fixation and photosynthesis in marine Trichodesmium. mSystems 4, 1–13 (2019).
Google Scholar 
77.Janson, S., Matveyev, A. & Bergman, B. The presence and expression of hetR in the non-heterocystous cyanobacterium Symploca PCC 8002. FEMS Microbiol. Lett. 168, 173–179 (1998).CAS 
PubMed 

Google Scholar 
78.Zhang, J. Y., Chen, W. L. & Zhang, C. C. hetR and patS, two genes necessary for heterocyst pattern formation, are widespread in filamentous nonheterocyst-forming cyanobacteria. Microbiology 155, 1418–1426 (2009).CAS 
PubMed 

Google Scholar 
79.Moore, J. K., Doney, S. C., Glover, D. M. & Fung, I. Y. Iron cycling and nutrient-limitation patterns in surface waters of the world ocean. Deep Sea Res. 2 Top. Stud. Oceanogr. 49, 463–507 (2001).
Google Scholar 
80.Chisholm, S. W. in Primary Productivity and Biogeochemical Cycles in the Sea (eds Falkowski, P. G. et al.) 213–237 (Springer, 1992).https://doi.org/10.1007/978-1-4899-0762-2_1281.Young, K. D. The selective value of bacterial shape. Microbiol. Mol. Biol. Rev. 70, 660–703 (2006).PubMed 
PubMed Central 

Google Scholar 
82.Lu, X. & Zhu, H. Tube-gel digestion: a novel proteomic approach for high-throughput analysis of membrane proteins. Mol. Cell Proteom. 4, 1948–1958 (2005).CAS 

Google Scholar 
83.Saito, M. A. et al. Multiple nutrient stresses at intersecting Pacific Ocean biomes detected by protein biomarkers. Science 345, 1173–1177 (2014).CAS 
PubMed 

Google Scholar 
84.McIlvin, M. R. & Saito, M. A. Online nanoflow two-dimension comprehensive active modulation reversed phase-reversed phase liquid chromatography high-resolution mass spectrometry for metaproteomics of environmental and microbiome samples. J. Proteome Res. 20, 4589–4597 (2021).CAS 
PubMed 

Google Scholar 
85.Lee, M. D. et al. Transcriptional activities of the microbial consortium living with the marine nitrogen-fixing cyanobacterium Trichodesmium reveal potential roles in community-level nitrogen cycling. Appl. Environ. Microbiol. 84, AEM.02026-17 (2017).
Google Scholar 
86.Zhang, Y., Wen, Z., Washburn, M. P. & Florens, L. Refinements to label-free proteome quantitation: how to deal with peptides shared by multiple proteins. Anal. Chem. 82, 2272–2281 (2010).CAS 
PubMed 

Google Scholar 
87.Gallien, S., Bourmaud, A., Kim, S. Y. & Domon, B. Technical considerations for large-scale parallel reaction monitoring analysis. J. Proteom. 100, 147–159 (2014).CAS 

Google Scholar 
88.Pino, L. K. et al. The skyline ecosystem: informatics for quantitative mass spectrometry proteomics. Mass Spectrom. Rev. 176, 139–148 (2019).
Google Scholar 
89.Held, N. A. et al. Mechanisms and heterogeneity of in situ mineral processing by the marine nitrogen fixer Trichodesmium revealed by single-colony metaproteomics. ISME Commun. https://doi.org/10.1038/s43705-021-00034-y (2021).90.White, A. E., Spitz, Y. H. & Letelier, R. M. Modeling carbohydrate ballasting by Trichodesmium spp. Mar. Ecol. Prog. Ser. 323, 35–45 (2006).
Google Scholar 
91.Morrison, F. A. An Introduction to Fluid Mechanics (Cambridge Univ. Press, 2013).92.Hunter, J. D. Matplotlib: a 2D graphics environment. Comput. Sci. Eng. 9, 90–95 (2007).
Google Scholar 
93.Virtanen, P. et al. SciPy 1.0: fundamental algorithms for scientific computing in Python. Nat. Methods 17, 261–272 (2020).CAS 
PubMed 
PubMed Central 

Google Scholar 
94.Hagberg, A. A., Schult, D. A. & Swart, P. J. Exploring network structure, dynamics, and function using NetworkX. In 7th Python Scientific Conference (SciPy 2008) 11–15 (2008). More

Study areaAll simulations were run at a 100 × 100 m resolution for the entire European Alps, which cover ~200,000 km². Elevations reach 4,810 m above sea level at the highest peak (Mont Blanc, elevational data were obtained from ref. 44). Mean annual temperature ranges from approximately −13 up to 16 °C and annual precipitation sums reach up to ~3,600 mm (climatic conditions were obtained from WorldClim45).Species dataTrue presences/absences were derived from complete species lists of 14,040 localized plots covering subalpine and alpine non-forest vegetation of the Alps, compiled from published46 and unpublished data sources. For more information see the supplementary information in ref. 21.Data on demographic rates as well as dispersal parameters were taken from ref. 21, Supplementary Table 1. Detailed values are provided in Supplementary Table 1.Environmental variablesCurrent climate dataMaps of current climatic conditions were generated on the basis of mean, minimum and maximum monthly temperature obtained from WorldClim45 and monthly precipitation sums derived from ref. 47 at a resolutions of 0.5 arcmin and 5 km, respectively. Temperature and precipitation data were downscaled to 100 m as described in ref. 21 and using ordinary kriging with elevation as covariable. As the reference periods of these two datasets do not match (temperature values represent averages for 1950–2000, while precipitation covers 1970–2005) temperature values were adapted using the E-OBS climate grids available online (www.ecad.eu/download/ensembles/download.php). We used these spatially refined temperature and precipitation grids to derive maps of mean annual temperature and mean annual precipitation sum. We chose only two climatic variables to keep models simple and, therefore, interpretation of results more straightforward. In fact, the climatic drivers of population performance and distribution can be more complex48 and vary among species, life history stages and vital rates49. However, as correlations between different descriptors of temperature (such as coldest month or warmest month temperature, Pearson correlation of 0.84) as well as between precipitation variables are high in the European Alps, we chose mean annual temperature and mean annual precipitation sum as they give the most basic description of how climatic conditions change over geographical and elevational gradients.Future climate dataMonthly time series of mean temperature as well as precipitation sums predicted for the twenty-first century were extracted from the Cordex data portal (http://cordex.org). We chose two IPCC5 scenarios from the RCP families representing mild (RCP 2.6) and severe (RCP 8.5) climate change to consider the uncertainty in the future climate predictions. Both scenarios were generated by Météo-France/Centre National de Recherches Météorologiques using the CNRM-ALADIN53 model, fed by output from the global circulation model CNRM-CM5 (ref. 50). The RCP 2.6 scenario assumes that radiative forcing reaches nearly 3 W m−2 (equal to 490 ppm CO2 equivalent) mid-century and will decrease to 2.6 W m−2 by 2100. In the RCP 8.5 scenario, radiative forcing continues to rise throughout the twenty-first century and reaches >8.5 W m−2 (equal to 1,370 ppm CO2 equivalent) in 210024.These time series were statistically downscaled (delta method) by (1) calculating differences (deltas) between monthly temperature and precipitation values hindcasted for the current climatic conditions (mean 1970–2005) and forecasted future values at the original spatial resolution of 11′; (2) spatially interpolating these differences to a resolution of 100 × 100 m using cubic splines and (3) adding them to the downscaled current climate data separately for each climatic variable21,36. Subsequently, we calculated running means (averaged over 35 years) for every tenth year (2030, 2040 through to 2080) for each climatic variable and finally derived, on the basis of the monthly data, mean annual temperature and mean annual precipitation sums for these decadal time steps. The application of 35-yr running means ensures that we use the same time interval for parameterization and prediction. Furthermore, for long-lived species such as alpine plants, running means over long time intervals appear most appropriate to characterize climatic niches33.Soil dataIn addition to the climatic data we used a map of the percentage of calcareous substrate within a cell (5′ longitude × 3′ latitude dissolved to 100 m resolution; further referred to as soil) as described in the supplementary information of ref. 21.Environmental suitability modellingWe parameterized logistic regression models (LRMs) with a logit link function using species presence/absence data as response and the three environmental (two bioclimatic and one soil) variables as predictors. All predictor variables entered the model as second-order polynomials in agreement with the standard unimodal niche concept.From the models, we also derived a threshold value to use for translating predicted probability of occurrence (as a measure of site suitability) into predicted presence or absence of each species at a site (called occurrence threshold, OT, henceforth). The threshold was defined such that it optimized the true skills statistic (TSS), a measure of predictive accuracy derived from comparing observed and predicted presence–absence maps51.Genetic model and niche partitioningSpecies-specific suitability curves for the annual mean temperature gradient derived from the LRMs were partitioned into suitability curves of ecologically different haplotypes of a species as defined by allelic variation in up to three diploid loci (Extended Data Fig. 6). The number of alleles was varied (n = 1, 2, 3, 5 and 10 alleles) as was the proportion of the entire species’ (temperature) niche covered by each haplotype. Models with more than one locus were run with diallelic loci, as otherwise computational efforts would have increased excessively (for each haplotype the number of seeds, juveniles and adults has to be stored and all seeds have to be distributed separately). Each combination of haplotype number and allelic niche size used in a particular simulation is further referred to as setting. Species-specific suitability curves along the other two dimensions (precipitation and soil) remained unpartitioned to ease interpretation of results. The implications of relaxing this assumption by modelling niche partitioning along different environmental gradients are hard to predict. Outcomes would probably depend on the direction and strength of individual specialization along these gradients, whether they are positively or negatively correlated, as well as on how both temperature and precipitation patterns will be affected by climate change. As a consequence, the patterns we found could be re-enforced, for example when the replacement of cold-adapted haplotypes within the species elevational range is further delayed, for example, because those haplotypes adapted to warmer conditions can cope less well with higher precipitation at higher elevations. Vice versa, maladaptation to the warming temperatures might be attenuated, for example, if temperature increase is paralleled by precipitation decrease and emerging drought stress. If, in this case, haplotypes from lower elevations can better cope with both higher temperatures and less water availability than those of median elevations, they may replace the latter faster at these median elevations and hence shorten the phase of maladaptation.Allelic effects were assumed to define the temperature optimum additively. Hence, the heterozygotes’ optimum is always exactly between the optima of the two corresponding homozygotes, corresponding to a codominant genetic model. Finally, all haplotypes corresponding to one setting were assumed to have constant (temperature) niche size, defined as a proportion (ω = 50%, 75% and 100%, for one haplotype only 100%) of the entire species’ (temperature) niche. The temperature niche was computed as the difference between the upper and lower temperature values at which the LRM-derived suitability curve predicted a suitability equal to OT (with precipitation and soil set to the respective optima of the species, also derived from the LRMs). To derive the same geographic distribution under current climatic conditions for each setting, the union of the niches of all haplotypes of one set has to approximate the niche of the single-species model (Extended Data Fig. 6). Note, however, that, the aspired equality of niches is impossible, as the niches resulting from a logistic regression with quadratic terms are always elliptic in shape. Therefore, the optima of the two extreme homozygotes (that is, those carrying haplotypes adapted to the coldest or warmest margins of the entire species’ niche) are fixed at: niche limits ± 1/2 × ω × niche size and all other optima are distributed between them at equal distances (Extended Data Fig. 6 gives a schematic illustration). As a consequence, the performance of the extreme haplotypes, both coldest and warmest, is modelled to be somewhat higher towards the cold and warm margins of the temperature niche, respectively, compared with the performance modelled for the species without intraspecific niche partitioning (compare the overlap of the black with the red and blue curves in Extended Data Fig. 6a). However, as haplotype number did not affect modelled range loss (compare Table 1), this marginal mismatch does not apparently impact our results. Homozygotes were ordered from the cold-adapted A1A1 up to the warm-adapted AnAn.Finally, the suitability curves of the genotypes were assumed to have the same value at their optimum as the species-specific suitability curve at this point (Extended Data Fig. 6).Artificial landscapesArtificial landscapes were defined as a raster of 50 × 112 cells (of 100 × 100 m). These rasters were homogeneous with respect to precipitation and soil carbon conditions which were set to the values optimal for each species according to the LRMs. With respect to temperature, by contrast, we implemented a gradient across the raster running from the minimum –9.1 °C to the maximum +3.8 °C temperature for which the LRM predicts values >OT across all six species. Buffering by 1 °C at both limits was done to avoid truncating simulation results. Further 4 °C at the lower limit were added to consider simulated temperature increase (below). The final temperature range represented by the raster ran from −14.1 to +4.8 °C. Temperature increased linearly along this axis by an increment of 0.171 °C per cell, derived from a 20° slope and a temperature decrease of 0.5 °C per 100 m of elevational change. Along the 50-cell breadth of the landscape, temperature values were kept constant. On the basis of these grids, we implemented a moderate and a severe climate change scenario, characterized by temperature increases of 2 and 4 °C, respectively, over 80 yr. Therefore, temperature of each raster cell increased annually by 0.025 and 0.05 °C, respectively.Simulations of spatiotemporal range dynamicsCATS21 is a spatially and temporarily explicit model operating on a two-dimensional grid (of 100 m mesh size in this case). CATS combines simulations of local species’ demography with species’ distribution models by scaling demographic rates relative to occurrence probabilities (suitabilities) predicted for the respective site or grid cell by the latter. Dispersal among grid cells is modelled as a combination of wind, exozoochoric and endozoochoric (that is, animal dispersal via attachment to the fur or ingestion and defecation, respectively) dispersal of plant seeds. Time proceeds in annual steps. The annually changing occurrence probabilities at each site affect the demography of local populations and hence, eventually, the number of seeds that are produced in each grid cell in the respective year. As a consequence, local populations grow or decrease, become extinct or establish anew and hence the species’ distribution across the whole study area changes as a function of the changing climate.Demographic modellingClimate-dependence of local demography was modelled by linking demographic rates (seed persistence, germination, survival, flowering frequency, seed yield and clonal reproduction) and carrying capacity to occurrence probabilities predicted by LRMs by means of sigmoidal functions. Furthermore, all rates were fixed at OT at a value ensuring stable population sizes; for more information see refs. 21,33. Demographic rates were confined by zero and a species-specific maximum value (Supplementary Table 1), which was assumed to be the same for all genotypes of a species. As a corollary, the demographic rates of all genotypes of a species are the same under optimal environmental conditions but their actual values for a particular site in a particular year differ due to different temperature optima of genotypes. In addition, germination, survival and clonal reproduction were modelled as density-dependent processes to account for intraspecific competition between genotypes. In our application, for all density-dependent functions, the species abundance is defined as the sum of all adult individuals in a given cell, regardless of their genotypes. Density dependence is commonly achieved by multiplying rates with (C – N)/C, where N is the species abundance and C is the (site- and genotype-specific) carrying capacity. This correction for density dependence causes the functions to drop to zero when N approximates C. To avoid the subsequent collapse of population sizes, we defined density-dependent rates as (C – N)/C × rate() + N/C × rate(OT), which ensures stable population sizes at densely populated sites occupied by only one genotype. To account for uncertainty in parameters of demographic rates, we assigned each species two value sets representing the upper and lower end of a plausible range of values on the basis of information derived from databases and own measurements (Supplementary Table 1).The simulations allowed for cross-pollination between genotypes. We used the relative amount of flowers (genotype-specific flowering frequency as defined by the sigmoid curve for the given suitability in the given raster cell for the given year × number of adults of that genotype in the population of that cell) to derive an estimate of the haplotype frequencies in the total pollen produced by the population within a grid cell. For the multiallelic case we allowed for recombination between loci with a recombination rate of 0.1%. These frequencies were set equal to the probability that particular haplotypes are transmitted to each year’s seed yield by pollination. Spatial pollen dispersal was accounted for in the following way: in each cell, 5% of the pollen involved in producing the annual seed yield, was assumed to stem from outside the respective raster cell. The proportions of different haplotypes in this 5% were derived from the overall pollen frequencies in all cells within a 700 m radius around the target cell (following estimates in ref. 52). Subsequently, produced seeds of each genotype were divided into resulting genotypes regarding the adult’s haplotype composition and the haplotype frequencies in the cells’ entire pollen load.Dispersal modellingFor wind dispersal of plant species we parameterized the analytical WALD kernel53 on the basis of measured seed traits and wind speed data from a meteorological station in the Central Alps of Austria. Exozoochorous and endozoochorous plant kernels were parameterized on the basis of correlated random walk simulations for the most frequent mammal dispersal vector in the study area, the chamois (Rupicapra rupicapra L.). For more details, see ref. 33. To account for uncertainties in species-specific dispersal rates, the proportion of seeds dispersed by the more far-reaching zoochorous kernels was assumed either as high (1–5%) or low (0.1–0.5%), setting upper and lower boundaries of a credible range of the dispersal ability of species.Simulation set up and simulation initializationTo test for the effects of climate change on genetic diversity in 2080, we ran CATS over the period 2000 to 2080 for each of the six study species across the entire Alps under a full factorial combination of (1) three niche sizes (50%, 75% and 100%); (2) six numbers of haplotypes (equal to two, three, five and ten alleles for one locus and four and eight for the diallelic two- and three-locus models, respectively); (3) three climate scenarios (current climate, RCP 2.6 and RCP 8.5); and (4) two sets of demographic and dispersal parameters. As a ‘control’ we also ran simulations for all climate scenarios and the two demographic and dispersal parameter sets for a setting with one genotype filling the whole niche of the species. To account for stochastic elements in CATS four replications were run for each combination of ‘treatments’.For simulations in artificial landscapes we used, instead of RCP 2.6 and RCP 8.5, ‘artificial’ climatic scenarios with an assumed warming of 2 and 4 °C, respectively, and no change in precipitation.All simulation runs were started with homozygotic individuals only. As initial distribution, for each simulation run each cell predicted to be environmentally suitable to the species (that is, occurrence probability of species >OT)—and within the real distribution range of the species28 (not relevant for simulations in artificial landscapes, of course)—was assumed to be occupied by an equal number of adults of each (homozygotic) genotype, with the total sum equal to the carrying capacity of the site. To accommodate this arbitrary within-cell genetic mixture of homozygotes (and missing heterozygotes) to actual local conditions we started simulations of range dynamics with a burn-in phase of 200 yr, run under constant current climatic conditions. To have a smooth transition from the burn-in phase under current climate (corresponding to the climate of the years 1970–2005; see current climate data) to future climate projections (starting with 2030) and to derive annual climate series, climate data were linearly interpolated between these two time intervals.Statistical analysisWe evaluated the contribution of climate scenario, haplotype number and haplotype niche size to overall species’ range change as well as the change in the frequency of the warm-adapted haplotype by means of linear models. In these models, log(range size in 2080/range size in 2000) and log(percentage of warm-adapted haplotype in 2080/percentage of warm-adapted haplotype in 2000), averaged over the four replicates and the two demographic and dispersal parameter sets, were the response variables. For the analysis of the change of the warm-adapted haplotype simulation settings with 100% niche size were ignored, as in this setting all haplotypes have the same temperature optimum (that is, neutral genetic variation). Approximate normality of residuals was confirmed by visual inspection.As indicator of the ‘topographic opportunity’ remaining to the species in the real world we calculated the area colonizable at elevations higher than those already occupied at the end of the simulation period. We therefore drew a buffer of 1 km around each cell predicted to be occupied in 2080 and then summed the area within these buffers at a higher elevation than the focal, occupied cell. Overlapping buffer areas were only counted once. This calculation was done for simulations conducted with one full-niche genotype per species only.Sensitivity analysisWe interpret the simulated relative decrease of warm-adapted haplotypes mainly as an effect of (1) the unrestricted expansion of cold-adapted haplotypes at the leading edge and (2) the resistance of the locally predominating haplotype that becomes increasingly maladapted with progressive climate warming, to recruitment by better-adapted haplotypes from below that are either rare or not present at all initially. However, the results achieved, and our conclusions, may be sensitive to assumptions about particular parameter values. Parameters that control the longevity of adult plants, and indirectly the rate of recruitment of new individuals, as well as those controlling gene flow via pollen (instead of seeds) may be particularly influential in this respect. We additionally ran simulations on artificial landscapes under alternative values of these parameters. In particular, we set the maximum age of plants to 10 yr instead of 100 yr and raised the proportion of locally produced pollen assumed to be transported up to 700 m to 10%. Both of these values are thus probably set to rather extreme levels: a maximum age of 10 yr is much shorter than the 30–50 yr assumed to be the standard age of (non-clonal) alpine plants31; and a cross-pollination rate between cells of 10% is high given that among the most important alpine pollinators only bumblebees are assumed to transport pollen >100 m regularly54,55. We add that we ran these additional simulations only in combination with the demographic species parameters set to high values because a short life span combined with low-level demographic parameters did not allow for stable populations of some species, even under current climatic conditions.For individual species, adapting plant age and cross-pollination rate between cells (Extended Data Fig. 7), did change the magnitude of loss of the warm-adapted haplotype. Nevertheless, for all of them the warm-adapted haplotype still became rarer when climate warmed and this effect increased with the level of warming. We are confident that our conclusions are qualitatively insensitive to variation of these parameters within a realistic range.Finally, in the simulations where we assumed a multilocus structure of the temperature niche, the recombination rate may also affect simulation results because it determines the rate by which new haplotypes can emerge. We also tested sensitivity of our simulations to doubling the recombination rate to 0.2%. Again, we found that a higher recombination rate had little qualitative effect on the results. Quantitatively, it resulted in a slightly pronounced relative decrease of the warmth-adapted haplotype in most species (Extended Data Fig. 8).Reporting SummaryFurther information on research design is available in the Nature Research Reporting Summary linked to this article. More

1.Science Plan and Implementation Strategy IGBP Report No. 53/IHDP Report No. 19 (Global Land Project, 2005).2.Millennium Ecosystem Assessment—Ecosystems and Human Well-Being: Desertification Synthesis Encyclopedia of the Anthropocene vols 1–5 (MEA, 2017).3.Huang, J., Yu, H., Guan, X., Wang, G. & Guo, R. Accelerated dryland expansion under climate change. Nat. Clim. Change 6, 166–171 (2015).Article 

Google Scholar 
4.Pimm, S. L. The complexity and stability of ecosystems. Nature 307, 321–326 (1984).Article 

Google Scholar 
5.Tilman, D. & Downing, J. A. Biodiversity and stability in grasslands. Nature 367, 363–365 (1994).Article 

Google Scholar 
6.Belnap, J. Surface disturbances: their role in acceleration desertification. Environ. Monit. Assess. 37, 38–57 (1995).Article 

Google Scholar 
7.Zhao, Y., Jia, R. L. & Wang, J. Towards stopping land degradation in drylands: water-saving techniques for cultivating biocrusts in situ. Land Degrad. Dev. 30, 2336–2346 (2019).Article 

Google Scholar 
8.Rodriguez-Caballero, E. et al. Dryland photoautotrophic soil surface communities endangered by global change. Nat. Geosci. 11, 185–189 (2018).CAS 
Article 

Google Scholar 
9.Coe, K. K. & Sparks, J. P. Physiology-based prognostic modeling of the influence of changes in precipitation on a keystone dryland plant species. Oecologia 176, 933–942 (2014).Article 

Google Scholar 
10.Ferrenberg, S., Tucker, C. L. & Reed, S. C. Biological soil crusts: diminutive communities of potential global importance. Front. Ecol. Environ. 15, 160–167 (2017).Article 

Google Scholar 
11.Belnap, J. & Gillette, D. A. Soil surface disturbance: impacts on potential wind erodibility of sand desert soils in SE Utah, USA. Land Degrad. Dev. 8, 355–362 (1997).Article 

Google Scholar 
12.Rutherford, W. A. et al. Albedo feedbacks to future climate via climate change impacts on dryland biocrusts. Sci. Rep. 7, 44188 (2017).13.Duniway, M. C. et al. Wind erosion and dust from US drylands: a review of causes, consequences, and solutions in a changing world. Ecosphere 10, e02650 (2019).14.Ferrenberg, S., Faist, A. M., Howell, A. & Reed, S. C. Biocrusts enhance soil fertility and Bromus tectorum growth, and interact with warming to influence germination. Plant Soil 429, 77–90 (2018).CAS 
Article 

Google Scholar 
15.Eldridge, D. J. et al. The pervasive and multifaceted influence of biocrusts on water in the world’s drylands. Glob. Change Biol. 26, 6003–6014 (2020).16.Ferrenberg, S., Reed, S. C. & Belnap, J. Climate change and physical disturbance cause similar community shifts in biological soil crusts. Proc. Natl Acad. Sci. USA 112, 12116–12121 (2015).CAS 
Article 

Google Scholar 
17.Reed, S. C. et al. Changes to dryland rainfall result in rapid moss mortality and altered soil fertility. Nat. Clim. Change 2, 752–755 (2012).CAS 
Article 

Google Scholar 
18.Concostrina-Zubiri, L. et al. Biological soil crusts across disturbance-recovery scenarios: effect of grazing regime on community dynamics. Ecol. Appl. 24, 1863–1877 (2014).CAS 
Article 

Google Scholar 
19.Weber, B., Bowker, M., Zhang, Y. & Belnap, J. in Biological Soil Crusts: An Organizing Principle in Drylands (eds Weber, B., Büdel, B. & Belnap, J.) 479–498 (Springer, 2016).20.Reynolds, J. F. et al. Global desertification: building a science for dryland development. Science 316, 847–851 (2007).CAS 
Article 

Google Scholar 
21.Berdugo, M. et al. Global ecosystem thresholds driven by aridity. Science 367, 787–790 (2020).CAS 
Article 

Google Scholar 
22.Bestelmeyer, B. T. et al. Analysis of abrupt transitions in ecological systems. Ecosphere 2, e03360 (2011).Article 

Google Scholar 
23.Bestelmeyer, B. T. et al. Desertification, land use, and the transformation of global drylands. Front. Ecol. Environ. 13, 28–36 (2015).Article 

Google Scholar 
24.Beisner, B., Haydon, D. & Cuddington, K. Alternative stable states in ecology. Front. Ecol. Environ. 1, 376–382 (2003).25.Fukami, T. & Nakajima, M. Community assembly: alternative stable states or alternative transient states? Ecol. Lett. 14, 973–984 (2011).Article 

Google Scholar 
26.Belnap, J. & Büdel, B. in Biological Soil Crusts: An Organizing Principle in Drylands (eds Weber, B., Büdel, B. & Belnap, J.) 305–320 (Springer, 2016).27.Belnap, J. & Warren, S. D. Measuring restoration success: a lesson from Patton’s tank tracks. Ecol. Bull. 79, 33 (1998).28.Belnap, J. & Elderidge, D. in Biological Soil Crusts: Structure, Function and Management (eds Belnap, J. & Lange, O. L.) 363–383 (Springer, 2001).29.Turner, M. G. Disturbance and landscape dynamics in a changing world. Ecology 91, 2833–2849 (2010).Article 

Google Scholar 
30.Scheffer, M. & Carpenter, S. R. Catastrophic regime shifts in ecosystems: linking theory to observation. Trends Ecol. Evol. 18, 648–656 (2003).Article 

Google Scholar 
31.Sala, O. E. & Lauenroth, W. K. Small rainfall events: an ecological role in semiarid regions. Oecologia 53, 301–304 (1982).32.Cayan, D. R. et al. Future dryness in the Southwest US and the hydrology of the early 21st century drought. Proc. Natl Acad. Sci. USA 107, 21271–21276 (2010).CAS 
Article 

Google Scholar 
33.Christensen, N. S., Wood, A. W., Nathalie, V., Lettenmaier, D. P. & Palmer, R. N. The effects of climate change on the hydrology and water resources of the Colorado river basin. Clim. Change 62, 337 (2004).Article 

Google Scholar 
34.Herrick, J. et al. Field soil aggregate stability kit for soil quality and rangeland health evaluations. Catena 44, 27–35 (2001).Article 

Google Scholar 
35.Escolar, C., Martínez, I., Bowker, M. A. & Maestre, F. T. Warming reduces the growth and diversity of biological soil crusts in a semi-arid environment: implications for ecosystem structure and functioning. Phil. Trans. R. Soc. B 367, 3087–3099 (2012).Article 

Google Scholar 
36.Scheffer, M. et al. Creating a safe operating space for iconic ecosystems: manage local stressors to promote resilience to global change. Science 347, 1317–1319 (2015).CAS 
Article 

Google Scholar 
37.Collins, S. L., Micheli, F. & Hartt, L. A method to determine rates and patterns of variability in ecological communities. Oikos 91, 285–293 (2000).Article 

Google Scholar 
38.Rillig, M. C. et al. The role of multiple global change factors in driving soil functions and microbial biodiversity. Science 366, 886–890 (2019).CAS 
Article 

Google Scholar 
39.IPCC. Climate Change 2014: Impacts, Adaptations, and Vulnerability (eds Field, C. B. et al.) (Cambridge Univ. Press, 2014).40.Mirzabaev, A. et al. in IPCC Special Report on Global Warming of 1.5 °C (eds Masson-Delmotte, V. et al.) Ch. 3 (WMO, 2018).41.Torres-Cruz, T. J. et al. Species-specific nitrogenase activity in lichen-dominated biological soil crusts from the Colorado Plateau, USA. Plant Soil 429, 113–125 (2018).CAS 
Article 

Google Scholar 
42.Omernik, J. M. & Griffith, G. E. Ecoregions of the conterminous United States: evolution of a hierarchical spatial framework. Environ. Manag. 54, 1249–1266 (2014).Article 

Google Scholar 
43.Kuske, C. R., Yeager, C. M., Johnson, S., Ticknor, L. O. & Belnap, J. Response and resilience of soil biocrust bacterial communities to chronic physical disturbance in arid shrublands. ISME J. 6, 886–897 (2011).Article 

Google Scholar 
44.Tucker, C. L., Ferrenberg, S. & Reed, S. C. Climatic sensitivity of dryland soil CO2 fluxes differs dramatically with biological soil crust successional state. Ecosystems 22, 15–32 (2018). https://doi.org/10.1007/s10021-018-0250-445.Cable, J. M. & Huxman, T. E. Precipitation pulse size effects on Sonoran Desert soil microbial crusts. Oecologia 141, 317–324 (2004).Article 

Google Scholar 
46.Karl, T. R., Knight, R. W. & Plummer, N. Trends in high-frequency climate variability in the twentieth century. Nature 377, 217–220 (1995).CAS 
Article 

Google Scholar 
47.Kunkel, K. E., Easterling, D. R., Redmond, K. & Hubbard, K. Temporal variations of extreme precipitation events in the United States: 1895–2000. Geophys. Res. Lett. 30, 1895–2000 (2003).Article 

Google Scholar 
48.Kim, J. A projection of the effects of the climate change induced by increased CO2 on extreme hydrologic events in the Western US. Clim. Change 68, 153–168 (2005).CAS 
Article 

Google Scholar 
49.Smith, S. J. et al. Climate change impacts for the conterminous USA: an integrated assessment part 1. Scenarios and context. Clim. Change 69, 7–25 (2005). https://doi.org/10.1007/1-4020-3876-3_250.Schwinning, S., Belnap, J., Bowling, D. R. & Ehleringer, J. R. Sensitivity of the Colorado Plateau to change: climate, ecosystems, and society. Ecol. Soc. 13, 28 (2008).51.Daly, C. et al. Physiographically sensitive mapping of climatological temperature and precipitation across the conterminous United States. Int. J. Climatol. 28, 2031–2064 (2008).Article 

Google Scholar 
52.Jonasson, S. The point intercept method for non-destructive estimation of biomass. Phytocoenologia 11, 385–388 (1983).Article 

Google Scholar 
53.R Core Team. R: A Language and Environment for Statistical Computing (R Foundation for Statistical Computing, 2019).54.Wickham, H. ggplot2: Elegant Graphics for Data Analysis (Springer-Verlag, 2016).55.Oksanen, A. J. et al. Vegan: Community Ecology Package. Rpackage version 2.5-7 https://CRAN.R-project.org/package=vegan (2020).56.Anderson, M. J. A new method for non-parametric multivariate analysis of variance. Austral Ecol. 26, 32–46 (2008).Article 

Google Scholar 
57.Venables, W. & Ripley, B. Modern Applied Statistics with S. (Springer, 2002).58.Lenth, R., Singmann, H., Love, J., Buerkner, P. & Herve, M. Package ‘emmeans’ https://github.com/rvlenth/emmeans (2018).59.Signorell, A. DescTools: Tools for Descriptive Statistics (2021).60.Hallett, L. M. et al. codyn: an R package of community dynamics metrics. Methods Ecol. Evol. 7, 1146–1151 (2016).61.Wood, S. N. Generalized Additive Models: An Introduction with R 2nd edn 1–476 (CRC/Taylor & Francis, 2017).62.Bürkner, P. C. brms: an R package for Bayesian multilevel models using Stan. J. Stat. Softw. 80, 1–28 (2017).63.Gelman, A. & Rubin, D. B. Inference from iterative simulation using multiple sequences. Stat. Sci. 7, 547–511 (1992).
Google Scholar 
64.Vehtari, A., Gelman, A. & Gabry, J. Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC. Stat. Comput. 27, 1413–1432 (2017).Article 

Google Scholar 
65.Modrák, M., Barrett, M., Weber, F. & Coronado, E. bayesplot: Plotting for Bayesian Models. R package version 1.8.0 https://mc-stan.org/bayesplot/ (2021).66.Makowski, D., Ben-Shachar, M. & Lüdecke, D. bayestestR: describing effects and their uncertainty, existence and significance within the Bayesian framework. J. Open Source Softw. 4, 1541 (2019).Article 

Google Scholar 
67.Phillips, M. L., Howell, A., Lauria, C. M., Belnap, J. & Reed, S. C. Data and software code from two long-term experiments (1996–2011 and 2005–2018) at three sites on the Colorado Plateau of North America (US Geological Survey, 2021); https://doi.org/10.5066/P9RUN1TP More

Environmental variable quantification along an altitudinal gradientThis study area included 22 sampling sites, and 66 samples, classified into three transects, namely Transect 1 (1047–1587 m), Transect 2 (876–3070 m), and Transect 3 (1602–2110 m). Significant differences in soil properties and plant parameters were observed along the three studied altitudinal transects (P  0.05%, while the remaining bacteria were merged into an “others” class. As shown in Fig. 1B, the proportion of Actinobacteria, Alphaproteobacteria and Gammaproteobacteria at each elevation was 45%, whereas Deltaproteabacteria, Acidobacteria_Subgroup_6, and Gemmatimonadetes were prevalent at low levels in most soil samples. At the genus level (Supplementary Fig. 1), 76 genera were detected in the research areas, with the dominant genera including norank_f_67-14_ o_Solirubrobacterales (5.72%), Rubrobacter (4.35%), Solirubrobacter (2.83%), Pseudonocardia (2.26%) and Bradyrhizobium (2.19%) and less than 0.01% of the bacterial genera were classified into others.Figure 1Bacterial community composition variations at the phylum (A) and class (B) levels in soil samples collected at different levels. These were done in R (v3.3.1, http://www.R-project.org).Full size imageBacterial community composition varies along elevation gradientsWe next sought to analyze the differences in relative bacterial abundance at the phylum level among Transects 1–3 (Fig. 2). Significant differences in the relative abundance of Actinobacteria, Proteobacteria, Acidobacteria, Verrucomicrobia, Firmicutes, and Rokubacteria were detected in samples from the different transects (Fig. 2A). The relative abundance of Actinobacteria and Firmicutes in Transect 1 (48.64% and 1.89%, respectively) was significantly higher than in Transect 2 (38.43% and 1.49%, respectively) and Transect 3 (39.63% and 0.98%, respectively) (P  More

1.Carter, D. O., Yellowlees, D. & Tibbett, M. Cadaver decomposition in terrestrial ecosystems. Naturwissenschaften 94(1), 12–24 (2007).ADS 
CAS 
PubMed 

Google Scholar 
2.Weathers, K. C., Strayer, D. L. & Likens, G. E. Fundamentals of Ecosystem Science (Academic Press, 2012).
Google Scholar 
3.Payne, J. A. A summer carrion study of the baby pig Sus scrofa Linnaeus. Ecology 46(5), 592–602 (1965).
Google Scholar 
4.Anderson, G. S., Cervenka, V. J., Haglund, W. & Sorg, M. Insects associated with the body: Their use and analyses. Adv. Forens. Taphonomy 2, 1 (2002).
Google Scholar 
5.Smith, K. G. A manual of forensic entomology. (1986).6.Tomberlin, J. K., Benbow, M. E., Tarone, A. M. & Mohr, R. M. Basic research in evolution and ecology enhances forensics. Trends Ecol. Evol. 26(2), 53–55 (2011).PubMed 

Google Scholar 
7.Benbow, M. E., Tomberlin, J. K. & Tarone, A. M. Carrion Ecology, Evolution, and Their Applications (CRC Press, 2015).
Google Scholar 
8.Wilson, E. E., Mullen, L. M. & Holway, D. A. Life history plasticity magnifies the ecological effects of a social wasp invasion. Proc. Natl. Acad. Sci. 106(31), 12809–12813 (2009).ADS 
CAS 
PubMed 
PubMed Central 

Google Scholar 
9.Pechal, J. L. et al. Field documentation of unusual post-mortem arthropod activity on human remains. J. Med. Entomol. 52(1), 105–108 (2015).PubMed 

Google Scholar 
10.Campobasso, C. P., Marchetti, D., Introna, F. & Colonna, M. F. Postmortem artifacts made by ants and the effect of ant activity on decompositional rates. Am. J. Forens. Med. Pathol. 30(1), 84–87 (2009).
Google Scholar 
11.Eubanks, M. D., Lin, C. & Tarone, A. M. The role of ants in vertebrate carrion decomposition. Food Webs 18, e00109 (2019).
Google Scholar 
12.Cornaby, B. W. Carrion reduction by animals in contrasting tropical habitats. Biotropica 2, 51–63 (1974).
Google Scholar 
13.Andrade-Silva, J., Pereira, E. K. C., Silva, O., Delabie, J. H. C. & Rebelo, J. M. M. Ants (Hymenoptera: Formicidae) associated with pig carcasses in an urban area. Sociobiology 62(4), 527–532 (2015).
Google Scholar 
14.Chin, H. C. et al. Ants (Hymenoptera: Formicidae) associated with pig carcasses in Malaysia. Trop. Biomed. 26(1), 106–109 (2009).
Google Scholar 
15.Prado Castro, C., García, M. D., Palma, C. & Martínez-Ibáñez, M. D. First report on sarcosaprophagous Formicidae from Portugal (Insecta: Hymenoptera). Annales de la Société entomologique de France 50(1), 51–58 (2014).
Google Scholar 
16.Neto-Silva, A., Dinis-Oliveira, R. J. & Prado e Castro, C.,. Diversity of the Formicidae (Hymenoptera) carrion communities in Lisbon (Portugal): Preliminary approach as seasonal and geographic indicators. Forens. Sci. Res. 3(1), 65–73 (2018).
Google Scholar 
17.Payne, J. A., King, E. W. & Beinhart, G. Arthropod succession and decomposition of buried pigs. Nature 219(5159), 1180–1181 (1968).ADS 
CAS 
PubMed 

Google Scholar 
18.Meyer, F., Monroe, M. D., Williams, H. N. & Goddard, J. Solenopsis invicta x richteri (Hymenoptera: Formicidae) necrophagous behavior causes post-mortem lesions in pigs which serve as oviposition sites for Diptera. Forens. Sci. Int. Rep. 2, 100067 (2020).
Google Scholar 
19.De Jong, G. D., Meyer, F. & Goddard, J. Relative roles of blow flies (Diptera: Calliphoridae) and invasive fire ants (Hymenoptera: Formicidae: Solenopsis spp.) in carrion decomposition. J. Med. Entomol. 58(3), 1074–1082 (2021).PubMed 

Google Scholar 
20.Early, M. & Goff, M. L. Arthropod succession patterns in exposed carrion on the island of O’ahu, Hawaiian Islands, USA. J. Med. Entomol. 23(5), 520–531 (1986).CAS 
PubMed 

Google Scholar 
21.Stoker, R. L., Grant, W. E. & Bradleigh Vinson, S. Solenopsis invicta (Hymenoptera: Formicidae) effect on invertebrate decomposers of carrion in central Texas. Environ. Entomol. 24(4), 817–822 (1995).
Google Scholar 
22.Ekanem, M. S. & Dike, M. C. Arthropod succession on pig carcasses in southeastern Nigeria. Papeis Avulsos de Zoologia 50, 561–570 (2010).
Google Scholar 
23.Lindgren, N. K., Bucheli, S. R., Archambeault, A. D. & Bytheway, J. A. Exclusion of forensically important flies due to burying behavior by the red imported fire ant (Solenopsis invicta) in southeast Texas. Forensic Sci. Int. 204(1–3), e1–e3 (2011).PubMed 

Google Scholar 
24.Pereira, E. K. C. et al. Solenopsis saevissima (Smith) (Hymenoptera: Formicidae) activity delays vertebrate carcass decomposition. Sociobiology 64(3), 369–372 (2017).
Google Scholar 
25.Dussutour, A. & Simpson, S. J. Description of a simple synthetic diet for studying nutritional responses in ants. Insectes Soc. 55(3), 329–333 (2008).
Google Scholar 
26.Csata, E. & Dussutour, A. Nutrient regulation in ants (Hymenoptera: Formicidae): A review. Myrmecol. News 29, 111–124 (2019).
Google Scholar 
27.Tschinkel, W. R. The Fire Ants (Belknap Press, 2013).
Google Scholar 
28.Paula, M. C. et al. Action of ants on vertebrate carcasses and blow flies (Calliphoridae). J. Med. Entomol. 53(6), 1283–1291 (2016).PubMed 

Google Scholar 
29.Brooks, M. E. et al. glmmTMB balances speed and flexibility among packages for zero-inflated generalized linear mixed modeling. R J. 9(2), 378–400 (2017).
Google Scholar 
30.Hartig, F. DHARMa: residual diagnostics for hierarchical (multi-level/mixed) regression models. R package version 0.2 4, (2019).31.Fox, J. & Weisberg, S. An R Companion to Applied Regression (Sage publications, 2018).
Google Scholar 
32.Hothorn, T., Bretz, F. & Westfall, P. Simultaneous inference in general parametric models. Biometr. J. 50(3), 346–363 (2008).MathSciNet 
MATH 

Google Scholar 
33.Porter, S. D., Bhatkar, A., Mulder, R., Vinson, B. S. & Clair, D. J. Distribution and density of polygyne fire ants (Hymenoptera: Formicidae) in Texas. J. Econ. Entomol. 84(3), 866–874 (1991).CAS 
PubMed 

Google Scholar 
34.Cook, S. C., Eubanks, M. D., Gold, R. E. & Behmer, S. T. Colony-level macronutrient regulation in ants: mechanisms, hoarding and associated costs. Anim. Behav. 79(2), 429–437 (2010).
Google Scholar 
35.Smith, C. R. & Tschinkel, W. R. Ant fat extraction with a Soxhlet extractor. Cold Spring Harbor Protocols 7, 5243 (2009).
Google Scholar 
36.Wills, B. D. et al. Effect of carbohydrate supplementation on investment into offspring number, size, and condition in a social insect. PLoS ONE 10(7), e0132440 (2015).PubMed 
PubMed Central 

Google Scholar 
37.Bockoven, A. A., Wilder, S. M. & Eubanks, M. D. Intraspecific variation among social insect colonies: persistent regional and colony-level differences in fire ant foraging behavior. PLoS ONE 10(7), e0133868 (2015).PubMed 
PubMed Central 

Google Scholar 
38.Bates, D., Maechler, M., Bolker, B. & Walker, S. Fitting Linear Mixed-Effects Models Using lme4. J. Stat. Softw. 67(1), 1–48 (2015).
Google Scholar 
39.Gavilanez-Slone, J. & Porter, S. D. Colony growth of two species of Solenopsis fire ants (Hymenoptera: Formicidae) reared with crickets and beef liver. Florida Entomol. 96(4), 1482–1488 (2013).
Google Scholar 
40.Sorensen, A. A., Busch, T. M. & Vinson, S. B. Factors affecting brood cannibalism in laboratory colonies of the imported fire ant, Solenopsis invicta Buren (Hymenoptera: Formicidae). J. Kansas Entomol. Soc. 2, 140–150 (1983).
Google Scholar 
41.Williams, D. F., Vander Meer, R. K. & Lofgren, C. S. Diet-induced nonmelanized cuticle in workers of the imported fire ant Solenopsis invicta Buren. Arch. Insect Biochem. Physiol. 4(4), 251–259 (1987).CAS 

Google Scholar 
42.Porter, S. D. Effects of diet on the growth of laboratory fire ant colonies (Hymenoptera: Formicidae). J. Kansas Entomol. Soc. 2, 288–291 (1989).
Google Scholar 
43.Bhatkar, A. & Whitcomb, W. H. Artificial diet for rearing various species of ants. Florida Entomol. 2, 229–232 (1970).
Google Scholar 
44.Porter, S. D., Valles, S. M. & Gavilanez-Slone, J. M. Long-term efficacy of two cricket and two liver diets for rearing laboratory fire ant colonies (Hymenoptera: Formicidae: Solenopsis invicta). Florida Entomol. 98(3), 991–993 (2015).
Google Scholar 
45.Arganda, S. et al. Parsing the life-shortening effects of dietary protein: Effects of individual amino acids. Proc. R. Soc. B Biol. Sci. 284(1846), 20162052 (2017).
Google Scholar 
46.Tschinkel, W. R. Sociometry and sociogenesis of colonies of the fire ant Solenopsis Invicta during one annual cycle. Ecol. Monogr. 63(4), 425–457 (1993).
Google Scholar 
47.Deslippe, R. J. & Savolainen, R. Sex investment in a social insect: The proximate role of food. Ecology 76(2), 375–382 (1995).
Google Scholar 
48.Rosenfeld, C. S. & Roberts, R. M. Maternal diet and other factors affecting offspring sex ratio: A review. Biol. Reprod. 71(4), 1063–1070 (2004).CAS 
PubMed 

Google Scholar 
49.Hasegawa, E. Sex allocation in the ant Camponotus (Colobopsis) nipponicus (Wheeler): II. The effect of resource availability on sex-ratio variability. Insectes Soc. 60(3), 329–335 (2013).
Google Scholar 
50.Knaden, M. & Graham, P. The sensory ecology of ant navigation: from natural environments to neural mechanisms. Annu. Rev. Entomol. 61, 63–76 (2016).CAS 
PubMed 

Google Scholar 
51.Liu, W., Longnecker, M., Tarone, A. M. & Tomberlin, J. K. Responses of Lucilia sericata (Diptera: Calliphoridae) to compounds from microbial decomposition of larval resources. Anim. Behav. 115, 217–225 (2016).
Google Scholar 
52.Tomberlin, J. K. et al. Indole: An evolutionarily conserved influencer of behavior across kingdoms. BioEssays 39(2), 1600203 (2017).
Google Scholar 
53.Frederickx, C. et al. Volatile organic compounds released by blowfly larvae and pupae: New perspectives in forensic entomology. Forensic Sci. Int. 219(1–3), 215–220 (2012).CAS 
PubMed 

Google Scholar 
54.Frederickx, C., Dekeirsschieter, J., Verheggen, F. J. & Haubruge, E. Host-habitat location by the parasitoid, Nasonia vitripennis Walker (Hymenoptera: Pteromalidae). J. Forensic Sci. 59(1), 242–249 (2014).CAS 
PubMed 

Google Scholar 
55.Schettino, M. et al. Response of a predatory ant to volatiles emitted by aphid-and caterpillar-infested cucumber and potato plants. J. Chem. Ecol. 43(10), 1007–1022 (2017).CAS 
PubMed 

Google Scholar 
56.Sawyer, S. J., Rusch, T. W., Tarone, A. M. & Tomberlin, J. K. Wing buzzing as a potential antipredator defense against an invasive predator. Food Webs 27, e00192 (2021).
Google Scholar 
57.Wells, J. D. & Greenberg, B. Effect of the red imported fire ant (Hymenoptera: Formicidae) and carcass type on the daily occurrence of postfeeding carrion-fly larvae (Diptera: Calliphoridae, Sarcophagidae). J. Med. Entomol. 31(1), 171–174 (1994).CAS 
PubMed 

Google Scholar  More

1.FAO. Superfruits: Myth or truth? in Proceedings International Symposium, Ho Chi Minh, Vietnam, 140 (2013).
2.Chamberlain, J., Darr, D. & Meinhold, K. Rediscovering the contributions of forest and trees to transition global food system. Forests 11, 1098. https://doi.org/10.3390/f11101098 (2020).Article 

Google Scholar 
3.Vanzani, P. et al. Wild mediterranean plants as traditional food: A valuable source of antioxidants. J. Food Sci. 76, 46–51 (2011).Article 

Google Scholar 
4.Genskowsky, E. et al. Determination of polyphenolic profile, antioxidant activity and antibacterial properties of maqui [Aristotelia chilensis (Molina) Stuntz] a Chilean blackberry. J. Sci. Food Agric. 96, 4235–4242 (2016).CAS 
Article 

Google Scholar 
5.Benedetti, S. Monografía de maqui, Aristotelia chilensis (Mol.) Stuntz 60 (Instituto Forestal, 2012).
Google Scholar 
6.Vogel, H., Razmilic, H., San Martin, I., Doll, U. & González, B. Plantas Medicinales Chilena. Experiencias de domesticación y cultivo de Boldo, Matico, Bailahuén, Canelo, Peumo y maqui. Editorial Universitaria de Talca, 192 (2005).7.Gironés-Vilaplana, A., Mena, P., García-Viguera, C. & Moreno, D. A novel beverage rich in antioxidant phenolics: Maqui berry (Aristotelia chilensis) and lemon juice. Food Sci. Tech. 47, 279–286 (2012).
Google Scholar 
8.Quispe-Fuentes, I., Vega-Gálvez, A., Vásquez, V., Uribe, E. & Astudillo, S. Mathematical modeling and quality properties of a dehydrated native Chilean berry. J. Food Process Eng. 40, 124–132 (2017).Article 

Google Scholar 
9.Fredes, C., Montenegro, G., Zoffoli, J., Gómez, M. & Robert, P. Polyphenol content and antioxidant activity of maqui during fruit development and maturation in central Chile. Chilean J. Agric. Res. 72, 582–589 (2012).Article 

Google Scholar 
10.Céspedes, C., El-Hafidi, M., Pavon, N. & Alarcon, J. Antioxidant and cardioprotective activities of phenolic extracts from fruits of Chilean blackberry Aristotelia chilensis (Elaeocarpaceae), Maqui. Food Chem. 107, 820–829 (2008).Article 

Google Scholar 
11.Céspedes, C., Alarcon, J., Avila, J. & Nieto, A. Anti-inflammatory activity of Aristotelia chilensis (stuntz) (Elaeocarpaceae). Boletín Latinoamericano y del Caribe de Plantas Medicinales y Aromáticas 9, 91–99 (2010).
Google Scholar 
12.Céspedes, C. et al. The chilean superfruit black-berry Aristotelia chilensis (Elaeocarpaceae), Maqui as mediator in inflammation-associated disorders. Food Chem. Toxicol. 108, 438–450 (2017).
13.Muñoz, O. et al. Chemical study and anti-inflammatory, analgesic and antioxidant activities of the leaves of Aristotelia chilensis (Mol.) Stuntz, Elaeocarpaceae. J. Pharm. Pharmacol. 63, 849–859 (2011).Article 

Google Scholar 
14.Rojo, L. et al. In vitro and in vivo anti-diabetic effects of anthocyanins from maqui berry (Aristotelia chilensis). Food Chem. 131, 387–396 (2012).CAS 
Article 

Google Scholar 
15.Zúñiga, G., Tapia, A., Arenas, A., Contreras, R. & Zuñiga-Libano, G. Phytochemistry and biological properties of Aristotelia chilensis a Chilean blackberry: A review. Phytochem. Rev. 16, 1081–1094. https://doi.org/10.1007/s11101-017-9533-1 (2017).CAS 
Article 

Google Scholar 
16.Vogel, H. et al. Maqui (Aristotelia chilensis): Morpho-phenological characterization to design high-yielding cultivation techniques. J. Appl. Res. Med. Aromat. Plants. 1, 123–133 (2014).
Google Scholar 
17.Liu, Y. & El-Kassaby, Y. Phenotypic plasticity of natural Populus trichocarpa populations in response to temporally environmental change in a common garden. BMC Evol. Biol. 19, 231. https://doi.org/10.1186/s12862-019-1553-6 (2019).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
18.Villemereuil, P., Gaggiotti, O., Mouterde, M. & Till-Bottraud, I. Common garden experiment in the genomic era: New perspectives and opportunities. Heredity 116, 249–254 (2016).Article 

Google Scholar 
19.Torres-Ruiz, J. et al. Genetic differentiation in functional traits among European sessile oak populations. Tree Physiol. 39, 1736–1749. https://doi.org/10.1093/treephys/tpz090 (2019).Article 
PubMed 

Google Scholar 
20.Sáenz-Romero, C., Kremer, A., Nagy, L., Kehlet, J. & Mátyás, C. Common garden comparison confirm inherited differences in sensitivity to climate change between forest tree species. PerrJ. 7, 6213. https://doi.org/10.7717/peerj.6213 (2019).Article 

Google Scholar 
21.Aspinwall, M. et al. Adaptation and acclimation both influence photosynthetic and respiratory temperature responses in Corymbia calophylla. Tree Physiol. 8, 1095–1112. https://doi.org/10.1093/treephys/tpx047 (2017).CAS 
Article 

Google Scholar 
22.Knutzen, F., Meier, I. & Leuschner, C. Does reduced precipitation trigger physiological and morphological drought adaptations in European beech (Fagus sylvatica L.)? Comparing provenances across a precipitation gradient. Tree Physiol. 35, 949–963. https://doi.org/10.1093/treephys/tpv057 (2015).CAS 
Article 
PubMed 

Google Scholar 
23.Mkwezalamba, I., Munthali, C. & Missanjo, E. Phenotypic variation in fruit morphology among provenances of Sclerocarya birrea (A. Rich.) Hochst. Int. J. Forestry Res. 1, 1–8. https://doi.org/10.1155/2015/735418 (2015).Article 

Google Scholar 
24.Sudrajat, D. Genetic variation of fruit, seed, and seedling characteristics among 11 populations of white Jabon in Indonesia. For. Sci. Tech. 12(1), 9–15. https://doi.org/10.1080/21580103.2015.1007896 (2016).Article 

Google Scholar 
25.Teklehaimanot, Z., Lanek, J. & Tomlinson, H. Provenance variation in morphology and leaflet anatomy of Parkia biglobosa and its relation to drought tolerance. Trees 13, 96–102. https://doi.org/10.1007/pl00009742 (1998).Article 

Google Scholar 
26.Åkerström, A., Jaakola, L., Bång, U. & Jäderlund, A. Effects of latitude-related factors and geographical origin on anthocyanidin concentrations in fruits of Vaccinium myrtillus L. (Bilberries). J. Agric. Food Chem. 58, 11939–11945. https://doi.org/10.1021/jf102407n (2010).CAS 
Article 
PubMed 

Google Scholar 
27.Lätti, A., Riihinen, K. & Kainulainen, P. Analysis of anthocyanin variation in wild populations of bilberry (Vaccinium myrtillus L.) in Finland. J. Agric. Food Chem. 56, 190–196. https://doi.org/10.1021/jf072857m (2008).CAS 
Article 
PubMed 

Google Scholar 
28.Uleberg, E. et al. Effects of temperature and photoperiod on yield and chemical composition of Northern and Southern Clones of Bilberry (Vaccinium myrtillus L.). J. Agric. Food Chem. 60, 10406–10414. https://doi.org/10.1021/jf302924m (2012).CAS 
Article 
PubMed 

Google Scholar 
29.Moya, M., González, B., Doll, U., Yuri, J. A. & Vogel, H. Different covers affect growth and development of three maqui clones (Aristotelia chilensis [Molina] Stuntz). J. Berry Res. 1, 1–10. https://doi.org/10.3233/jbr-180377 (2019).Article 

Google Scholar 
30.Cona, M. et al. New polymorphic nuclear microsatellites from Aristotelia chilensis (Mol.) Stuntz (Elaeocarpaceae). Chilean J. Agri. Res. 80, 153–160. https://doi.org/10.4067/S0718-58392020000200153 (2020).Article 

Google Scholar 
31.Hamrick, J. Response of forest trees to global environmental changes. For. Ecol. Manag. 197, 323–335. https://doi.org/10.1016/j.foreco.2004.05.023 (2004).Article 

Google Scholar 
32.Salgado, P., Prinz, K., Finkeldey, R., Ramírez, C. & Vogel, H. Genetic variability of Aristotelia chilensis (“maqui”) based on AFLP and chloroplast microsatellite markers. Gen. Resour. Crop Evol. 64, 2083–2091 (2017).CAS 
Article 

Google Scholar 
33.Holderegger, R., Kamm, U. & Gugerli, F. Adaptive vs. neutral genetic diversity: Implications for landscape genetics. Landsc. Ecol. 21, 797–807. https://doi.org/10.1007/s10980-005-5245-9 (2006).Article 

Google Scholar 
34.O’Brien, E., Mazanex, R. & Krauss, S. Provenance variation of ecologically important traits of forest trees: implications for restoration. J. Appl. Ecol. 44, 583–593. https://doi.org/10.1111/j.1365-2664.2007.01313.x (2007).Article 

Google Scholar 
35.Singleton, V. & Rossi, J. Colorimetry of total phenolics withphosphomolybdic-phosphotungstic acid reagents. Am. J. Enol. Vitic. 16, 144–158 (1965).CAS 

Google Scholar 
36.Giusti, M. & Wrolstad, R. Current protocols in food analytical chemistry. In Current Protocols in Food Analytical Chemistry (eds Wrolstad, R. et al.) F1.2.1-F1.2.13 (Wiley, 2001).
Google Scholar 
37.González, B., Vogel, H., Razmilic, I. & Wolfram, E. Polyphenol, anthocyanin and antioxidant content in different parts of maqui fruits (Aristotelia chilensis) during ripening and conservation treatments after harvest. Ind. Crops Prod. 76, 158–165. https://doi.org/10.1016/j.indcrop.2015.06.038 (2015).CAS 
Article 

Google Scholar 
38.Winn, M., Araman, P. & Lee, S-M. UrbanCrowns: An assessment and monitoring tool for urban trees. Gen. Tech. Rep. SRS-135. Asheville, NC: U.S. Department of Agriculture, Forest Service, Southern Research Station, 10 (2011).39.Welham, S., Cullis, B., Gogel, B., Gilmour, A. & Thompson, R. Prediction in linear mixed models. Aust. N. Z. J. Stat. 46, 325–347 (2004).MathSciNet 
Article 

Google Scholar 
40.Bastías, A. et al. Identification and characterization of microsatellite loci in Maqui (Aristotelia chilensis (Molina) Stuntz) using next-generation sequencing (NGS). PLoS ONE 11, e0159825. https://doi.org/10.1371/journal.pone.0159825 (2016).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
41.Espinoza, S. et al. Influence of provenance origin on the early performance of two sclerophyllous Mediterranean species established in burned drylands. Sci. Rep. 11, 6212. https://doi.org/10.1038/s41598-021-85599-3 (2021).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
42.Vander Mijnsbrugge, K., Bischoff, A. & Smith, B. A question of origin: Where and how to collect seed for ecological restoration. Basic Appl. Ecol. 11, 300–311. https://doi.org/10.1016/j.baae.2009.09.002 (2010).Article 

Google Scholar 
43.Gao, S. B. et al. Phenotypic plasticity vs. local adaptation in quantitative traits differences of Stipa grandis in semi-arid steppe, China. Sci. Rep. 8, 3148. https://doi.org/10.1038/s41598-018-21557-w (2018).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
44.Lusk, C. & Del Pozo, A. Survival and growth of seedlings of 12 Chilean rainforest trees in two light environments: Gas exchange and biomass distribution correlates. Aust. Ecol. 27, 173–182. https://doi.org/10.1046/j.1442-9993.2002.01168.x (2002).Article 

Google Scholar 
45.Brito, C., Bown, H., Fuentes, J., Franck, N. & Perez-Quezada, J. Mesophyll conductance constrains photosynthesis in three common sclerophyllous species in Central Chile. Rev. Chilena de Historia Natural. https://doi.org/10.1186/s40693-014-0008-0 (2014).Article 

Google Scholar 
46.Prado, C. & Damascos, M. Gas exchange and leaf specific mass of different foliar cohorts of the wintergreen shrub Aristotelia chilensis (Mol.) Stuntz (Eleocarpaceae) fifteen days before the flowering and the fall of the old cohort. Braz. Arch. Biol. Tech. 44, 277–282 (2001).Article 

Google Scholar 
47.Repetto-Giavalli, F., Cavieres, L. & Simonetti, J. Respuestas foliares de Aristotelia chilensis (Molina) Stuntz (Elaeocarpaceae) a la fragmentación del bosque maulino. Revista Chilena Hist. Nat. 80, 469–477 (2007).
Google Scholar 
48.Bustan, A. et al. Fruit load governs transpiration of olive trees. Tree Physiol. 36, 380–391. https://doi.org/10.1093/treephys/tpv138 (2016).Article 
PubMed 
PubMed Central 

Google Scholar 
49.Wünsche, J. & Lakso, A. Apple tree physiology—Implications for orchard and tree management. Compact Fruit Tree 33, 82–88 (2000).
Google Scholar 
50.Kelc, D., Vindis, P., Lakota, M. Measurements of Photosynthesis and Transpiration on Apple Trees, Chapter 18 in DAAAM International Scientific Book 2015. in (ed. Katalinic, B.), 199–208. (DAAAM International, 2015). https://doi.org/10.2507/daaam.scibook.2015.18. (ISBN 978-3-902734-05-1, ISSN 1726–9687).51.Lortie, C. & Aarssen, L. The specialization hypothesis for phenotypic plasticity in plants. Int. J. Plant Sci. 157, 484–487. https://doi.org/10.1086/297365 (1996).Article 

Google Scholar 
52.Nemeskéri, E. & Helyes, L. Physiological responses of selected vegetable crop species to water stress. Agronomy 9, 447. https://doi.org/10.3390/agronomy9080447 (2019).CAS 
Article 

Google Scholar 
53.Tian, M., Yu, G., He, N. & Hou, J. Leaf morphological and anatomical traits from tropical to temperate coniferous forests: Mechanisms and influencing factors. Sci. Rep. https://doi.org/10.1038/srep19703 (2016).Article 
PubMed 
PubMed Central 

Google Scholar 
54.Poorter, H., Niinemets, Ü., Poorter, L., Wright, I. J. & Villar, R. Causes and consequences of variation in leaf mass per area (LMA): A meta-analysis. New Phytol. 182, 565–588. https://doi.org/10.1111/j.1469-8137.2009.02830.x (2009).Article 
PubMed 

Google Scholar 
55.Allegro, G., Pastore, C., Valentini, G. & Filippetti, I. The evolution of phenolic compounds in Vitis vinifera L. red berries during ripening: Analysis and role on wine sensory—A review. Agronomy 11, 999. https://doi.org/10.3390/agronomy11050999 (2021).CAS 
Article 

Google Scholar 
56.Chagné, D. et al. Genetic and environmental control of fruit maturation, dry matter and firmness in apple (Malus × domestica Borkh.). Hortic. Res. 1, 14046. https://doi.org/10.1038/hortres.2014.46 (2014).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
57.Gashu, K. et al. Temperature shift between vineyards modulates berry phenology and primary metabolism in a varietal collection of wine grapevine. Front. Plant Sci. 11, 588739. https://doi.org/10.3389/fpls.2020.588739 (2020).Article 
PubMed 
PubMed Central 

Google Scholar 
58.Suter, B., Destrac Irvine, A., Gowdy, M., Dai, Z. & van Leeuwen, C. Adapting wine grape ripening to global change requires a multi-trait approach. Front. Plant Sci. 12, 624867. https://doi.org/10.3389/fpls.2021.624867 (2021).Article 
PubMed 
PubMed Central 

Google Scholar 
59.Nesmith, D. Fruit development period of several Southern Highbush Blueberry Cultivars. Int. J. Fruit Sci. 12, 249–255. https://doi.org/10.1080/15538362.2011.619430 (2012).Article 

Google Scholar 
60.Romero-Román, M. et al. Native species facing climate changes: Response of Calafate Berries To Low Temperature and UV radiation. Foods. 10, 196. https://doi.org/10.3390/foods10010196 (2021).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
61.Cabrera, S., Bozzo, S. & Fuenzalida, H. Variations in UV radiation in Chile. J. Photochem. Photobiol. 28, 137–142 (1995).CAS 
Article 

Google Scholar 
62.Ebel, R. C., Proebsting, E. L. & Evans, R. G. Deficit irrigation to control vegetative growth in apple and monitoring fruit growth to schedule irrigation. HortScience 30, 1229–1232. https://doi.org/10.21273/hortsci.30.6.1229 (1995).Article 

Google Scholar 
63.Fereres, E. & Soriano, M. A. Deficit irrigation for reducing agricultural water use. J. Exp. Bot. 58(2), 147–159. https://doi.org/10.1093/jxb/erl165 (2007).CAS 
Article 
PubMed 

Google Scholar 
64.Barnuud, N., Zerihun, A., Gibberd, M. & Bates, B. Berry composition and climate: Responses and empirical models. Inter. J. Biometeor. 58, 1207–1223. https://doi.org/10.1007/s00484-013-0715-2 (2014).ADS 
Article 

Google Scholar 
65.Spinardi, A., Cola, G., Gardana, C. & Mignani, I. Variation of anthocyanin content and profile throughout fruit development and ripening of highbush blueberry cultivars grown at two different altitudes. Front. Plant Sci. 10, 1045. https://doi.org/10.3389/fpls.2019.01045 (2019).Article 
PubMed 
PubMed Central 

Google Scholar 
66.Stevenson, D. & Scalzo, J. Anthocyanin composition and content of blueberries from around the world. J. Berry Res. 2, 179–189. https://doi.org/10.3233/JBR-2012-038 (2012).CAS 
Article 

Google Scholar 
67.Zarrouk, O. et al. Grape ripening is regulated by deficit irrigation/elevated temperatures according to cluster position in the canopy. Front. Plant Sci. https://doi.org/10.3389/fpls.2016.01640 (2016).Article 
PubMed 
PubMed Central 

Google Scholar 
68.Prange, R. K. & DeEll, J. R. Preharvest factors affecting postharvest quality of berry crops. HortScience 32, 824–830. https://doi.org/10.21273/hortsci.32.5.824 (1997).Article 

Google Scholar 
69.Mignard, O., Beguería, S., Reig, G. & Fonti, C. Genetic origin and climate determine fruit quality and antioxidant traits on apple (Malus × domestica Borkh). Sci. Hortic. 285, 110142. https://doi.org/10.1016/j.scienta.2021.110142 (2021).CAS 
Article 

Google Scholar 
70.González-Villagra, J., Rodrigues-Salvador, A., Nunes-Nesi, A., Cohen, J. & Reyes-Díaz, M. Age-related mechanism and its relationship with secondary metabolism and abscisic acid in Aristotelia chilensis plants subjected to drought stress. Plant Physiol. Biochem. 124, 136–145. https://doi.org/10.1016/j.plaphy.2018.01.010 (2018).CAS 
Article 
PubMed 

Google Scholar 
71.Calderan, A. et al. Managing moderate water deficit increased anthocyanin concentration and proanthocyanidin galloylation in “Refošk” grapes in Northeast Italy. Agric. Water Manage. 246, 106684. https://doi.org/10.1016/j.agwat.2020.106684 (2021).Article 

Google Scholar 
72.Yáñez, M., Seiler, J. & Fox, T. Crown physiological responses of loblolly pine clones and families to silvicultural intensity: Assessing the effect of crown ideotype. For. Ecol. Manage. 398, 25–36. https://doi.org/10.1016/j.foreco.2017.05.002 (2017).Article 

Google Scholar  More

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