Ecology

1.Tapia, J., Audry, S., Townley, B. & Duprey, J. L. Geochemical background, baseline and origin of contaminants from sediments in the mining-impacted Altiplano and Eastern Cordillera of Oruro, Bolivia. Geochemistry 12, 3–20. https://doi.org/10.1144/1467-7873/10-RA-049 (2012).CAS 
Article 

Google Scholar 
2.Sarricolea, E. P. & Romero, H. Variabilidad y cambios climáticos observados y esperados en el Altiplano del norte de Chile. Revista de Geografía Norte Grande. 62, 169–183 (2015).Article 

Google Scholar 
3.Garreaud, R., Vuille, M. & Clement, C. A. The climate of the Altiplano: observed current conditions and mechanisms of past changes. Palaeogeogr. Palaeoclimatol. Palaeoecol. 194(1–3), 5–22. https://doi.org/10.1016/S0031-0182(03)00269-4 (2003).Article 

Google Scholar 
4.Vuille, M. & Keiming, F. Interannual variability of summertime convective cloudiness and precipitation in the central andes derived from ISCCP-B3 data. J. Clim. 17(17), 3334–3348. https://doi.org/10.1175/15200442(2004)017%3c3334:IVOSCC%3e2.0.CO;2 (2004).ADS 
Article 

Google Scholar 
5.Cerveny, R. Present climates of South America. In Climates of the Southern Continents: Present, Past and Future (ed. Hobbs, J. E.) 107–135 (Wiley, Chichester, 1998).
Google Scholar 
6.Garreaud, R. & Aceituno, P. Interannual rainfall variability over the South American Altiplano. J. Clim. 14(12), 2779–2789. https://doi.org/10.1175/1520-0442(2001)014%3c2779:IRVOTS%3e2.0.CO;2 (2001).ADS 
Article 

Google Scholar 
7.Coronel, J., Declerck, S. & Brendonck, L. High-altitude peatland temporary pools in Bolivia house a high cladoceran diversity. Wetlands 27(4), 1166–1174. https://doi.org/10.1672/0277-5212(2007)27[1166:HPTPIB]2.0.CO;2 (2007).Article 

Google Scholar 
8.Dorador, C., Vila, I., Witzel, K. P. & Imhoff, J. F. Bacterial and archaeal diversity in high altitude wetlands of the Chilean Altiplano. Fundam. Appl. Limnol. 182(2), 135–159. https://doi.org/10.1127/1863-9135/2013/0393 (2013).CAS 
Article 

Google Scholar 
9.Garcia, E. & Otto, M. Caracterización ecohidrológica de humedales alto andinos usando imágenes de satélite multitemporales en la cabecera de cuenca del río Santa, Ancash, Perú. Ecología Aplicada 14(2), 115–125 (2013).
Google Scholar 
10.Buytaert, W., Camacho, F. C. & Tobón, C. Potential impacts of climate change on the environmental services of humid tropical alpine regions. Glob. Ecol. Biogeogr. 20(1), 19–33. https://doi.org/10.1111/j.1466-8238.2010.00585.x (2011).Article 

Google Scholar 
11.Hribljan, J. A. et al. Carbon storage and long-term rate of accumulation in high-altitude Andean peatlands of Bolivia. Mires Peat 15(12), 1–14 (2015).
Google Scholar 
12.Yager, K. et al. Dimensiones socioecológicas del cambio del paisaje pastoral andino: puente entre el conocimiento ecológico tradicional y el análisis de imágenes satelitales en Sajama Parque Nacional, Bolivia. Cambio ambiental regional 17, 27–37 (2019).
Google Scholar 
13.Urrutia, R. & Vuille, M. Climate change projections for the tropical Andes using a regional climate model: temperature and precipitation simulations for the end of the 21st century. J. Geophys. Res. 114, D02108. https://doi.org/10.1029/2008JD011021 (2009).ADS 
Article 

Google Scholar 
14.Buytaert, W. et al. Uncertainties in climate change projections and regional downscaling in the tropical Andes: implications for water resources management. Hydrol. Earth. Syst. Sci. 14, 1247–1258. https://doi.org/10.5194/hess-14-1247-2010 (2010).ADS 
Article 

Google Scholar 
15.Rabatel, A. et al. Current state of glaciers in the tropical Andes: a multi-century perspective on glacier evolution and climate change. Cryosphere 7, 81–102. https://doi.org/10.5194/tc-7-81-2013 (2013).ADS 
Article 

Google Scholar 
16.Prieto, G. et al. A mass sacrifice of children and camelids at the Huanchaquito-Las Llamas site, Moche Valley, Peru . PLoS ONE 14(3), e0211691. https://doi.org/10.1371/journal.pone.0211691 (2019).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
17.Babidge, S., Kalazich, F., Prieto, M. & Yager, K. That’s the problem with that lake; it changes sides’: mapping extraction and ecological exhaustion in the Atacama. J. Political Ecol. 26(1), 738–760. https://doi.org/10.2458/v26i1.23169 (2019).Article 

Google Scholar 
18.Prieto, M., Fragkou, M. & Calderón, M. Water policy and management in Chile. In The Wiley Encyclopedia of Water: Science, Technology, and Society (ed. Strickland, C.) 2589–2600 (Wiley, New York, 2020).
Google Scholar 
19.Fritz, S. C., Baker, P. A., Tapia, P., Spanbauer, T. & Westover, K. Evolution of the Lake Titicaca basin and its diatom flora over the last 370,000 years. Palaeogeogr. Palaeoclim. Palaeoecol. 317–318, 93–103 (2012).ADS 
Article 

Google Scholar 
20.Cohen, S. C. Scientific drilling and biological evolution in ancient lakes: lessons learned and recommendations for the future. Hydrobiologia 682(1), 3–25. https://doi.org/10.1007/s10750-010-0546-7 (2012).CAS 
Article 

Google Scholar 
21.Tapia, P. M., Fritz, S. C., Baker, P. A., Seltzer, G. A. & Dunbar, R. B. A Late Quaternary diatom record of tropical climatic history from Lake Titicaca (Peru and Bolivia). Palaeogeogr. Palaeoclimatol. Palaeoecol. 194, 1–3. https://doi.org/10.1016/S0031-0182(03)00275-X (2003).Article 

Google Scholar 
22.Vining, B. R., Steinman, B. A., Abbott, M. B. & Woods, A. Paleoclimatic and archaeological evidence from Lake Suches for highland Andean refugia during the arid middle-Holocene. The Holocene 29(2), 328–344. https://doi.org/10.1177/0959683618810405 (2019).ADS 
Article 

Google Scholar 
23.Fritz, S. C., Baker, P. A., Tapia, P. & Garland, J. Spatial and temporal variation in cores from Lake Titicaca, Bolivia/Peru during the last 13,000 years. Quat. Int. 158(1), 23–29. https://doi.org/10.1016/j.quaint.2006.05.014 (2006).Article 

Google Scholar 
24.Hernández, A. et al. Biogeochemical processes controlling oxygen and carbon isotopes of diatom silica in Late Glacial to Holocene lacustrine rhythmites. Palaeogeogr. Palaeoclimatol. Palaeoecol. 299(3–4), 413–425. https://doi.org/10.1016/j.palaeo.2010.11.020 (2012).Article 

Google Scholar 
25.Placzek, C. et al. Climate in the dry central Andes over Geologic, millennial, and interannual timescales. Ann. Mo. Bot. Gard. 96(3), 386–397. https://doi.org/10.3417/2008019 (2009).Article 

Google Scholar 
26.Cerda, M. et al. A new 20th century lake sedimentary record from the Atacama Desert/Chile reveals persistent PDO (Pacific Decadal Oscillation) impact. J. S. Am. Earth Sci. 95, 102302. https://doi.org/10.1016/j.jsames.2019.102302 (2019).Article 

Google Scholar 
27.Aránguiz-Acuña, A. et al. Aquatic community structure as sentinel of recent environmental changes unraveled from lake sedimentary records from the Atacama Desert, Chile . PLoS ONE 15(2), e0229453. https://doi.org/10.1371/journal.pone.0229453 (2020).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
28.Flores-Varas, A. et al. Ascotán and Carcote salt flats as sensors of humidity fluctuations and anthropic impacts in the transition zone of the Andean Altiplano. J. S. Am. Earth Sci. 105, 102934. https://doi.org/10.1016/j.jsames.2020.102934 (2021).CAS 
Article 

Google Scholar 
29.Maher, B. A. & Taylor, R. M. Formation of ultrafine-grained magnetite in soils. Nature 336, 368–370 (1988).ADS 
CAS 
Article 

Google Scholar 
30.Dearing, J. A. et al. Frequency-dependent susceptibility measurements of environmental materials. Geophys. J. Int. 124, 228–240. https://doi.org/10.1111/j.1365-246X.1996.tb06366.x (1996).ADS 
Article 

Google Scholar 
31.Evans, M. & Heller, F. Environmental magnetism: principles and applications of enviromagnetics. Int. Geophys. 86, 202 (2003).
Google Scholar 
32.Pizarro, H. et al. The origin of the magnetic record in Eocene-Miocene coarse-grained sediments deposited in hyper-arid/arid conditions: examples from the Desert. Palaeogeogr. Palaeoclimatol. Palaeoecol. 516, 322–335. https://doi.org/10.1016/j.palaeo.2018.12.009 (2019).Article 

Google Scholar 
33.Risacher, F., Alonso, H. & Salazar, C. Geoquímica de Aguas en Cuencas cerradas: I, II y III Regiones-Chile. Volumen III Estudio de Cuencas de la II Región. (Convenio de Cooperación DGA-UCN-IRD. S.I.T. Nº51, 1999).34.Tapia, R. et al. Glacial differences of Southern Ocean Intermediate Waters in the Central South Pacific. Quat. Sci. Rev. 208, 105–117. https://doi.org/10.1016/j.quascirev.2019.01.016 (2019).ADS 
Article 

Google Scholar 
35.Horne, D. J. Life-cycles of podocopid Ostracoda – a review (with particular reference to marine and brackish-water species). In Applications of Ostracoda. Proceedings of the Eighth International Symposium on Ostracoda (ed. Maddocks, R.) 581–590 (University of Houston, Texas, 1983).
Google Scholar 
36.Cohen, A. C. & Morin, J. G. Patterns of reproduction in ostracodes; a review. J. Crust. Biol. 10(2), 184–211. https://doi.org/10.2307/1548480 (1990).Article 

Google Scholar 
37.Mesquita-Joanes, F., Smith, A. J. & Viehberg, F. A. The ecology of ostracoda across levels of biological organization from individual to ecosystem. J. Quat. Sci. 17, 15–35. https://doi.org/10.1016/B978-0-444-53636-5.00002-0 (2012).Article 

Google Scholar 
38.McLay, C. L. The population biology of Cyprinotus carolinensis and Herpetocypris reptans (Crustacea, Ostracoda). Can. J. Zool. 56(5), 1170–1179. https://doi.org/10.1139/z78-161 (1978).ADS 
Article 

Google Scholar 
39.Hamouda, S. A., Sames, B., Mohammed, A. & Bensalah, M. First record of non-marine ostracods from the Paleogene “hamadian deposits” of Méridja area, west of Bechar (southwestern Algeria). Annales de Paléontologie 104(1), 27–44. https://doi.org/10.1016/j.annpal.2017.12.001 (2018).Article 

Google Scholar 
40.Bergue, C. T., Maranhao, M. D. S. A. S. & Fauht, G. Paleolimnological inferences based on Oligocene ostracods (Crustacea: Ostracoda) from Tremembé Formation. Southeast Brazil. An. Acad. Bras. Cienc. 87(3), 1531–1544. https://doi.org/10.1590/0001-3765201520140366 (2015).Article 
PubMed 

Google Scholar 
41.Sylvestre, F., Servant-Vildary, S. & Roux, M. Diatom-based ionic concentration and salinity models from the south Bolivian Altiplano (15–23°S). J. Paleolimnol. 25, 279–295 (2001).ADS 
Article 

Google Scholar 
42.Nunnery, J. A., Fritz, S. C., Baker, P. A. & Selenbien, W. Lake-level variability in Salar de Coipasa, Bolivia during the past ∼40,000 yr. Quat. Res. https://doi.org/10.1017/qua.2018.108 (2018).Article 

Google Scholar 
43.Herrera, C. et al. Investigaciones hidrogeológicas en la laguna Tuyajto perteneciente a la Reserva Nacional de los Flamencos (Atacama, Chile). Bol. Geol. Min. 130(4), 789–806. https://doi.org/10.21701/bolgeomin.130.4.011 (2019).Article 

Google Scholar 
44.Houston, J. Variability of precipitation in the Atacama Desert: its causes and hydrological impact. Q. J. R. Meteorol. Soc. 26(15), 2181–2198. https://doi.org/10.1002/joc.1359 (2006).Article 

Google Scholar 
45.Herrera, C. et al. Groundwater flow in a closed basin with a saline shallow lake in a volcanic area: Laguna Tuyajto, northern Chilean Altiplano of the Andes. Sci. Total Environ. 541, 303–318. https://doi.org/10.1016/j.scitotenv.2015.09.060 (2016).ADS 
CAS 
Article 
PubMed 

Google Scholar 
46.Munk, L. A., Boutt, D. F., Hynek, S. A. & Moran, B. J. Hydrogeochemical fluxes and processes contributing to the formation of lithium-enriched brines in a hyper-arid continental basin. Chem. Geol. 493, 37–57. https://doi.org/10.1016/j.chemgeo.2018.05.013 (2018).ADS 
CAS 
Article 

Google Scholar 
47.Godfrey, L. & Álvarez-Amado, F. Volcanic and Saline Lithium Inputs to the Salar de Atacama. Minerals 10(2), 201. https://doi.org/10.3390/min10020201 (2020).CAS 
Article 

Google Scholar 
48.Marazuela, M. A., Ayora, C., Vázquez-Suñé, E., Olivella, S. & García-Gil, A. Hydrogeological constraints for the genesis of the extreme lithium enrichment in the Salar de Atacama (NE Chile): A thermohaline flow modelling approach. Sci. Total Environ. 739, 139959. https://doi.org/10.1016/j.scitotenv.2020.139959 (2020).ADS 
CAS 
Article 
PubMed 

Google Scholar 
49.Bobst, A. L. et al. A 106 ka paleoclimate record from drill core of the Salar de Atacama, northern Chile. Palaeogeogr. Palaeoclimatol. Palaeoecol. 173(1–2), 21–42. https://doi.org/10.1016/S0031-0182(01)00308-X (2001).Article 

Google Scholar 
50.Baspineiro, C. F., Franco, J. & Flexer, V. Potential water recovery during lithium mining from high salinity brines. Sci. Total Environ. 720, 137523. https://doi.org/10.1016/j.scitotenv.2020.137523 (2020).ADS 
CAS 
Article 
PubMed 

Google Scholar 
51.Marazuela, M. A., Vázquez-Suñé, E., Ayora, C. & García-Gil, A. Towards more sustainable brine extraction in salt flats: Learning from the Salar de Atacama. Sci. Total Environ. 703, 135605. https://doi.org/10.1016/j.scitotenv.2019.135605 (2020).ADS 
CAS 
Article 
PubMed 

Google Scholar 
52.Babidge, S. Sustaining ignorance: the uncertainties of groundwater and its extraction in the Salar de Atacama, northern Chile. J. R. Anthropol. Inst. 25(1), 83–102. https://doi.org/10.1111/1467-9655.12965 (2018).Article 

Google Scholar 
53.Sonter, L. J., Dade, M. C., Watson, J. E. M. & Valenta, R. K. Renewable energy production will exacerbate mining threats to biodiversity. Nat. Commun. 11, 4174. https://doi.org/10.1038/s41467-020-17928-5 (2020).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
54.Stahl, A. T., Fremier, A. K. & Cosens, B. A. Mapping legal authority for terrestrial conservation corridors along streams. Conserv. Biol. 34(4), 943–955. https://doi.org/10.1111/cobi.13484 (2020).Article 
PubMed 
PubMed Central 

Google Scholar 
55.García, M., Prieto, M. & Kalazich, F. The protection of the mountain ecosystems of the Southern Central Andes: tensions between Aymara herding practices and conservation policies. Eco. Mont. 13(1), 22–30 (2021).
Google Scholar 
56.Vila, T. Geología de los depósitos salinos andinos, provincia de Antofagasta, Chile. Revista de Geología de Chile 2, 41–55 (1975).
Google Scholar 
57.CIREN. Catastro Agrícola, estudio de metodología para la realización y actualización de catastro agrícola regional en base a la utilización de tecnología geoespacial 1–35 (CCIRA, Atacama, 2013).
Google Scholar 
58.Villagrán, C., Kalin-Arroyo, M. T. & Marticorena, C. Efectos de la destización en la distribución de la flora andina de Chile. Rev. Chil. Hist. Nat. 56, 137–157 (1983).
Google Scholar 
59.CONAF. Actualización Plan de Manejo Participativo Reserva Nacional Los Flamencos, Region de Antofagasta (2008).60.Núñez, L., Grosjean, M. & Cartajena, I. Ocupaciones humanas y paleoambientes en la Puna de Atacama (Universidad Católica del Norte-Taraxacum, Antofagasta, 2005).
Google Scholar 
61.Los Ostracodos, M. P. VI. 4f. In El lago Titicaca, síntesis del conocimiento limnológico actual (eds Dejoux, C. & Iltis, A.) 345–352 (Orstom, New Caledonia, 1991).
Google Scholar 
62.Karanovic, I. Recent Freshwater Ostracods of the World, Crustacea, Ostracoda, Podocopida (Springer , Berlin, 2012).
Google Scholar 
63.Palacios-Fest, M. R., Cusminsky, G. C. & McGlue, M. M. Late Quaternary lacustrine ostracods (Ostracoda, Crustacea) and charophytes (Charophyta, Charales) from the Puna Plateau, Argentina. Micropaleontology 35, 66–78 (2016).
Google Scholar 
64.Brandão, S. N., Angel, M. V., Karanovic, I., Perrier, V. & Meidla, T. World Ostracoda Database. http://www.marinespecies.org/ostracoda/aphia.php?p=taxdetails&id=1091 on 2020–01–15 (2018).65.Fatela, F. & Taborda, R. Confidence limits of species proportions in microfossil assemblages. Mar. Micropaleontol. 45(2), 169–174. https://doi.org/10.1016/S0377-8398(02)00021-X (2002).ADS 
Article 

Google Scholar 
66.Díaz, C. P. & Maidana, N. I. Diatomeas de los Salares Atacama y Punta Negra, II Región-Chile (Centro de Ecología Aplicada , La Reina, 2005).
Google Scholar 
67.Diatoms of North America. The source for diatom identification and ecology. https://diatoms.org (2019).68.Hammer, Ø. & Harper, D. Paleontological Data Analysis (Blackwell Publishing , Hoboken, 2006).
Google Scholar 
69.Oksanen, J. et al. vegan: Community Ecology Package. Ordination methods, diversity analysis and other functions for community and vegetation ecologists. Version 2.5–1. URL https://CRAN.R-project.org/package=vegan (2019).70.McArdle, B. H. & Anderson, M. J. Fitting multivariate models to community data: a comment on distance-based redundancy analysis. Ecology 82(1), 290–297. https://doi.org/10.1890/0012-9658(2001)082[0290:FMMTCD]2.0.CO;2 (2001).Article 

Google Scholar  More

Study areaThe study was conducted along the Agi-gawa River, a tributary of the Kiso-gawa River system in central Japan (35°23 42″–35°26 49″N, 137°25 12″–137°28 01″E; Fig. 1), with the Agi-gawa Dam (110 km from the river mouth, 418 m a.s.l.). The Agi-gawa River is a 3rd to 4th-order river with a naturally sand-rich bed derived from weathered granite that characterizes the local geology36. The Agi-gawa Dam (35°25 32″N, 137°25 55″E) had begun operations in 1990; it is a 102 m high rockfill dam with a catchment area of 82 km2, a storage capacity of 4.8 × 107 m3, a mean depth of ~ 45 to 50 m at the dam site, and a hydraulic residence time of 71 days. Although three small sub-dams at the upstream end of the impoundment trap particulates, the sediment speed in the reservoir has been 1,000,000 m3 for 24 years. The dam serves multiple purposes, including flood control, industrial and urban water supply, and the maintenance of baseflow. Further information on the Agi-gawa Dam is available in Katano et al.37.Figure 1The study area shows six study reaches in three stream segments along the Agi-gawa River and Iinuma-gawa Stream, Gifu Prefecture, Japan. Gray circles denote reaches, which are numbered from upstream to downstream within each segment: UD1 and UD2 are upstream of the dam, DD1 and DD2 are downstream of the dam, whilst TR1 and TR2 are in the tributary. The two black circles denote the sediment replenished reaches (S1 and S2). The three small rectangles at the upstream ends of the impoundment are sub-dams, constructed to reduce the inputs of particulates to the impoundment. This map is based on the Digital Topographic Map 25,000 published by Geospatial Information Authority of Japan.Full size imageSediment replenishment and sampling sitesSediment replenishment was undertaken 0.8 and 1.8 km downstream of the Agi-gawa Dam (S1 and S2, Fig. 1) on February 16 and 27, 2005. A total of 1,200 m3 of sediment (D50 ≈ 0.6 mm; mainly sand) was mined from the upstream sub-dams and transported to S1 and S2. We estimate that this constituted 0.086% of the annual sedimentation in the Agi-gawa Dam (e.g., in 2007, replenished sediment per year × 100/sedimentation in the reservoir). The sediment (800 and 400 m3) was replenished at high-flow banks in both sites. The replenished sediment was gradually washed during the high flows at the end of June (visual observation by dam administrators) (Fig. 2). We confirmed that this replenished sediment remained on both banks in March, and no sediments remained on both banks in early July.Figure 2(a) Precipitation (mm·d) (b) mean inflow to the impoundment per day (m3·s-1); and (c) mean outflow from the Agi-gawa Dam per day (m3·s-1). The vertical broken line indicates the study period. Note that the y-axes for (b) and (c) have a logarithmic scale.Full size imageField sampling was conducted twice between March 15 and 18, 2005, prior to sediment flushing and between August 22 and 24, 2005, following sediment flushing [7 weeks after the end of the sediment drift out (Fig. 2)]. The later sampling date was scheduled to investigate the continuous effects (i.e., not immediate effects) of replenished sediment on the riverbed environment and macroinvertebrate assemblages before the replenished sediment had completely been transported further downstream from S1 and S2.Three study segments (length: 1–2 km each) were selected: (1) upstream of the dam and impounded area (UD); (2) downstream of the dam (DD); and (3) in the tributary (TR). These sites were along a 6.0 km stretch of the Agi-gawa River and a 1.0 km stretch of the Iinuma-gawa Stream (catchment area = 24 km2); the latter is a tributary that flows into the Agi-gawa River 2.7 km, downstream of the dam (Fig. 1, Table 1). Each segment contained two study reaches (six reaches in total), and each study reach was 160 m long with several pool–riffle sequences; all reaches were > 300 m apart. DD1 and DD2 were located immediately downstream of the sediment-displaced banks (S1 and S2; 100 m upstream of DD1 and DD2, respectively). Measurements at the two reaches within the same segment were completed on the same day, and the reaches were surveyed in an upstream direction. The dominant land use along the study area was paddy fields, with sparse riparian forest.Table 1 General characteristics of the three study segments and two seasons.Full size tableAlthough the most suitable reference site for DD is the DD prior to dam construction, we were unable to investigate the site prior to dam construction. Therefore, we treated the reference sites as sites that were less affected by the dam than DD on the present day. Katano et al.37 indicated that the difference between the TR and UD sites was smaller than that between DD and UD/TR sites in terms of biota and geology. However, UD was characterized by a wider channel and higher discharge than TR, due to differences in their catchment areas (Table 1). As we did not have a definitive reference, we treated both UD and TR as reference sites (see “Statistical analysis” section). Therefore, how DD in March and DD in August is different from UD and TR can be interpreted as the effect of sediment reduction.Physical environment and water qualitySix riffles were selected at each study reach, and a sampling location (50 × 50 cm quadrat) was established in the mid-channel area of each riffle. Prior to invertebrate sampling, physical environmental variables were measured.Substrate coarseness was measured by gently floating a Plexiglas observation box (50 × 50 × 10 cm deep) divided into four grid squares (25 × 25 cm) on the surface water such that the grid had projected onto the streambed. The size of the substrate material was coded based on the intermediate-axis length: 1 = sand (particles  16 mm) and sieved through a 0.25 mm mesh sieve. Sieved samples and substrate material smaller than pebbles were mixed in a container and preserved in 5% formalin in the field.The material in each container was later divided into two size fractions using 1-and 0.25 mm mesh sieves. To simplify the sorting process, all material retained in the 0.25 mm sieve was mixed and divided into 2n subsamples (maximum n = 32) using a splitter (Idea Co., Tokyo, Japan), following the method described by Vinson and Hawkins43. All macroinvertebrates in subsamples in the 1 mm sieve were counted and identified to the lowest taxonomic level possible, usually to genus or species level using the taxonomic keys of Kawamura and Ueno44, Merritt and Cummins34, Kathman and Brinkhurst45, Kawai and Tanida35, and Torii46.Macroinvertebrate taxa were also classified into five functional feeding groups (FFGs) according to Kawamura and Ueno44, Merritt and Cummins34, Kathman and Brinkhurst45, Kawai and Tanida (2005)35, and Torii46. FFGs were defined as collector-filterers, collector-gatherers, predators, scrapers, and shredders. If a species belonged to ≥ 2 FFGs, the number of individuals was apportioned across the FFGs. We also counted the number of burrowers (#burrowers), inorganic case-bearing caddisflies (#ICB), and net-spinners (#net spinners) of macroinvertebrate assemblages according Kawamura and Ueno44, Merritt and Cummins34, Kathman and Brinkhurst45, Kawai and Tanida35, and Torii46 (see Supplementary Table S1). This classification was carried out as such life-habit traits are important for surviving in a regulated river containing reduced quantities of sand and gravel on the riverbed37. The Chironomidae family was excluded in the life-habit analysis as they consist of various life forms. Once all invertebrates were removed, dry mass (mg m−2) and ash-free dry mass (AFDM, mg m−2) of benthic coarse particulate organic matter (BCPOM,  > 1 mm), and benthic fine particulate organic matter (BFPOM,  0.25 mm) were obtained by drying in an oven at 60 °C for 1 day and combusting in a muffle furnace at 550 °C for 4 h. BCPOM and BFPOM were calculated based on the difference between the dry mass and the AFDM.The total number of invertebrate individuals and the AFDM of BFPOM in each sample were estimated by multiplying by the corresponding 2n value. The number of taxa and density of invertebrates in each sample were calculated as the sum of the values in both size fractions. Additionally, we determined Shannon’s diversity index (H), Simpson’s evenness index, and the percentage of Ephemeroptera, Plecoptera, and Trichoptera (%EPT)47. A sample from UD2 in March had been lost and therefore could not be included in the analyses.Periphyton was sampled from cobbles adjacent to each sampling location. Periphyton was removed from a 5 × 5 cm area on the upper surface of each cobble with a toothbrush. Each sample was placed in a separate container with 250 mL of water. Within 24 h of sample collection, a subsample of the well-mixed content in each container was filtered using a glass-fiber filter (GF/C; Whatman Co., Maidstone, UK). Each filter was placed in a separate vial with 20 mL of 99.5% ethanol and stored in a dark refrigerator at 4 °C for 24 h. The extracted pigments were measured using a spectrophotometer (U-1800; Shimadzu Co., Kyoto, Japan), following the method of Lorenzen48.Analysis of case materials of an inorganic case-bearing caddisflyWe compared the particle size structure of replenished sediment, riverbed sediment, and case materials for case-bearing caddisfly. The replenished sediment was directly sampled in a 1 L polyethylene jar at the upstream replenished bank (S1) on March 16, 2005 (Fig. 1). Riverbed sediment was sampled at two stations; 100 m upstream of S1, and 100 m upstream of DD1 between August 22 and 24, 2005. At each station of the river, a metallic narrow cup (200 mL) with a lid was pushed into a vacancy between the cobbles, which had been randomly selected, and fine sediments (up to small gravel) in the vacancies were sampled by closing the lid underwater. Sampling was carried out three times (i.e., three different vacancies in the cobbles), and subsamples were pooled for measurement. The replenished and riverbed sediment was combusted at 550 °C for 2 h in a muffle furnace to remove organic contamination. Combusted samples were separated with eight sieves with a mesh size range of 0.075–9.5 mm (JIS A 1204). Each fraction was weighed, and the grain size accumulation curve of each type of sediment and its D50 were obtained.In a macroinvertebrate sample at DD1 between August 22 and 24, 2005, ten individuals from two case-bearing caddisfly larvae, Glossosoma sp. and Gumaga orientalis, which were prevalent at DD1 during this period (see Results), were randomly selected from the formalin-fixed sample. The case was carefully removed from the larvae and combusted as described above for the replenished and riverbed sediment. The number of case material grains was measured using a dissection microscope.Statistical analysesWe described results based on two main assumptions: (1) the DD in March is the dam-affected reach (cf. unregulated reaches UD and TR), and (2) the changes in DD from March to August were mainly a result of sediment replenishment. In the statistical analyses, the p criterion (⍺) was set at 0.05.To consider the effects of the segment, replicate reach, and season on variables, nested multivariate analysis of variance (MANOVA) was used to test whether any measured variables at the riffle scale differed between segments (UD, DD, and TR). Three segments and two replicate reaches were nested within each season (March and August) and segment (i.e., Season/Segment/Reach), whereby measurements within each reach were treated as subsamples. In the MANOVA, we also consider the interactions of the variances to interpret the interactions among the sampling segments and seasons to consider the independent effects on the factors.To perform MANOVA, we assumed that temporal variability was greater than spatial variability within each reach for variables measured over 24 h (e.g., water quality), and the opposite would hold true for variables measured only once (e.g., macroinvertebrates). Therefore, subsamples within each reach were either spatially or temporally replicated, depending on the variable type. Temporal replicates (four samples collected every 6 h) were treated as a repeated factor (time factor). A nested MANOVA was used for variables quantified once at each location (e.g., macroinvertebrates), and nested repeated-measures MANOVA (rm-MANOVA) were used for variables quantified over a 24 h period at each reach (e.g., water quality). When a significant difference was detected by MANOVA with non-significant interactions, each variable was tested separately with a nested ANOVA for variable groups once at each location or the nested rm-ANOVA for repeated-measured variables, as appropriate for the particular variable. The risk of inflating Type 1 errors for the ANOVA was reduced using Bonferroni adjustments.These MANOVA and ANOVA tests were conducted with R version 3.6.049. The residuals of each variable in each MANOVA and ANOVA model were verified using the Shapiro–Wilk normality test prior to analyses, and normality was improved using arcsine(x) or log (x + 1) transformation when appropriate.Tukey’s multiple comparison test in a one-way ANOVA model (Season/Segment/Reach) was used for comparisons between segments. Any significant changes in values for variables from UD to DD were interpreted as the effects of the dam based on the assumption that conditions in UD and DD were similar prior to dam construction; this was because replenished sediment had not been supplied in March (see before). However, UD may be unsuitable as a reference site compared with TR as the former may be at least partly affected by the dam. This may particularly be the case for benthic invertebrates, such as the interruption of the upstream flight of adult females50. Therefore, UD and TR were treated as reference sites for reservoir and tributary effects, respectively. This was because both were unaffected by the dam, and sediment replenishment as tributaries may function as sites for resource recovery for the dam-affected mainstem of the river37,51,52, despite differing watershed areas. Therefore, the similarity of variables between the TR and UD sites was statistically confirmed such that they could be treated as reference sites. As such, the recovery from March to August could reliably demonstrate the effect of sediment replenishment. For example, although the value at DD differed from that at TR and/or UD in March, it was similar to that at UD and/or TR in August.Multivariate analyses were conducted using the R “vegan” package version 2.5.6 to compare invertebrate assemblage structures between segments. Bray–Curtis coefficients based on species abundance were used to calculate a dissimilarity matrix, and dissimilarities between UD and DD, and between TR and DD in each season were tested using two-way ANOVA and Tukey post-hoc tests.Macroinvertebrate assemblage organization in relation to environmental gradients was analyzed using redundancy analysis (RDA) with the “rda” function of “vegan” package. This was because the preliminary analysis using detrended correspondence analysis (DCA) showed that the gradient lengths of DCA were More

Experimental study and protein domain analysisInsectsHosts and parasitoids were maintained as previously described25. Five Leptopilina heterotoma (Hymenoptera: Figitidae) populations were used for experiments: a population from Japan (Sapporo), two populations from the United Kingdom (1: Whittlesford; 2: Great Shelford) and two populations from Belgium (1: Wilsele; 2: Eupen). Information on collection sites, including GPS coordinates, can be found in25.Determination of host fat contentD. simulans and D. melanogaster hosts were allowed to lay eggs during 24 h in glass flasks containing ~ 50 mL standard medium25. After two days, developing larvae were sieved and ~ 200 were larvae placed in a Drosophila tube containing ~ 10 mL medium. Seven days after egg laying, newly formed pupae were frozen at – 18 °C, after which fat content was determined as described in25, where dry weight before and after neutral fat extraction was used to calculate absolute fat amount (in μg) for each host. The host pupal stage was chosen for estimating fat content, because at this point the host ceases to feed, while the parasitoid starts consuming the entire host36. All data were analysed using R Project version 3.4.360. Fat content of hosts was compared using a one-way ANOVA with host species as fixed factor.Manipulation of host fat contentTo generate leaner D. melanogaster hosts, we adapted our standard food medium25 to contain 100 times less (0.5 g) sugar per litre water. Manipulating sugar content did not alter the structure of the food medium, thus maintaining similar rearing conditions, with the exception of sugar content. Fat content of leaner and fatter D. melanogaster hosts was determined and analysed as described above.Fat synthesis quantification of wasp populationsMated female L. heterotoma were allowed to lay eggs on host fly larvae collected as described above with ad libitum access to honey as a food source until death. Honey consists of sugars and other carbohydrates that readily induce fat synthesis. After three weeks, adult offspring emergence was monitored daily and females were haphazardly placed in experimental treatments. Females were either killed at emergence (to measure teneral lipid reserves) or after feeding for 7 days on honey. Wasps were frozen at − 18 °C after completion of experiments. Fat content was determined as described above for hosts. The ability for fat synthesis was then determined by comparing mean fat levels of recently emerged compared to fed individuals, similar to procedures described in15,25,28. An increase in fat levels after feeding is indicative of active fat synthesis; equal or lower fat levels suggest fat synthesis did not take place. Each population tested on D. melanogaster or D. simulans represented an independent dataset that was analysed separately, as in Visser et al. 201825, because we are interested in the response of each population on each host species. We used T-tests when data was normally distributed and variances equal, log-transformed data for non-normal data, and a Welch’s t-test when variances were unequal. We corrected for multiple testing using Benjamini and Hochberg’s False Discovery Rate61.Fat synthesis quantification using a familial design and GC–MS analysesTo tease apart the effect of wasp genotype and host environment, we used a split-brood design where the offspring of each mother developed on lean D. simulans or fat D. melanogaster hosts in two replicated experiments (experiment 1 and 2). In both experiments, mothers were allowed to lay eggs in ~ 200 2nd to 3rd instar host larvae of one species for four days, after which ~ 200 host larvae of the other species were offered during four days. The order in which host larvae were presented was randomized across families. Following offspring emergence, daughters were allocated into two treatment groups: a control where females were fed a mixture of honey and water (1:2 w/w) or a treatment group fed a mixture of honey and deuterated water (Sigma Aldrich) (1:2 w/w; stable isotope treatment) for 7 days. Samples were prepared for GC–MS as described in 28. Incorporation of up to three deuterium atoms can be detected, but percent incorporation is highest when only 1 deuterium atom is incorporated. As incorporation of a single atom unequivocally demonstrates active fat synthesis, we only analysed percent incorporation (in relation to the parent ion) for the abundance of the m + 1 ion. Percent incorporation was determined for five fatty acids, C16:1 (palmitoleic acid), C16:0 (palmitate), C18:2 (linoleic acid), C18:1 (oleic acid), and C18:0 (stearic acid), and the internal standard C17:0 (margaric acid). Average percent incorporation for C17:0 was 19.4 (i.e. baseline incorporation of naturally occurring deuterium) and all values of the internal standard remained within 3 standard deviations of the mean (i.e. 1.6). Percent incorporation of control samples was subtracted from treatment sample values to correct for background levels of deuterium (i.e. only when more deuterium is incorporated in treatment compared to controls fatty acids are actively being synthesized). For statistical analyses, percent incorporation was first summed for C16:1, C16:0, C18:2, C18:1 and C18:0 to obtain overall incorporation levels, as saturated C16 and C18 fatty acids are direct products of the fatty acid synthesis pathway (that can subsequently be desaturated).Data (presented in Fig. 1) was analysed by means of a linear mixed effects model (GLMM, lme4 package) with host (lean D. simulans and fat D. melanogaster) and experiment (conducted twice) as fixed effect, family nested within population (Japan, United Kingdom 1 and 2, Belgium 1 and 2) as random factor, and percentage of incorporation of stable isotopes as dependent variable (log transformed; n = 138). Non-significant terms (i.e., experiment) were sequentially removed from the model to obtain the minimal adequate model as reported in Table 2. When referring to “families,” we are referring to the comparison of daughters of singly inseminated females, which (in these haplodiploid insects) share 75% of their genome.Identification of functional acc and fas genes in distinct parasitoid speciesTo obtain acc and fas nucleotide sequences for L. clavipes, G. legneri, P. maculata and A. bilineata, we used D. melanogaster mRNA ACC transcript variant A (NM_136498.3 in Genbank) and FASN1-RA (FBtr0077659 in FlyBase) and blasted both sequences against transcripts of each parasitoid (using the blast function available at http://www.parasitoids.labs.vu.nl62,63). Each nucleotide sequence was then entered in the NCBI Conserved Domain database64 to determine the presence of all functional protein domains. All sequences were then translated using the Expasy translate tool (https://web.expasy.org/translate/), where the largest open reading frame was selected for further use and confirming no stop codons were present. Protein sequences were then aligned using MAFFT v. 7 to compare functional amino acid sequences between all species (Supplementary files 1 and 2)65.Simulation studyWe consider the general situation where phenotypic plasticity is only sporadically adaptive and ask the question whether and under what circumstances plasticity can remain functional over long evolutionary time periods when the regulatory processes underlying plasticity are gradually broken down by mutations. We consider a regulatory mechanism that switches on or off a pathway (like fat synthesis) in response to environmental conditions (e.g., host fat content).Fitness considerationsWe assume that the local environment of an individual is characterized by two factors: fat content F and nutrient content N, where nutrients represent sugars and other carbohydrates that can be used to synthesize fat. Nutrients are measured in units corresponding to the amount of fat that can be synthesized from them. We assume that fitness (viability and/or fecundity) is directly proportional to the amount of fat stored by the individual. When fat synthesis is switched off, this amount is equal to F, the amount of fat in the environment. When fat synthesis is switched on, the amount of fat stored is assumed to be (N – c + (1 – k)F). This expression reflects the following assumptions: (i) fat is synthesized from the available nutrients, but this comes at a fitness cost c; (ii) fat can still be absorbed from the environment, but at a reduced rate ((1 – k)). It is adaptive to switch on fat synthesis if (N – c + (1 – k)F) is larger than F, or equivalently if (F < tfrac{1}{k}(N - c)).The right-hand side of this inequality is a straight line, which is illustrated by the blue line in Fig. 4. The three boxes in Fig. 4 illustrate three types of environmental conditions. Red box low-fat environments. Here, (F < tfrac{1}{k}(N - c)) is always satisfied, implying that fat synthesis should be switched on constitutively. Yellow box high-fat environments. Here, (F > tfrac{1}{k}(N – c)), implying that fat synthesis should be switched off constitutively.

Orange box intermediate-fat environments. Here, fat synthesis should be plastic and switched on if for the given environment (N, F) the fat content is below the blue line and switched off otherwise.

Figure 4Environmental conditions encountered by the model organisms. For a given combination of environmental nutrient content N and environmental fat content F, it is adaptive to switch on fat synthesis if (N, F) is below the blue line (corresponding to (F < tfrac{1}{k}(N - c))) and to switch it off otherwise. The three boxes illustrate three types of environment: a low-fat environment (red) where fat synthesis should be switched on constitutively; a high-fat environment (yellow) where fat synthesis should be switched off constitutively; and an intermediate-fat environment (orange) where a plastic switch is selectively favoured.Full size imageThe simulations reported here were all run for the parameters (k = tfrac{1}{2}{text{ and }}c = tfrac{1}{4}). We also investigated many other combinations of these parameters; in all cases, the results were very similar to those reported in Fig. 3.Gene regulatory networks (GRN)In our model, the switching device was implemented by an evolving gene regulatory network (as in van Gestel and Weissing66). The simulations shown in Fig. 3 of the main text are based on the simplest possible network that consists of two receptor nodes (sensing the fat and the nutrient content in the local environment, respectively) and an effector node that switches on fat synthesis if the combined weighted input of the two receptor nodes exceeds a threshold value T and switches it off otherwise. Hence, fat synthesis is switched on if (w_{F} F + w_{N} N > T) (and off otherwise). The GRN is characterized by the weighing factors (w_{F} {text{ and }}w_{N}) and the threshold T. These parameters are transmitted from parents to offspring, and they evolve subject to mutation and selection. We also considered alternative network structures (all with two receptor nodes and one effector node, but with a larger number of evolvable weighing factors67, and obtained very similar results, see below).For the simple GRN described above, the switching device is 100% adaptive when the switch is on (i.e., (w_{F} F + w_{N} N > T)) if (F < tfrac{1}{k}(N - c)) and off otherwise. A simple calculation yields that this is the case if: (w_{N} > 0{, }w_{F} = – k{kern 1pt} w_{N} {text{ and }}T = c{kern 1pt} w_{N}).Evolution of the GRNFor simplicity, we consider an asexual haploid population with discrete, non-overlapping generations and fixed population size (N = 10,000). Each individual has several gene loci, each locus encoding one parameter of the GRN. In case of the simple network described above, there are three gene loci, each with infinitely many alleles. Each individual harbours three alleles, which correspond to the GRN parameters (w_{F} {, }w_{N} {text{ and }}T), and hence determine the functioning of the genetic switch. In the simulations, each individual encounters a randomly chosen environment ((N{, }F)). Based on its (genetically encoded) GRN, the individual decides on whether to switch on or off fat synthesis. If synthesis is switched on, the individual’s fitness is given by (N – c + (1 – k)F); otherwise its fitness is given by F. Subsequently, the individuals produce offspring, where the number of offspring produced is proportional to the amount of fat stored by an individual. Each offspring inherits the genetic parameters of its parent, subject to mutation. With probability μ (per locus) a mutation occurs. In such a case the parental value (in case of a simple network: the parent’s allelic value (w_{F} {, }w_{N} {text{ or }}T)) is changed to a mutated value ((w_{F} { + }delta {, }w_{N} { + }delta {text{ or }}T + delta)), where the mutational step size δ is drawn from a normal distribution with mean zero and standard deviation σ. In the reported simulations, we chose (mu = 0.001) and (sigma = 0.1). The speed of evolution is proportional to (mu cdot sigma^{2}), implying that the rate of change in Fig. 3 (both the decay of plasticity and the rate of regaining adaptive plasticity) are positively related to μ and σ.Preadaptation of the GRNsStarting with a population with randomly initialized alleles for the GRN parameters, we first let the population evolve for 10,000 generations in the intermediate-fat environment (the orange box in Fig. 4). In all replicate simulations, a “perfectly adapted switch” (corresponding to (w_{N} > 0{, }w_{F} = – k{kern 1pt} w_{N} {text{ and }}T = c{kern 1pt} w_{N})) evolved, typically within 1,000 generations. Still, the evolved GRNs differed across replicates, as they evolved different values of (w_{N} > 0). These evolved networks were used to seed the populations in the subsequent “decay” simulations.Evolutionary decay of the GRNsFor the decay experiments reported in Fig. 3 of the main text, we initiated a large number of monomorphic replicate populations with one of the perfectly adapted GRNs from the preadaptation phase. These populations were exposed for an extended period of time (1,000,000 generations) to a high-fat environment (the yellow box in Fig. 4), where all preadapted GRNs switched off fat synthesis. However, in some scenarios, the environmental conditions changed back sporadically (with probability q) to the intermediate-fat environment (the orange box in Fig. 4), where it is adaptive to switch on fat metabolism in 50% of the environmental conditions (when (N, F) is below the blue line in Fig. 4). In Fig. 3, we report on the changing rates (q = 0.0) (no changing back; red), (q = 0.001) (changing back once every 1,000 generations; purple), and (q = 0.01) (changing back once every 100 generations; pink). When such a change occurred, the population was exposed to the intermediate-fat environment for t generations (Fig. 3 is based on t = 3).Throughout the simulation, the performance of the network was monitored every 100 generations as follows: 100 GRNs were chosen at random from the population, and each of these GRNs was exposed to 100 randomly chosen environmental conditions from the intermediate-fat environment (orange box in Fig. 4). From this, we could determine the average percentage of “correct” decisions (where the network should be switched on if and only if (F < tfrac{1}{k}(N - c)). 1.0 means that the GRN is still making 100% adaptive decisions; 0.5 means that the GRN only makes 50% adaptive decision, as would be expected by a random GRN or a GRN that switches the pathway constitutively on or off. This measure for performance in the “old” intermediate-fat environment was determined for 100 replicate simulations per scenario and plotted in Fig. 3 (mean ± standard deviation).Evolving robustness of the GRNsThe simulations in Fig. 3 are representative for all networks and parameters considered. Whenever (q = 0.0), the performance of the regulatory switch eroded in evolutionary time, but typically at a much lower rate in case of the more complex GRNs. Whenever (q = 0.01), the performance of the switch went back to levels above 90% and even above 95% for the more complex GRNs. Even for (q = 0.001), a sustained performance level above 75% was obtained in all cases.Intriguingly, in the last two scenarios the performance level first drops rapidly (from 1.0 to a much lower level, although this drop is less pronounced in the more complex GRNs) and subsequently recovers to reach high levels again. Apparently, the GRNs have evolved a higher level of robustness, a property that seems to be typical for evolving networks8. For the simple GRN studied in Fig. 3, this outcome can be explained as follows. The initial network was characterized by the genetic parameters (w_{N} > 0{, }w_{F} = – k{kern 1pt} w_{N} {text{ and }}T = c{kern 1pt} w_{N}) (see above), where (w_{N}) was typically a small positive number. In the course of evolutionary time, the relation between the three evolving parameters remained approximately the same, but (w_{N}) (and with it the other parameters) evolved to much larger values. This automatically resulted in an increasingly robust network, since mutations with a given step size distribution affect the performance of a network much less when the corresponding parameter is large in absolute value.Costs of plasticityPhenotypically plastic organisms can incur different types of costs68. In our simple model, we only consider the cost of phenotype-environment mismatching, that is, the costs of expressing the ‘wrong’ phenotype in a given environment. When placed in a high-fat environment, the preadapted GRNs in our simulations take the ‘right’ decision to switch off fat metabolism. Accordingly, they do not face any costs of mismatching. Yet, the genetic switch rapidly decays (as indicated in Fig. 3 by the rapid drop in performance when tested in an intermediate-fat environment), due to the accumulation of mutations.It is not unlikely that there are additional fitness costs of plasticity, such as the costs for the production and maintenance of the machinery underlying plasticity68. In the presence of such constitutive costs, plasticity will be selected against when organisms are living in an environment where only one phenotype is optimal (as in the high- and low-fat environments in Fig. 4). This would obviously affect the evolutionary dynamics in Fig. 3, but the size of the effect is difficult to judge, as the constitutive costs of plasticity are notoriously difficult to quantify. In case of the simple switching device considered in our model, we consider the constitutive costs of plasticity as marginal, but these costs might be substantial in other scenarios. More

1.Conant, R. T. et al. Temperature and soil organic matter decomposition rates—synthesis of current knowledge and a way forward. Glob. Change Biol. 17, 3392–3404 (2011).ADS 
Article 

Google Scholar 
2.Yost, J. L. & Hartemink, A. E. How deep is the soil studied—an analysis of four soil science journals. Plant Soil 452, 5–18 (2020).CAS 
Article 

Google Scholar 
3.Jobbágy, E. G. & Jackson, R. B. The vertical distribution of soil organic carbon and its relation to climate and vegetation. Ecol. Appl. 10, 423–436 (2000).Article 

Google Scholar 
4.Dove, N. C. et al. Continental-scale patterns of extracellular enzyme activity in the subsoil: an overlooked reservoir of microbial activity. Environ. Res. Lett. 15, 1040a1 (2020).CAS 
Article 

Google Scholar 
5.Hicks Pries, C. E., Castanha, C., Porras, R. C. & Torn, M. S. The whole-soil carbon flux in response to warming. Science 355, 1420–1423 (2017).ADS 
CAS 
PubMed 
Article 

Google Scholar 
6.Taş, N. et al. Impact of fire on active layer and permafrost microbial communities and metagenomes in an upland Alaskan boreal forest. ISME J. 8, 1904–1919 (2014).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
7.Brewer, T. E. et al. Ecological and genomic attributes of novel bacterial taxa that thrive in subsurface soil horizons. mBio 10, e01318–e01319 (2019).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
8.Sulman, B. N. et al. Multiple models and experiments underscore large uncertainty in soil carbon dynamics. Biogeochemistry 141, 109–123 (2018).CAS 
Article 

Google Scholar 
9.Romero-Olivares, A. L., Allison, S. D. & Treseder, K. K. Soil microbes and their response to experimental warming over time: a meta-analysis of field studies. Soil Biol. Biochem. 107, 32–40 (2017).CAS 
Article 

Google Scholar 
10.Melillo, J. M. et al. Long-term pattern and magnitude of soil carbon feedback to the climate system in a warming world. Science 358, 101–105 (2017).ADS 
CAS 
PubMed 
Article 

Google Scholar 
11.Melillo, J. M. et al. Soil warming, carbon–nitrogen interactions, and forest carbon budgets. Proc. Natl Acad. Sci. USA 108, 9508–9512 (2011).ADS 
CAS 
PubMed 
Article 

Google Scholar 
12.Crowther, T. W. et al. Quantifying global soil carbon losses in response to warming. Nature 540, 104–108 (2016).ADS 
CAS 
PubMed 
Article 

Google Scholar 
13.Feng, X., Simpson, A. J., Wilson, K. P., Dudley Williams, D. & Simpson, M. J. Increased cuticular carbon sequestration and lignin oxidation in response to soil warming. Nat. Geosci. 1, 836–839 (2008).ADS 
CAS 
Article 

Google Scholar 
14.Dove, N. C., Stark, J. M., Newman, G. S. & Hart, S. C. Carbon control on terrestrial ecosystem function across contrasting site productivities: the carbon connection revisited. Ecology 100, e02695 (2019).PubMed 
Article 

Google Scholar 
15.Frey, S. D., Lee, J., Melillo, J. M. & Six, J. The temperature response of soil microbial efficiency and its feedback to climate. Nat. Clim. Change 3, nclimate1796 (2013).Article 
CAS 

Google Scholar 
16.Olander, L. P. & Vitousek, P. M. Regulation of soil phosphatase and chitinase activity by N and P availability. Biogeochemistry 49, 175–190 (2000).CAS 
Article 

Google Scholar 
17.DeAngelis, K. M. et al. Long-term forest soil warming alters microbial communities in temperate forest soils. Front. Microbiol. 6, 104–115 (2015).18.Cheng, L. et al. Warming enhances old organic carbon decomposition through altering functional microbial communities. ISME J. 11, 1825–1835 (2017).PubMed 
PubMed Central 
Article 

Google Scholar 
19.Billings, S. A. & Ballantyne, F. How interactions between microbial resource demands, soil organic matter stoichiometry, and substrate reactivity determine the direction and magnitude of soil respiratory responses to warming. Glob. Change Biol. 19, 90–102 (2013).ADS 
Article 

Google Scholar 
20.Rasmussen, C., Torn, M. S. & Southard, R. J. Mineral assemblage and aggregates control carbon dynamics in a California conifer forest. Soil Sci. Soc. Am. J. 69, 1711–1721 (2005).ADS 
CAS 
Article 

Google Scholar 
21.Jones, D. L. et al. Microbial competition for nitrogen and carbon is as intense in the subsoil as in the topsoil. Soil Biol. Biochem. 117, 72–82 (2018).CAS 
Article 

Google Scholar 
22.Pold, G. et al. Carbon use efficiency and its temperature sensitivity covary in soil bacteria. mBio 11, e02293–19 (2020).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
23.Geyer, K. M., Dijkstra, P., Sinsabaugh, R. & Frey, S. D. Clarifying the interpretation of carbon use efficiency in soil through methods comparison. Soil Biol. Biochem. 128, 79–88 (2018).24.Schmidt, M. W. I. et al. Persistence of soil organic matter as an ecosystem property. Nature 478, 49–56 (2011).ADS 
CAS 
PubMed 
Article 

Google Scholar 
25.Vitousek, P. M., Porder, S., Houlton, B. Z. & Chadwick, O. A. Terrestrial phosphorus limitation: mechanisms, implications, and nitrogen–phosphorus interactions. Ecol. Appl. 20, 5–15 (2010).PubMed 
Article 

Google Scholar 
26.Sullivan, B. W. et al. Assessing nutrient limitation in complex forested ecosystems: alternatives to large‐scale fertilization experiments. Ecology 95, 668–681 (2014).PubMed 
Article 

Google Scholar 
27.Hart, S. C., Firestone, M. K. & Paul, E. A. Decomposition and nutrient dynamics of ponderosa pine needles in a Mediterranean-type climate. Can. J. Forest Res. 22, 306–314 (1992).CAS 
Article 

Google Scholar 
28.Dijkstra, P. et al. Effect of temperature on metabolic activity of intact microbial communities: evidence for altered metabolic pathway activity but not for increased maintenance respiration and reduced carbon use efficiency. Soil Biol. Biochem. 43, 2023–2031 (2011).CAS 
Article 

Google Scholar 
29.Don, A., Rödenbeck, C. & Gleixner, G. Unexpected control of soil carbon turnover by soil carbon concentration. Environ. Chem. Lett. 11, 407–413 (2013).CAS 
Article 

Google Scholar 
30.Rodriguez-R, L. M. & Konstantinidis, K. T. Nonpareil: a redundancy-based approach to assess the level of coverage in metagenomic datasets. Bioinformatics 30, 629–635 (2014).CAS 
PubMed 
Article 

Google Scholar 
31.Rodriguez-R, L. M. & Konstantinidis, K. T. Estimating coverage in metagenomic data sets and why it matters. ISME J. 8, 2349–2351 (2014).PubMed 
PubMed Central 
Article 

Google Scholar 
32.Lombard, V., Golaconda Ramulu, H., Drula, E., Coutinho, P. M. & Henrissat, B. The carbohydrate-active enzymes database (CAZy) in 2013. Nucleic Acids Res. 42, D490–D495 (2014).CAS 
Article 

Google Scholar 
33.Levasseur, A., Drula, E., Lombard, V., Coutinho, P. M. & Henrissat, B. Expansion of the enzymatic repertoire of the CAZy database to integrate auxiliary redox enzymes. Biotechnol. Biofuels 6, 41 (2013).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
34.Benoit, I. et al. Degradation of different pectins by fungi: correlations and contrasts between the pectinolytic enzyme sets identified in genomes and the growth on pectins of different origin. BMC Genomics 13, 321 (2012).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
35.Tveit, A. T., Urich, T. & Svenning, M. M. Metatranscriptomic analysis of arctic peat soil microbiota. Appl. Environ. Microbiol. 80, 5761–5772 (2014).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
36.Huergo, L. F. & Dixon, R. The emergence of 2-Oxoglutarate as a master regulator metabolite. Microbiol. Mol. Biol. Rev. 79, 419–435 (2015).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
37.Bai, X. et al. Expression of a β-mannosidase from Paenibacillus polymyxa A-8 in Escherichia coli and characterization of the recombinant enzyme. PLoS ONE 9, e111622 (2014).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
38.Sorensen, J. W., Dunivin, T. K., Tobin, T. C. & Shade, A. Ecological selection for small microbial genomes along a temperate-to-thermal soil gradient. Nat. Microbiol. 4, 55–61 (2019).CAS 
PubMed 
Article 

Google Scholar 
39.Pold, G., Grandy, A. S., Melillo, J. M. & DeAngelis, K. M. Changes in substrate availability drive carbon cycle response to chronic warming. Soil Biol. Biochem. 110, 68–78 (2017).CAS 
Article 

Google Scholar 
40.Yue, H. et al. The microbe-mediated mechanisms affecting topsoil carbon stock in Tibetan grasslands. ISME J. 9, 2012–2020 (2015).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
41.Xue, K. et al. Warming alters expressions of microbial functional genes important to ecosystem functioning. Front. Microbiol. 7, 668–681 (2016).PubMed 
PubMed Central 

Google Scholar 
42.Malik, A. A. et al. Defining trait-based microbial strategies with consequences for soil carbon cycling under climate change. ISME J. 14, 1–9 (2020).CAS 
PubMed 
Article 

Google Scholar 
43.Johnston, E. R. et al. Responses of tundra soil microbial communities to half a decade of experimental warming at two critical depths. Proc. Natl Acad. Sci. USA 116, 15096–15105 (2019).ADS 
CAS 
PubMed 
Article 

Google Scholar 
44.Soong, J. L. et al. Microbial carbon limitation: the need for integrating microorganisms into our understanding of ecosystem carbon cycling. Glob. Change Biol. 26, 1953–1961 (2020).ADS 
Article 

Google Scholar 
45.Chapin, F. S., Matson, P. A. & Vitousek, P. Principles of Terrestrial Ecosystem Ecology (Springer Science & Business Media, 2011).46.Woodcroft, B. J. et al. Genome-centric view of carbon processing in thawing permafrost. Nature 560, 49–54 (2018).ADS 
CAS 
PubMed 
Article 

Google Scholar 
47.Oliverio, A. M., Bradford, M. A. & Fierer, N. Identifying the microbial taxa that consistently respond to soil warming across time and space. Glob. Change Biol. 23, 2117–2129 (2017).ADS 
Article 

Google Scholar 
48.Pold, G. et al. Long-term warming alters carbohydrate degradation potential in temperate forest soils. Appl. Environ. Microbiol. 82, 6518–6530 (2016).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
49.Morrison, E. W. et al. Warming alters fungal communities and litter chemistry with implications for soil carbon stocks. Soil Biol. Biochem. 132, 120–130 (2019).ADS 
CAS 
Article 

Google Scholar 
50.Bai, W., Wang, G., Xi, J., Liu, Y. & Yin, P. Short-term responses of ecosystem respiration to warming and nitrogen addition in an alpine swamp meadow. Eur. J. Soil Biol. 92, 16–23 (2019).CAS 
Article 

Google Scholar 
51.Tilman, D. Niche tradeoffs, neutrality, and community structure: a stochastic theory of resource competition, invasion, and community assembly. Proc. Natl Acad. Sci. USA 101, 10854–10861 (2004).ADS 
CAS 
PubMed 
Article 

Google Scholar 
52.Sharrar, A. M. et al. Bacterial secondary metabolite biosynthetic potential in soil varies with phylum, depth, and vegetation type. mBio 11, e00416–e00420 (2020).PubMed 
PubMed Central 
Article 

Google Scholar 
53.Qi, J. et al. Drying-wetting cycles: effect on deep soil carbon. Soil Syst. 2, 3 (2018).CAS 
Article 

Google Scholar 
54.Butcher, K. R., Nasto, M. K., Norton, J. M. & Stark, J. M. Physical mechanisms for soil moisture effects on microbial carbon-use efficiency in a sandy loam soil in the western United States. Soil Biol. Biochem. 150, 107969 (2020).CAS 
Article 

Google Scholar 
55.Zhou, W. P., Shen, W. J., Li, Y. E. & Hui, D. F. Interactive effects of temperature and moisture on composition of the soil microbial community. Eur. J. Soil Biol. 68, 909–918 (2017).CAS 

Google Scholar 
56.Spohn, M., Klaus, K., Wanek, W. & Richter, A. Microbial carbon use efficiency and biomass turnover times depending on soil depth—Implications for carbon cycling. Soil Biol. Biochem. 96, 74–81 (2016).CAS 
Article 

Google Scholar 
57.Saifuddin, M., Bhatnagar, J. M., Segrè, D. & Finzi, A. C. Microbial carbon use efficiency predicted from genome-scale metabolic models. Nat. Commun. 10, 1–10 (2019).CAS 
Article 

Google Scholar 
58.Geyer, K. M., Kyker-Snowman, E., Grandy, A. S. & Frey, S. D. Microbial carbon use efficiency: accounting for population, community, and ecosystem-scale controls over the fate of metabolized organic matter. Biogeochemistry 127, 173–188 (2016).CAS 
Article 

Google Scholar 
59.Dove, N. C., Safford, H. D., Bohlman, G. N., Estes, B. L. & Hart, S. C. High-severity wildfire leads to multi-decadal impacts on soil biogeochemistry in mixed-conifer forests. Ecol. Appl. 30, e02072 (2020).PubMed 
Article 

Google Scholar 
60.Bradford, M. A. et al. Thermal adaptation of soil microbial respiration to elevated temperature. Ecol. Lett. 11, 1316–1327 (2008).PubMed 
Article 

Google Scholar 
61.Walters, W. et al. Improved bacterial 16S rRNA gene (V4 and V4-5) and fungal internal transcribed spacer marker gene primers for microbial community surveys. mSystems 1, e00009–e00015 (2016).PubMed 
Article 

Google Scholar 
62.Caporaso, J. G. et al. QIIME allows analysis of high-throughput community sequencing data. Nat. Methods 7, 335–336 (2010).CAS 
PubMed 
PubMed Central 
Article 

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

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

Google Scholar 
65.Zhang, Z., Schwartz, S., Wagner, L. & Miller, W. A greedy algorithm for aligning DNA sequences. J. Comput. Biol. 7, 203–214 (2000).CAS 
PubMed 
Article 

Google Scholar 
66.Nguyen, N. H. et al. FUNGuild: An open annotation tool for parsing fungal community datasets by ecological guild. Fungal Ecol. 20, 241–248 (2016).Article 

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

Google Scholar 
68.Menzel, P., Ng, K. L. & Krogh, A. Fast and sensitive taxonomic classification for metagenomics with Kaiju. Nat. Commun. 7, 11257 (2016).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
69.Hyatt, D. et al. Prodigal: prokaryotic gene recognition and translation initiation site identification. BMC Bioinformatics 11, 119 (2010).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
70.Buchfink, B., Xie, C. & Huson, D. H. Fast and sensitive protein alignment using DIAMOND. Nat. Methods 12, 59–60 (2015).CAS 
PubMed 
Article 

Google Scholar 
71.Eddy, S. R. Accelerated profile HMM searches. PLoS Comput. Biol. 7, e1002195 (2011).ADS 
MathSciNet 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
72.Li, D. et al. MEGAHIT v1.0: A fast and scalable metagenome assembler driven by advanced methodologies and community practices. Methods 102, 3–11 (2016).CAS 
PubMed 
Article 

Google Scholar 
73.Bushnell, B. BBMap: A Fast, Accurate, Splice-Aware Aligner. No. LBNL-7065E (Ernest Orlando Lawrence Berkeley National Laboratory, 2014).74.Nurk, S., Meleshko, D., Korobeynikov, A. & Pevzner, P. A. metaSPAdes: a new versatile metagenomic assembler. Genome Res. 27, 824–834 (2017).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
75.Xue, Y., Jonassen, I., Øvreås, L. & Taş, N. Metagenome-assembled genome distribution and key functionality highlight importance of aerobic metabolism in Svalbard permafrost. FEMS Microbiol Ecol. 96, fiaa057 (2020).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
76.Wu, Y.-W., Simmons, B. A. & Singer, S. W. MaxBin 2.0: an automated binning algorithm to recover genomes from multiple metagenomic datasets. Bioinformatics 32, 605–607 (2016).CAS 
PubMed 
Article 

Google Scholar 
77.Kang, D. D., Froula, J., Egan, R. & Wang, Z. MetaBAT, an efficient tool for accurately reconstructing single genomes from complex microbial communities. PeerJ 3, e1165 (2015).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
78.Sieber, C. M. K. et al. Recovery of genomes from metagenomes via a dereplication, aggregation and scoring strategy. Nat. Microbiol. 3, 836–843 (2018).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
79.Parks, D. H., Imelfort, M., Skennerton, C. T., Hugenholtz, P. & Tyson, G. W. CheckM: assessing the quality of microbial genomes recovered from isolates, single cells, and metagenomes. Genome Res. 25, 1043–1055 (2015).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
80.Jain, C., Rodriguez-R, L. M., Phillippy, A. M., Konstantinidis, K. T. & Aluru, S. High throughput ANI analysis of 90K prokaryotic genomes reveals clear species boundaries. Nat. Commun. 9, 5114 (2018).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
81.Bowers, R. M. et al. Minimum information about a single amplified genome (MISAG) and a metagenome-assembled genome (MIMAG) of bacteria and archaea. Nat. Biotechnol. 35, 725–731 (2017).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
82.Chaumeil, P.-A., Mussig, A. J., Hugenholtz, P. & Parks, D. H. GTDB-Tk: a toolkit to classify genomes with the Genome Taxonomy Database. Bioinformatics 36, 1925–1927 (2019).
Google Scholar 
83.Aziz, R. K. et al. The RAST server: rapid annotations using subsystems technology. BMC Genomics 9, 1–15 (2008).Article 
CAS 

Google Scholar 
84.Arkin, A. P. et al. KBase: The United States Department of Energy Systems Biology Knowledgebase. Nat. Biotechnol. 36, 566–569 (2018).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
85.Zhou, Z., Tran, P., Liu, Y., Kieft, K. & Anantharaman, K. METABOLIC: a scalable high-throughput metabolic and biogeochemical functional trait profiler based on microbial genomes. bioRxiv https://doi.org/10.1101/761643 (2019).86.Emiola, A. & Oh, J. High throughput in situ metagenomic measurement of bacterial replication at ultra-low sequencing coverage. Nat. Commun. 9, 1–8 (2018).ADS 
CAS 
Article 

Google Scholar 
87.Binkley, D. & Hart, S. C. in Advances in Soil Science (ed. Stewart, B. A.) 57–112 (Springer New York, 1989).88.Fox, J. & Weisberg, S. An R companion to applied regression Ch.4. (Sage, 2019).89.Bates, D., Maechler, M., Bolker, B. & Walker, S. Fitting linear mixed-effects models using lme4. J. Stat. Softw. 67, 1–48 (2015).Article 

Google Scholar 
90.McMurdie, P. J. & Holmes, S. phyloseq: An R package for reproducible interactive analysis and graphics of microbiome census data. PLoS ONE 8, e61217 (2013).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
91.Love, M. I., Huber, W. & Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15, 550 (2014).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
92.Oksanen, J. et al. vegan: Community Ecology Package. https://cran.r-project.org/web/packages/vegan/index.html (2013).93.Lozupone, C. & Knight, R. UniFrac: a new phylogenetic method for comparing microbial communities. Appl. Environ. Microbiol. 71, 8228–8235 (2005).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
94.Bruns, T. D., White, T. J. & Taylor, J. W. Fungal molecular systematics. Annu. Rev. Ecol. Evol. Syst. 22, 525–564 (1991).Article 

Google Scholar  More

Solving the problem of engineering water shortage is key to ensure water resource security in the karst area It can be seen from the subsystems of the indices sorted by the absolute MIV that the engineering water shortage subsystem had the greatest impact on water resource security in the karst area, which is the main reason to promote its transformation.The water shortage in karst areas is caused by poor natural conditions and inadequate engineering conditions, that is, “engineering water shortage”. It is a serious problem in the Guizhou karst area. The main reasons are as follows. First, the karst hydrogeological and geomorphic conditions, with high mountains and deep rivers, make Guizhou a water shortage area. Second, the karst area is rich in water resources, but it is difficult to develop and utilize these resources. Inter annual variations of rainfall are not significant, but there are large differences within the year, which can easily lead to seasonal drought. Third, the layout of water conservancy projects such as water retention, water storage, and water transfer is unreasonable or insufficient, resulting in conditions of shortage of irrigation and the inadequacy of drinking water for people and livestock. Therefore, the Guizhou karst area has become an area of water shortage, especially engineering water shortage. This is the main bottleneck restricting the coordinated development of the region’s social economy and ecology.Water conservancy projects can determine the diversion and allocation of water resources across time and district to achieve reasonable allocation, efficient utilization, and protection. This indicates the need for higher requirements for engineering water storage and improving water resource utilization efficiency. Therefore, the construction of water conservancy projects is key to ensure future water resource security.The modes of development and utilization of water resources are also significant in the karst area In the past 15 years, Guizhou Province has attached great importance to the development and utilization of water resources. The subsystems of water resource carrying capacity and vulnerability in the Guizhou karst area have risen steadily, which has improved water resource security. However, the development and utilization of water resources will cause changes in the quantity and structure of water usage. This has both optimization and constraints on regional development. Therefore, the geological, hydrological, and hydrogeological characteristics of the karst area must be investigated. The development and utilization of water resources in the karst area should involve appropriate technologies or methods in accordance with these different hydrogeological structures. Geology, geomorphology, rainwater, distributions of farmland and residences, and hydrogeological structures in the karst area are the major factors to consider for solving water shortages in this area35. Rain collection, underground reservoirs, a decentralized water supply and runoff gathering are significant modes of development in the karst area.The situation of water resource security in karst area of Guizhou is gradually getting better This is achieved through water conservation projects and technological measures for water resource exploitation, utilization, projection, and reasonable allocation and control. Meanwhile, Guizhou achieves the security of regional water resource utilization and development through adjusting the regional economic pattern, water resource utilization technology, and so forth.From 2001 to 2006, the status of water resource security was serious, and there was a moderate warning level. At that time, the industrialization of Guizhou province was developing rapidly, and the construction of water conservancy and other infrastructure was also advancing rapidly. Increased attention was given to soil erosion, desertification, water resource pollution, and other problems. Despite high water consumption, the water environment was gradually improving. However, rapid economic and social development has exceeded the carrying capacity of the water resources during this period. Some problems persist in the study area, such as inadequacy of urban sewage treatment facilities, outdated water conservancy facilities, and insufficient prevention of environmental pollution. Urban water pollution treatment facilities and garbage treatment facilities are seriously outdated and cannot meet the requirements of urban development and water environmental protection. These problems have led to a low starting point for water resource security utilization in Guizhou Province. Although the situation has been improved and alleviated year by year, it is still in a moderate warning level, and the water resource security situation is still severe.After reaching the critical safety level in 2007, the water resource security of Guizhou Province declined slightly in 2009 and 2013, although a critical safety level was maintained; the safety level further deteriorated to a moderate warning level in 2011. This deterioration occurred because Guizhou suffered its worst drought in a century from 2009 to 2011, and another drought in 2013. According to the information provided by single indices, the treatment rate of urban waste water, proportion of water supply for water lifting and diversion projects, qualifying rate of water environment function zones, qualifying rate of industrial waste water, degree of development and utilization of groundwater, and density of large and medium-sized reservoirs all showed increasing trends year by year or showed relatively high levels. In contrast, the indices of irrigation water consumption per unit area, above moderate rocky desertification area ratio, water consumption per ten thousand yuan GDP, and water consumption per ten thousand yuan industrial output decreased year by year. All of these indices played a driving role in water utilization and water resource security in the study area. Although the once-in-a-century drought reduced the amount of water, Guizhou Province improved the utilization rate of water resources in the dry years, which alleviated the impact of the reduction of water resources to a certain extent, and allowed the water resource security in the study area to barely maintain the critical safety level. This finding is consistent with previous research conclusions: the engineering water shortage subsystem had largest effect on water resource security in the karst area, whereas the water quantity subsystem had the least influence.It can be inferred that the requirements for ensuring water resource security in the karst area are a good economic development model, environmental protection, pollution control, and improvement of basic water conservancy facilities. These measures can be conducive to actively coping with the impact of abnormal climate changes on the utilization of water resources. More

1.Weiner, J. Physiological limits to sustainable energy budgets in birds and mammals: ecological implications. Trends Ecol. Evol. 7, 384–388 (1992).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
2.Mcnab, B. K. Food habits, energetics, and the population biology of mammals. Am. Nat. 116, 106–124 (1980).Article 

Google Scholar 
3.Hovey, F. W. & Harestad, A. S. Estimating effects of snow on shrub availability for black-tailed deer in southwestern British Columbia. Wildl. Soc. Bull. 20, 308–313 (1992).
Google Scholar 
4.Post, E. & Stenseth, N. Climatic variability, plant phenology, and northern ungulates. Ecology 80, 1322–1339 (1999).Article 

Google Scholar 
5.Moen, A. N. Seasonal changes in heart rates, activity, metabolism, and forage intake of white-tailed deer. J. Wildl. Manag. 42, 715–738 (1978).Article 

Google Scholar 
6.Holand, Ø., Mysterud, A., Wannag, A. & Linnell, J. D. C. Roe deer in northern environments: physiology and behaviour. In The European Roe Deer: Biology of Success (eds Andersen, R. et al.) 117–137 (Scandinavian University Press, 1998).
Google Scholar 
7.Foromozov, A. N. Snow Cover as an Integral Factor of the Environment and Its Importance in the Ecology of Mammals and Birds (The University of Alberta, 1963).
Google Scholar 
8.Cagnacci, F. et al. Partial migration in roe deer: migratory and resident tactics are end points of a behavioural gradient determined by ecological factors. Oikos 120, 1790–1802 (2011).Article 

Google Scholar 
9.Dussault, C., Courtois, R., Ouellet, J.-P. & Girard, I. Space use of moose in relation to food availability. Can. J. Zool. 83, 1431–1437 (2005).Article 

Google Scholar 
10.Mysterud, A. & Sæther, B.-E. Climate change and implications for the future distribution and management of ungulates in Europe. In Ungulate Management in Europe: Problems and Practices (eds Putman, R. et al.) 349–375 (Cambridge University Press, 2011).
Google Scholar 
11.Parmesan, C. Ecological and evolutionary responses to recent climate change. Annu. Rev. Ecol. Evol. Syst. 37, 637–669 (2006).Article 

Google Scholar 
12.Scherrer, S. C., Wüthrich, C., Croci-Maspoli, M., Weingartner, R. & Appenzeller, C. Snow variability in the Swiss Alps 1864–2009. Int. J. Climatol. 33, 3162–3173 (2013).Article 

Google Scholar 
13.Milner, J. M., van Beest, F. M., Schmidt, K. T., Brook, R. K. & Storaas, T. To feed or not to feed? Evidence of the intended and unintended effects of feeding wild ungulates. J. Wildl. Manag. 78, 1322–1334 (2014).Article 

Google Scholar 
14.Ossi, F. et al. Plastic response by a small cervid to supplemental feeding in winter across a wide environmental gradient. Ecosphere 8, e01629 (2017).Article 

Google Scholar 
15.Putman, R. & Staines, B. W. Supplementary winter feeding of wild red deer Cervus elaphus in Europe and North America: justifications, feeding practice and effectiveness. Mamm. Rev. 34, 285–306 (2004).Article 

Google Scholar 
16.Cagnacci, F., Boitani, L., Powell, R. A. & Boyce, M. S. Animal ecology meets GPS-based radiotelemetry: a perfect storm of opportunities and challenges. Philos. Trans. R. Soc. B Biol. Sci. 365, 2157–2162 (2010).Article 

Google Scholar 
17.Peters, W. et al. Migration in geographic and ecological space by a large herbivore. Ecol. Monogr. 87, 297–320 (2017).Article 

Google Scholar 
18.Morellet, N. et al. Seasonality, weather and climate affect home range size in roe deer across a wide latitudinal gradient within Europe. J. Anim. Ecol. 82, 1326–1339 (2013).PubMed 
Article 
PubMed Central 

Google Scholar 
19.Johnson, D. H. The comparison of usage and availability measurements for evaluating resource preference. Ecology 61, 65–71 (1980).Article 

Google Scholar 
20.Ossi, F., Gaillard, J. M., Hebblewhite, M. & Cagnacci, F. Snow sinking depth and forest canopy drive winter resource selection more than supplemental feeding in an alpine population of roe deer. Eur. J. Wildl. Res. 61, 111–124 (2015).Article 

Google Scholar 
21.Mysterud, A. & Østbye, E. Bed-site selection by European roe deer (Capreolus capreolus) in southern Norway during winter. Can. J. Zool. 73, 924–932 (1995).Article 

Google Scholar 
22.Ramanzin, M., Sturaro, E. & Zanon, D. Seasonal migration and home range of roe deer (Capreolus capreolus) in the Italian eastern Alps. Can. J. Zool. 85, 280–289 (2007).Article 

Google Scholar 
23.Endrizzi, S., Gruber, S., Dall’Amico, M. & Rigon, R. GEOtop 2.0: simulating the combined energy and water balance at and below the land surface accounting for soil freezing, snow cover and terrain effects. Geosci. Model. Dev. 7, 2831–2857 (2014).Article 
ADS 

Google Scholar 
24.Cohen, J. A coefficient of agreement for nominal scales. Educ. Psychol. Meas. 20, 37–46 (1960).Article 

Google Scholar 
25.Thomson, A. M. et al. RCP 4.5: a pathway for stabilization of radiative forcing by 2100. Clim. Change 109, 77–94 (2011).CAS 
Article 
ADS 

Google Scholar 
26.Riahi, K. et al. RCP 8.5—a scenario of comparatively high greenhouse gas emissions. Clim. Change 109, 33–57 (2011).CAS 
Article 
ADS 

Google Scholar 
27.Thomas, C. D. Climate, climate change and range boundaries. Divers. Distrib. 16, 488–495 (2010).Article 

Google Scholar 
28.Penteriani, V. et al. Evolutionary and ecological traps for brown bears Ursus arctos in human-modified landscapes. Mamm. Rev. 48, 180–193 (2018).Article 

Google Scholar 
29.Sorensen, A., van Beest, F. M. & Brook, R. K. Impacts of wildlife baiting and supplemental feeding on infectious disease transmission risk: a synthesis of knowledge. Prev. Vet. Med. 113, 356–363 (2014).PubMed 
Article 
PubMed Central 

Google Scholar 
30.Mysterud, A., Viljugrein, H., Solberg, E. J. & Rolandsen, C. M. Legal regulation of supplementary cervid feeding facing chronic wasting disease. J. Wildl. Manag. 83, 1667–1675 (2019).Article 

Google Scholar 
31.Ceacero, F. et al. Benefits for dominant red deer hinds under a competitive feeding system: food access behavior, diet and nutrient selection. PLoS ONE 7, e32780 (2012).CAS 
PubMed 
PubMed Central 
Article 
ADS 

Google Scholar 
32.Beever, E. A. et al. Behavioral flexibility as a mechanism for coping with climate change. Front. Ecol. Environ. 15, 299–308 (2017).Article 

Google Scholar 
33.Loe, L. E. et al. Behavioral buffering of extreme weather events in a high-Arctic herbivore. Ecosphere 7, e01374 (2016).Article 

Google Scholar 
34.Sih, A., Ferrari, M. C. O. & Harris, D. J. Evolution and behavioural responses to human-induced rapid environmental change. Evol. Appl. 4, 367–387 (2011).PubMed 
PubMed Central 
Article 

Google Scholar 
35.Radchuk, V. et al. Adaptive responses of animals to climate change are most likely insufficient. Nat. Commun. 10, 3109 (2019).PubMed 
PubMed Central 
Article 
ADS 
CAS 

Google Scholar 
36.Mysterud, A. Still walking on the wild side? Management actions as steps towards ‘semi-domestication’ of hunted ungulates. J. Appl. Ecol. 47, 920–925 (2010).Article 

Google Scholar 
37.Felton, A. M. et al. Interactions between ungulates, forests, and supplementary feeding: the role of nutritional balancing in determining outcomes. Mamm. Res. 62, 1–7 (2017).Article 

Google Scholar 
38.Ricci, S. et al. Impact of supplemental winter feeding on ruminal microbiota of roe deer Capreolus capreolus. Wildl. Biol. 2019, wlb.00572 (2019).Article 

Google Scholar 
39.Lone, K. et al. Living and dying in a multi-predator landscape of fear: roe deer are squeezed by contrasting pattern of predation risk imposed by lynx and humans. Oikos 123, 641–651 (2014).Article 

Google Scholar 
40.Chapron, G. et al. Recovery of large carnivores in Europe’s modern human-dominated landscapes. Science (80-) 346, 1517–1519 (2014).CAS 
Article 
ADS 

Google Scholar 
41.Milanesi, P., Breiner, F. T., Puopolo, F. & Holderegger, R. European human-dominated landscapes provide ample space for the recolonization of large carnivore populations under future land change scenarios. Ecography (Cop.) 40, 1359–1368 (2017).Article 

Google Scholar 
42.Pascual-Rico, R. et al. Is diversionary feeding a useful tool to avoid human-ungulate conflicts? A case study with the aoudad. Eur. J. Wildl. Res. 64, 1–7 (2018).Article 

Google Scholar 
43.van Beest, F. M., Loe, L. E., Mysterud, A. & Milner, J. M. Comparative space use and habitat selection of moose around feeding stations. J. Wildl. Manag. 74, 219–227 (2010).Article 

Google Scholar 
44.Jerina, K. Roads and supplemental feeding affect home-range size of Slovenian red deer more than natural factors. J. Mamm. 93, 1139–1148 (2012).Article 

Google Scholar 
45.Ranc, N. et al. Preference and familiarity mediate spatial responses of a large herbivore to experimental manipulation of resource availability. Scientific Reports 10, 11946 (2020). 46.Brown, R. D. & Robinson, D. A. Northern Hemisphere spring snow cover variability and change over 1922–2010 including an assessment of uncertainty. Cryosphere 5, 219–229 (2011).Article 
ADS 

Google Scholar 
47.Schloss, C. A., Nuñez, T. A. & Lawler, J. J. Dispersal will limit ability of mammals to track climate change in the Western Hemisphere. Proc. Natl. Acad. Sci. U. S. A. 109, 8606–8611 (2012).CAS 
PubMed 
PubMed Central 
Article 
ADS 

Google Scholar 
48.Gurarie, E. et al. A framework for modelling range shifts and migrations: asking when, whither, whether and will it return. J. Anim. Ecol. 86, 943–959 (2017).PubMed 
Article 

Google Scholar 
49.Rivrud, I. M. et al. Leave before it’s too late: anthropogenic and environmental triggers of autumn migration in a hunted ungulate population. Ecology 97, 1058–1065 (2016).PubMed 
Article 
PubMed Central 

Google Scholar 
50.Courtois, R., Dussault, C., Potvin, F. & Daigle, G. Habitat selection by moose (Alces alces) in clear-cut landscapes. Alces 38, 177–192 (2002).
Google Scholar 
51.Gilbert, S. L., Hundertmark, K. J., Person, D. K., Lindberg, M. S. & Boyce, M. S. Behavioral plasticity in a variable environment: snow depth and habitat interactions drive deer movement in winter. J. Mamm. 98, 246–259 (2017).Article 

Google Scholar 
52.Chevin, L. M., Lande, R. & Mace, G. M. Adaptation, plasticity, and extinction in a changing environment: towards a predictive theory. PLoS Biol. 8, e1000357 (2010).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
53.Bauer, S. & Hoye, B. J. Migratory animals couple biodiversity and ecosystem functioning worldwide. Science (80-) 344, 1242552 (2014).CAS 
Article 

Google Scholar 
54.Mason, T. H. E., Stephens, P. A., Apollonio, M. & Willis, S. G. Predicting potential responses to future climate in an alpine ungulate: Interspecific interactions exceed climate effects. Glob. Change Biol. 20, 3872–3882 (2014).Article 
ADS 

Google Scholar 
55.Carnevali, L., Pedrotti, L., Riga, F. & Toso, S. Banca dati ungulati: Status, distribuzione, consistenza, gestione e prelievo venatorio delle popolazioni di ungulati in Italia. Rapporto 2001–2005 Vol. 117 (Biologia e Conservazione della Fauna, 2009).
Google Scholar 
56.Provincia Autonoma di Trento. Analisi delle consistenze e dei prelievi di ungulati, tetraonidi e coturnice. Stagione Venatoria 2018 (Provincia Autonoma di Trento, 2018).
Google Scholar 
57.Rockel, B., Will, A. & Hense, A. The regional climate model COSMO-CLM (CCLM). Meteorol. Z. 17, 347–348 (2008).Article 

Google Scholar 
58.Boyce, M. S. & McDonald, L. L. Relating populations to habitats using resource selection functions. Trends Ecol. Evol. 14, 268–272 (1999).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
59.Boyce, M. S., Vernier, P. R., Nielsen, S. E. & Schmiegelow, F. K. A. Evaluating resource selection functions. Ecol. Modell. 157, 281–300 (2002).Article 

Google Scholar 
60.Benoit, T. & Achraf, E. suncalc: compute sun position, sunlight phases, moon position and lunar phase. R package version 0.5.0. https://cran.r-project.org/package=suncalc (2019).61.DeCesare, N. J. et al. Transcending scale dependece in identifying habitat with resource selection functions. Ecol. Appl. 22, 1068–1083 (2012).PubMed 
Article 
PubMed Central 

Google Scholar 
62.Kendall, M. A new measure of rank correlation. Biometrika 30, 81–89 (1938).MATH 
Article 

Google Scholar 
63.Cohen, J. Weighted kappa: nominal scale agreement with provision for scaled disagreement or partial credit. Psychol. Bull. 70, 213–220 (1968).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
64.Gamer, M., Lemon, J., Fellows, I. & Singh, P. irr: various coefficients of interrater reliability and agreement. R package version 0.84.1. https://cran.r-project.org/package=irr (2019).65.Lele, S. R., Keim, J. L. & Solymos, P. ResourceSelection: resource selection (probability) functions for use-availability data. R package version 0.3-5. https://cran.r-project.org/package=ResourceSelection (2019).66.Bivand, R., Keitt, T. & Rowlingson, B. rgdal: bindings for the ‘Geospatial’ Data Abstraction Library. R package version 1.4-8. https://cran.r-project.org/package=rgdal (2019).67.McLeod, A. I. Kendall: Kendall rank correlation and Mann-Kendall trend test. R package version 2.2. https://cran.r-project.org/package=Kendall (2011).68.Bright Ross, J. G., Peters, W., Ossi, F., Moorcroft P. R., Cordano, E., Eccel, E., Bianchini, F., Ramanzin, M., and Cagnacci, F. Datasets for “Climate change and anthropogenic food manipulation interact in shifting the distribution of a large herbivore at its altitudinal range limit.” https://doi.org/10.5281/zenodo.4637674 (2021). More

1.Alexandra-Maria, K. et al. Importance of pollinators in changing landscapes for world crops. Proc. R. Soc. B Biol. Sci. 274, 303–313 (2007).Article 

Google Scholar 
2.Winfree, R., Williams, N. M., Gaines, H., Ascher, J. S. & Kremen, C. Wild bee pollinators provide the majority of crop visitation across land-use gradients in New Jersey and Pennsylvania, USA. J. Appl. Ecol. 45, 793–802 (2008).Article 

Google Scholar 
3.Brosi, B. J. & Briggs, H. M. Single pollinator species losses reduce floral fidelity and plant reproductive function. Proc. Natl. Acad. Sci. 110, 13044–13048 (2013).CAS 
Article 
ADS 

Google Scholar 
4.Potts, S. G. et al. Global pollinator declines: trends, impacts and drivers. Trends Ecol. Evol. 25, 345–353 (2010).Article 

Google Scholar 
5.Cameron, S. A. et al. Patterns of widespread decline in North American bumble bees. Proc. Natl. Acad. Sci. 108, 662–667 (2011).CAS 
Article 
ADS 

Google Scholar 
6.Koh, I. et al. Modeling the status, trends, and impacts of wild bee abundance in the United States. Proc. Natl. Acad. Sci. 113, 140–145 (2016).CAS 
Article 
ADS 

Google Scholar 
7.Cameron, S. A. & Sadd, B. M. Global trends in bumble bee health. Annu. Rev. Entomol. 65, 209–232 (2020).CAS 
Article 

Google Scholar 
8.Murray, T. E., Kuhlmann, M. & Potts, S. G. Conservation ecology of bees: populations, species and communities. Apidologie 40, 211–236 (2009).Article 

Google Scholar 
9.Michener, C. D. The Bees of the World (Johns Hopkins University Press, Baltimore, 2007).
Google Scholar 
10.Milam, J. et al. Validating morphometrics with DNA barcoding to reliably separate three cryptic species of bombus cresson (Hymenoptera: Apidae). Insects 11, 669 (2020).Article 

Google Scholar 
11.Williams, P. H. et al. Widespread polytypic species or complexes of local species? Revising bumblebees of the subgenus Melanobombus world-wide (Hymenoptera, Apidae, Bombus). Eur. J. Taxon. 719, 1–120 (2020).
Google Scholar 
12.Drew, L. W. Are we losing the science of taxonomy? As need grows, numbers and training are failing to keep up. Bioscience 61, 942–946 (2011).Article 

Google Scholar 
13.Portman, Z. M., Bruninga-Socolar, B. & Cariveau, D. P. The state of bee monitoring in the United States: A call to refocus away from bowl traps and towards more effective methods. Ann. Entomol. Soc. Am. 113, 337–342 (2020).Article 

Google Scholar 
14.Valan, M., Makonyi, K., Maki, A., Vondráček, D. & Ronquist, F. Automated taxonomic identification of insects with expert-level accuracy using effective feature transfer from convolutional networks. Syst. Biol. 68, 876–895 (2019).Article 

Google Scholar 
15.Gratton, C. & Zuckerberg, B. Citizen science data for mapping bumble bee populations, in Novel Quantitative Methods in Pollinator Ecology and Management (2019).16.MacPhail, V. J., Gibson, S. D., Hatfield, R. & Colla, S. R. Using Bumble Bee Watch to investigate the accuracy and perception of bumble bee (Bombus spp.) identification by community scientists. PeerJ 8, e9412 (2020).Article 

Google Scholar 
17.Weeks, P. J. D., Gauld, I. D., Gaston, K. J. & O’Neill, M. A. Automating the identification of insects: a new solution to an old problem. Bull. Entomol. Res. 87, 203–211 (1997).Article 

Google Scholar 
18.Schröder, S. et al. The new key to bees: Automated identification by image analysis of wings. in The Conservation Link Between Agriculture and Nature (eds. Kevan, P. & Imperatriz-Fonseca, V.) 209–216 (Ministry of Environment, 2002).19.MacLeod, N., Benfield, M. & Culverhouse, P. Time to automate identification. Nature 467, 154–155 (2010).CAS 
Article 
ADS 

Google Scholar 
20.Fuentes, A., Yoon, S., Kim, S. C. & Park, D. S. A robust deep-learning-based detector for real-time tomato plant diseases and pests recognition. Sensors 17, 2022 (2017).Article 

Google Scholar 
21.Motta, D. et al. Application of convolutional neural networks for classification of adult mosquitoes in the field. PLoS ONE 14, e0210829 (2019).CAS 
Article 

Google Scholar 
22.Bojarski, M. et al. End to end learning for self-driving cars. arXxiv:1604.07316 (2016).23.Anthimopoulos, M., Christodoulidis, S., Ebner, L., Christe, A. & Mougiakakou, S. Lung pattern classification for interstitial lung diseases using a deep convolutional neural network. IEEE Trans. Med. Imaging 35, 1207–1216 (2016).Article 

Google Scholar 
24.Liu, Z., Gao, J., Yang, G., Zhang, H. & He, Y. Localization and classification of paddy field pests using a saliency map and deep convolutional neural network. Sci. Rep. 6, 20410 (2016).CAS 
Article 
ADS 

Google Scholar 
25.Martineau, M., Raveaux, R., Chatelain, C., Conte, D. & Venturini, G. Effective training of convolutional neural networks for insect image recognition. In Advanced Concepts for Intelligent Vision Systems, pp 426–437 (eds Blanc-Talon, J. et al.) (Springer International Publishing, Cham, 2018).
Google Scholar 
26.Marques, A. C. R. et al. Ant genera identification using an ensemble of convolutional neural networks. PLoS ONE 13, e0192011 (2018).Article 

Google Scholar 
27.Williams, P. H., Thorp, R. W., Richardson, L. L. & Colla, S. R. Bumble Bees of North America: An Identification Guide (Princeton University Press, Princeton, 2014).
Google Scholar 
28.He, K., Zhang, X., Ren, S. & Sun, J. Deep residual learning for image recognition. 2016 IEEE Conf. Comput. Vis. Pattern Recognit. CVPR 770–778 (2015).29.Zagoruyko, S. & Komodakis, N. Wide residual networks. arXxiv:1605.07146 (2017).30.Szegedy, C., Vanhoucke, V., Ioffe, S., Shlens, J. & Wojna, Z. Rethinking the inception architecture for computer vision. arXxiv:1512.00567 (2015).31.Tan, M. et al. MnasNet: Platform-aware neural architecture search for mobile. arXxiv:1807.11626 (2019).32.Deng, J. et al. ImageNet: A large-scale hierarchical image database. in 2009 IEEE Conference on Computer Vision and Pattern Recognition 248–255 (2009).33.Hernández-García, A. & König, P. Further advantages of data augmentation on convolutional neural networks. arXxiv:1906.11052 11139, 95–103 (2018).34.Fard, F. S., Hollensen, P., Mcilory, S. & Trappenberg, T. Impact of biased mislabeling on learning with deep networks. in 2017 International Joint Conference on Neural Networks (IJCNN) 2652–2657 (2017).35.Clare, J. D. J., Townsend, P. A. & Zuckerberg, B. Generalized model-based solutions to false positive error in species detection/non-detection data. Ecology 102, e03241 (2021).Article 

Google Scholar 
36.Clare, J. D. J. et al. Making inference with messy (citizen science) data: when are data accurate enough and how can they be improved?. Ecol. Appl. 29, e01849 (2019).Article 

Google Scholar 
37.Tian, Z. et al. Discriminative CNN via metric learning for hyperspectral classification. in IGARSS 2019 – 2019 IEEE International Geoscience and Remote Sensing Symposium 580–583 (2019).38.Nazki, H., Yoon, S., Fuentes, A. & Park, D. S. Unsupervised image translation using adversarial networks for improved plant disease recognition. Comput. Electron. Agric. 168, 105117 (2020).Article 

Google Scholar 
39.Wäldchen, J. & Mäder, P. Machine learning for image based species identification. Methods Ecol. Evol. 9, 2216–2225 (2018).Article 

Google Scholar 
40.Woodard, S. H. et al. Towards a U.S. national program for monitoring native bees. Biol. Conserv. 252, 108821 (2020).Article 

Google Scholar 
41.Wagner, D. L. Insect declines in the anthropocene. Annu. Rev. Entomol. 65, 457–480 (2020).CAS 
Article 

Google Scholar 
42.Montgomery, G. A. et al. Is the insect apocalypse upon us? How to find out. Biol. Conserv. 241, 108327 (2020).Article 

Google Scholar 
43.Høye, T. T., Mann, H. M. R. & Bjerge, K. Camera-based monitoring of insects on green roofs. DCE – Natl. Cent. Environ. Energy 18 (2020).44.Ärje, J. et al. Automatic image-based identification and biomass estimation of invertebrates. Methods Ecol. Evol. 11, 922–931 (2020).Article 

Google Scholar 
45.Norouzzadeh, M. S. et al. Automatically identifying, counting, and describing wild animals in camera-trap images with deep learning. Proc. Natl. Acad. Sci. 115, E5716–E5725 (2018).CAS 
Article 

Google Scholar 
46.Hallmann, C. A. et al. More than 75 percent decline over 27 years in total flying insect biomass in protected areas. PLoS ONE 12, e0185809 (2017).Article 

Google Scholar 
47.Ghisbain, G. et al. Substantial genetic divergence and lack of recent gene flow support cryptic speciation in a colour polymorphic bumble bee (Bombus bifarius) species complex. Syst. Ecol. 45, 635–652 (2020).
Google Scholar  More

Thank you for visiting nature.com. You are using a browser version with limited support for CSS. To obtain
the best experience, we recommend you use a more up to date browser (or turn off compatibility mode in
Internet Explorer). In the meantime, to ensure continued support, we are displaying the site without styles
and JavaScript. More