Philippa Kaur

Bethany, J., Giraldo-Silva, A., Nelson, C., Barger, N. N. & Garcia-Pichel, F. Optimizing the production of nursery-based biological soil crusts for restoration of arid land soils. Appl. Environ. Microbiol. 85(15), e00735-e819 (2019).CAS 
PubMed 
PubMed Central 

Google Scholar 
Belnap, J. & Gardner, J. S. Soil microstructure in soils of the colorado plateau—The role of the cyanobacterium Microcoleus-vaginatus. Gt. Basin Nat. 53(1), 40–47 (1993).
Google Scholar 
Belnap, J. Factors influencing nitrogen fixation and nitrogen release in biological soil crusts. Ecol. Stud. Biol. Soil Crusts Struct. Funct. Manag. 150, 241–261 (2001).
Google Scholar 
Cameron, R. E. & Blank, G. B. Desert algae: Soil crusts and diaphanous substrata as algal habitats. Tech. Rep. Jet Propul. Lab. Calif. Technol. 32–971, 1–41 (1966).
Google Scholar 
Friedmann EI, Galun M. Desert algae lichens and fungi. in Desert Biology (Brown Jr, G.W. eds). Vol. 2. 165–212. (Illus Academic Press, Inc., 1974).Maier, S., Tamm, A., Wu, D.A.-O., Caesar, J., Grube, M., & Weber, B.A.-O. Photoautotrophic Organisms Control Microbial Abundance, Diversity, and Physiology in Different Types of Biological Soil Crusts. (1751–7370 (electronic)).Cable, J. M. & Huxman, T. E. Precipitation pulse size effects on Sonoran Desert soil microbial crusts. Oecologia 141(2), 317–324 (2004).ADS 
PubMed 

Google Scholar 
Evans, R. D. & Johansen, J. R. Microbiotic crusts and ecosystem processes. Crit. Rev. Plant Sci. 18(2), 183–225 (1999).
Google Scholar 
Thompson, J. N. et al. Frontiers of ecology. Bioscience 51(1), 15–24 (2001).
Google Scholar 
Warren, S. D., Rosentreter, R. & Pietrasiak, N. Biological soil crusts of the Great Plains: A review. Rangel Ecol. Manag. 1(78), 213–219 (2021).
Google Scholar 
Warren, S. D. et al. Biological soil crust response to late season prescribed fire in a Great Basin Juniper Woodland. Rangel. Ecol. Manag. 68(3), 241–247 (2015).
Google Scholar 
Thomas, A. D., Hoon, S. R. & Linton, P. E. Carbon dioxide fluxes from cyanobacteria crusted soils in the Kalahari. Appl. Soil Ecol. 39, 254–263 (2008).
Google Scholar 
Williams, A. J., Buck, B. J. & Beyene, M. A. Biological soil crusts in the Mojave Desert, USA: Micromorphology and pedogenesis. Soil Sci. Soc. Am. J. 76(5), 1685–1695 (2012).ADS 
CAS 

Google Scholar 
Belnap, J. & Lange, O. L. Ecological studies: Biological soil crusts: Structure, function, and management. Ecol. Stud. Biol. Soil Crusts Struct. Funct. Manag. 150, 1–503 (2001).
Google Scholar 
Jordan, W. R. I. Restoration ecology: A synthetic approach to ecological research. Rehabil. Damaged Ecosyst. 2, 373–384 (1995).
Google Scholar 
Nash, T. H. et al. Photosynthetic patterns of Sonoran desert lichens.1. Environmental considerations and preliminary field-measurements. Flora 172(4), 335–345 (1982).
Google Scholar 
St. Clair, L. L., Johansen, J. R. & Rushforth, S. R. Lichens of soil crust communities in the Intermountain Area of the western United States. Gt Basin Nat. 53(1), 5 (1993).
Google Scholar 
Bowker, M. A., Belnap, J. & Miller, M. E. Spatial modeling of biological soil crusts to support rangeland assessment and monitoring. Rangel. Ecol. Manag. 59(5), 519–529 (2006).
Google Scholar 
Mayland, H. F., McIntosh, T. H. & Fuller, W. H. Fixation of isotopic nitrogen on a semiarid soil by algal crust organisms. Soil Sci. Soc. Am. Proc. 30(1), 56 (1966).ADS 
CAS 

Google Scholar 
McIlvanie, S. K. Grass seedling establishment, and productivity—Overgrazed vs. protected range soils. Ecology 23(2), 228–231 (1942).
Google Scholar 
Webb, R. H. & Wilshire, H. G. Environmental Effects of Off-Road Vehicles : Impacts and Management in Arid Regions (Springer, 1983).
Google Scholar 
Zobel, D. & Antos, J. A decade of recovery of understory vegetation buried by volcanic tephra from Mount St. Helens. Ecol. Monogr. 1, 67 (1997).
Google Scholar 
Condon, L. & Pyke, D. Resiliency of biological soil crusts and vascular plants varies among morphogroups with disturbance intensity. Plant Soil. 12, 433 (2020).
Google Scholar 
Harper, K., & Marble, J. A role for nonvascular plants in management of arid and semiarid rangelands. in Vegetation Science Applications for Rangeland Analysis and Management [Internet] (Tueller, P.T., ed.). Handbook of Vegetation Science. Vol. 14. 135–169. https://doi.org/10.1007/978-94-009-3085-8_7. (Springer, 1988). Evans, R. D. & Belnap, J. Long-term consequences of disturbance on nitrogen dynamics in an arid ecosystem. Ecology 80(1), 150–160 (1999).
Google Scholar 
Sheridan, R. P. Impact of emissions from coal-fired electricity generating facilities on N2-fixing lichens. Bryologist 82(1), 54–58 (1979).CAS 

Google Scholar 
Henriksson, L. E. & Dasilva, E. J. Effects of some inorganic elements on nitrogen-fixation in blue-green-algae and some ecological aspects of pollution. Z. Allg. Mikrobiol. 18(7), 487–494 (1978).CAS 
PubMed 

Google Scholar 
Freebury, C. Lichens and lichenicolous fungi of Grasslands National Park (Saskatchewan, Canada). Opusc Philolichenum 13, 102–121 (2009).
Google Scholar 
Szyja, M. et al. Neglected but potent dry forest players: ecological role and ecosystem service provision of biological soil crusts in the human-modified Caatinga. Front. Ecol. Evol. (Internet). https://doi.org/10.3389/fevo.2019.00482 (2019).Article 

Google Scholar 
Rosentreter, R. Biological soil of crusts of North American drylands: Cryptic diversity at risk. in Reference Module in Earth Systems and Environmental Sciences [Internet]. https://www.sciencedirect.com/science/article/pii/B9780128211397000738 (Elsevier, 2021). Kranz, C. N., McLaughlin, R. A., Johnson, A., Miller, G. & Heitman, J. L. The effects of compost incorporation on soil physical properties in urban soils—A concise review. J. Environ. Manag. 261, 110209 (2020).
Google Scholar 
Barberán, A. et al. Continental-scale distributions of dust-associated bacteria and fungi. Proc. Natl. Acad. Sci. 112(18), 5756 (2015).ADS 
PubMed 
PubMed Central 

Google Scholar 
Kaye, J. P., Groffman, P. M., Grimm, N. B., Baker, L. A. & Pouyat, R. V. A distinct urban biogeochemistry?. Trends Ecol. Evol. 21(4), 192–199 (2006).PubMed 

Google Scholar 
Pavao-Zuckerman, M. A. The nature of urban soils and their role in ecological restoration in cities. Restor. Ecol. 16(4), 642–649 (2008).
Google Scholar 
Pouyat, R., Groffman, P., Yesilonis, I. & Hernandez, L. Soil carbon pools and fluxes in urban ecosystems. Environ. Pollut. 116, S107–S118 (2002).CAS 
PubMed 

Google Scholar 
Behzad, H., Mineta, K., & Gojobori, T. Global Ramifications of Dust and Sandstorm Microbiota. (1759–6653 (electronic)).Warren, S., Clair, L. & Leavitt, S. Aerobiology and passive restoration of biological soil crusts. Aerobiologia 3, 35 (2021).
Google Scholar 
Hall, S. J. et al. Urbanization alters soil microbial functioning in the Sonoran Desert. Ecosystems 12(4), 654–671 (2009).CAS 

Google Scholar 
Ball, B. A. & Guevara, J. A. The nutrient plasticity of moss-dominated crust in the urbanized Sonoran Desert. Plant Soil. 389(1–2), 225–235 (2015).CAS 

Google Scholar 
Allen, C. D. Monitoring environmental impact in the Upper Sonoran lifestyle: A new tool for rapid ecological assessment. Environ. Manag. 43(2), 346–356 (2009).ADS 

Google Scholar 
Evans, R. A. & Love, R. M. The step-point method of sampling: A practical tool in range research. J. Range Manag. 10(5), 208–212 (1957).
Google Scholar 
Coulloudon, B., & National Applied Resource Sciences C. Sampling Vegetation Attributes Interagency Technical Reference [Internet]. http://www.blm.gov/nstc/library/pdf/samplveg.pdf. (Bureau of Land Management : National Business Center, 1999). Faithfull, N. T. Methods in agricultural chemical analysis: A practical handbook. Methods Agric. Chem. Anal. Pract. Handb. 1–22, 1–266 (2002).
Google Scholar 
Kuske, C. R., Yeager, C. M., Johnson, S., Ticknor, L. O. & Belnap, J. Response and resilience of soil biocrust bacterial communities to chronic physical disturbance in arid shrublands. ISME J. 6(4), 886–897 (2012).CAS 
PubMed 

Google Scholar 
Lorenz, K. & Lal, R. Biogeochemical C and N cycles in urban soils. Environ. Int. 35(1), 1–8 (2009).CAS 
PubMed 

Google Scholar 
Chamizo, S., Canton, Y., Lazaro, R., Sole-Benet, A. & Domingo, F. Crust composition and disturbance drive infiltration through biological soil crusts in semiarid ecosystems. Ecosystems 15(1), 148–161 (2012).
Google Scholar 
Kidron, G. J. & Gutschick, V. P. Soil moisture correlates with shrub-grass association in the Chihuahuan Desert. CATENA 107, 71–79 (2013).
Google Scholar 
Kidron, G. J., Monger, H. C., Vonshak, A. & Conrod, W. Contrasting effects of microbiotic crusts on runoff in desert surfaces. Geomorphology 15(139), 484–494 (2012).ADS 

Google Scholar 
Berdugo, M., Soliveres, S. & Maestre, F. T. Vascular plants and biocrusts modulate how abiotic factors affect wetting and drying events in drylands. Ecosystems 17, 1242 (2014).CAS 

Google Scholar 
Maestre, F. T. et al. Changes in biocrust cover drive carbon cycle responses to climate change in drylands. Glob Change Biol. 19, 3835 (2013).ADS 

Google Scholar 
Valenzuela, A. et al. Aerosol radiative forcing during African desert dust events (2005–2010) over southeastern Spain. Atmos. Chem. Phys. 12(21), 10331–10351 (2012).ADS 
CAS 

Google Scholar 
Kaya, S., Basar, U. G., Karaca, M. & Seker, D. Z. Assessment of urban heat islands using remotely sensed data. Ekoloji 21(84), 107–113 (2012).
Google Scholar 
Demmigadams, B. et al. Effect of high light on the efficiency of photochemical energy-conversion in a variety of lichen species with green and blue-green phycobionts. Planta 180(3), 400–409 (1990).CAS 

Google Scholar 
Gauslaa, Y. & Rikkinen, J. What’s behind the pretty colours? A study on the photobiology of lichens. Nord. J. Bot. 17(5), 556–556 (1995).
Google Scholar 
Garciapichel, F. & Castenholz, R. W. Characterization and biological implications of scytonemin, a cyanobacterial sheath pigment. J. Phycol. 27(3), 395–409 (1991).CAS 

Google Scholar 
Garcia-Pichel, F. & Castenholz, R. W. The role of sheath pigments in the adaptation of terrestrial cyanobacteria to near UV radiation. J. Phycol. 27(3 SUPPL), 24–24 (1991).
Google Scholar 
McDonnell, M. J. et al. Ecosystem processes along an urban-to-rural gradient. Urban Ecosyst. 1(1), 21–36 (1997).
Google Scholar 
Pavao-Zuckerman, M. A. & Byrne, L. B. Scratching the surface and digging deeper: Exploring ecological theories in urban soils. Urban Ecosyst. 12(1), 9–20 (2009).
Google Scholar 
Pavao-Zuckerman, M. A. Urban greenscape, soils, and ecosystem functioning in a semi-arid urban ecosystem. J. Nematol. 41(4), 369–370 (2009).
Google Scholar 
Collins, S. L. et al. Pulse dynamics and microbial processes in aridland ecosystems. J. Ecol. 96(3), 413–420 (2008).
Google Scholar 
Noy-Meir, I. Desert ecosystems environment and producers. In Annual Review on Ecology System (Johnston Richard, F. ed.). Vol. 4. 25–51. (Illus Map Annu Rev Inc, 1973). More

Jean Combes’s love of nature as a child led her to note the signs of starting spring. Her long-term records are now part of a vital growing citizen science dataset that starkly shows how climate change is shifting the timing of the natural world.For people living in colder parts of the world, watching for the first signs of spring — from the opening of snowdrops and daffodils, to birds building their nests, to the return of bees and butterflies — is a common winter pastime. Jean Combes has not just been watching out for these signs, but also recording them, ever since she was a child. Taking note of the earliest emergence of leaves in springtime — first as a child of 11 years, and then continuously from the age of 20 years — Jean has now collected one of the longest continuous datasets of spring leaf-out time in the UK (see also Correspondence by Vitasse et al.). These almost 75 years of data show a clear shift that corroborates shifts now acknowledged for diverse species around the world: springtime is coming earlier, and the patterns of advance match the global trends in the changing climate. Jean’s naturalist endeavours have already earned her high honours in the form of an OBE (Order of the British Empire), and recognition of her own work is mirrored in a growing recognition of the vital role of citizen scientists in tracking the signs of our rapidly changing world.
This is a preview of subscription content More

System dynamics methodThe SD method applies systemic processing to simulate complex non-linear dynamics and feedback. Systemic processing resorts to various tools to simulate complex system behavior and performance24. Systems evolve through states, which change with flows. An example of a state variable is water storage in the study of lakes. The SD method simulates changes in system states driven by flows and various feedbacks25.This work employs the SD method to simulate storage change in Lake Urmia in one historical period (1957–2005) and two future periods (2021–2050 and 2051–2080). The lake’s water volume is the state variable, which is governed by inflows (precipitation, surface water inflows, and groundwater inflows) and outflows (evaporation, leakage, and surface water outflows). The lake’s mass balance equation is expressed as:$$S_{t + 1} = intlimits_{t}^{t + 1} {[I_{s} – O_{s} ]ds + S_{t} }$$
(1)
where St+1 , St, Is, and Os denote the lake’s storage at time t + 1, the lake’s storage at time t, the inflow rate to the lake at time s (units of volume/time), and the outflow rate from the lake at time s (units of volume/time), respectively.The SD method employs the Euler and Runge Kutta methods for the solution of differential equations. The software STELLA, Vensim, Powersim, and Dynamo feature SD solvers26. This work applies the widely-used Vensim software27.Climate changeThe data sets needed for modeling Lake Urmia’s storage over the two future periods were generated after simulating the lake’s water balance during the historical period. HADCM3, a coupled atmosphere–ocean general circulation model’s (AOGCM) climate projections were used to generate precipitation and surface temperature projections over the future periods. The AOGCM data at coarse spatial scales were downscaled to the regional scale suitable for lake storage simulation. The commonly used downscaling methods are statistic and dynamic in nature28,29. This works applies the delta-change downscaling method, in which monthly temperature and precipitation differences between the future and historical are calculated by29:$$Delta T_{t} = overline{T}_{GCM,fut,t} – overline{T}_{GCM,hist,t}$$
(2)
$$Delta P_{t} = overline{P}_{GCM,fut,t} – overline{P}_{GCM,hist,t}$$
(3)
where ∆Tt denotes the difference in long-term average temperatures simulated by HADCM3 for the future ((overline{T}_{GCM,fut,t})) and historical ((overline{T}_{GCM,hist,t})) periods in month t (°C); ∆Pt represents the difference in long-term average precipitations simulated by HADCM3 for the future ((overline{P}_{GCM,fut,t})) and historical ((overline{P}_{GCM,hist,t})) periods in month t (mm). Then, ∆Tt and ∆Pt are applied to project the future downscaled data as follows29:$$T_{t} = T_{obs,t} + , Delta T_{t}$$
(4)
$$P_{t} = P_{obs,t} { + }Delta P_{t}$$
(5)
where Tobs,t, and Pobs,t denote respectively the observed temperature (°C) and precipitation (mm) in month t in the baseline period; and Tt and Pt are the downscaled temperature (°C) and precipitation (mm) in month t of the future period, respectively. Delta-change downscaling is a simple yet efficient option when it comes to spatial downscaling of climate change projections (e.g.30,31,32). The gist of this method is to replicate the changing patterns that are projected by the atmospheric ocean general circulation models (AOGCMs) to generate the climate change patterns of hydro-climatic variables on a regional scale. As such, one would simply compute the relative changes in the long-term variations of the variable that is projected by the models within the baseline and future timeframes. These relative changing patterns would be applied to the historical data to project the impact of climate change on a local scale.Rainfall-runoff modelingThe IHACRES (identification of unit hydrographs and component flows from rainfall, evapotranspiration and streamflow) model is herein applied to simulate runoff from precipitation. Ashofteh et al.33 implemented the IHACRES model to investigate the effects of climate change on reservoir performance in agricultural water supply. Ashofteh et al.34 evaluated the probability of flood occurrence in future periods with IHACRES.The IHACRES model includes a non-linear loss module and a linear unit hydrograph module. The non-linear loss module converts the observed rainfall into the effective rainfall, after which the linear unit hydrograph module converts the effective rainfall into the simulated streamflow35. Here, precipitation rk in time step k is converted to effective precipitation uk through the non-linear loss module employing a catchment wetness index sk:$$u_{k} = , s_{k} times , r_{k}$$
(6)
The effective precipitation is converted to the surface runoff in time step k with the linear unit hydrograph module. The parameters of this model can be set through a thorough grid numeric search and trial-and-error. Perhaps, one of the major advantages of the IHACRES model over other commonly-used rainfall-runoff models is its minimal input data requirement (i.e., air temperature and precipitation)31,35.The other alternative for hydrologic simulation is to use data-driven models. Here, the multilayer perceptron (MLP), a variety of the artificial neural network (ANN) method, was also used to simulate runoff. This model consists of an inlet layer, one or several middle (hidden) layer(s), and an output layer. All of the neurons of a layer are connected to the ones in the next layer, forming a network with complete connections. The primary parameters in modeling the neural network of MLP are: (1) the number of neurons in each layer, (2) the number of layers in the network, and (3) the forcing functions. A regular MLP neural network has three layers36. The first and the third layers are respectively the system inputs and outputs. The middle layer consists of neurons that perform calculations on the inputs. Choosing the number of layers in a neural network is made by trial and error37. From a hydrological simulation standpoint the main idea behind this model is to create a suitable artificial neural network that is capable of accurately converting a set of hydro-climatic variables such as precipitation and temperature as input data into streamflow values. It should be noted that, like most data-driven models, the process of opting for a proper neural network architecture (i.e., selecting the number of layers, number of neurons, and the forcing function) is, for the most part, a trial-and-error procedure.One must objectively evaluate the performance of the hydrological models in order to opt for the setting of a suitable parameter. The root mean square error (RMSE), coefficient of determination (R2), and mean absolute error (MAE) are herein employed to assess the performance of the rainfall-runoff model. They are respectively calculated as follows:$$RMSE = sqrt {frac{{sumlimits_{t = 1}^{N} {(x_{t} – y_{t} )^{2} } }}{N}}$$
(7)
$$R^{2} = left( {frac{{sumnolimits_{t = 1}^{N} {(x_{t} – overline{x} ).(y_{t} – overline{y} )} }}{{sqrt {sumnolimits_{t = 1}^{N} {(x_{t} – overline{x} )^{2} } } .sqrt {sumnolimits_{t = 1}^{N} {(y_{t} – overline{y} )^{2} } } }}} right)^{2}$$
(8)
$$MAE = frac{{sumnolimits_{t = 1}^{N} {left| {x_{t} – y_{t} } right|} }}{N}$$
(9)
where xt , yt, and N denote the simulated value in time step t; the observed value in time step t; and the number data values, respectively. Large errors have a disproportionately large effect on RMSE or MAE.Performance criteriaVarious quantitative measures can be used to assess the performance of water resources systems under different strategies. When it comes to water resources planning and management, perhaps, some of the most common performance criteria are the probability-based performance criteria (PBPC) (i.e., reliability, vulnerability, and resiliency)31,38. In this context, reliability represents the probability of successful functioning of a system; resiliency measures the probability of successful functioning following a system failure; lastly, vulnerability is the severity of failure during an operation horizon39,40. The basic idea behind a performance evaluation attribute is to provide a quantitative measure to describe and assess the performance of a system. In the context of water resources planning and management, these measures have proven time and again that they can be reliable options to evaluate a set of strategic management options objectively (see, e.g.40,41,42,43, and44, just to name a few).Operating policyAny water resources system requires something called the “rule curve,” which determines how water is allocated in a given situation45. A common and effective rule curve when it comes to operation of water resource systems is the standard operation policy (SOP). SOP is a simple, and perhaps best-known real-time operation policy in water resources planning and management46. The core principle here is to minimize the water shortage at the current time step with no conservation policy (e.g., hedging rules) in place. The SOP, as a standard rule curve, determines how the operator acts to control a system at any given state of a reservoir47,48. This rule curve is established as an attempt to balance various water demands including but not limited to flood control, hydropower, water supply, and recreation49. A SOP operating system attempts to release water to meet a water demand at the current time, with no regard to the future. Thus, according to the SOP’s principle, the decision-makers, first allocate the available water to meet the demand of the stakeholder with the highest priority. After this first water demand is fully satisfied, the available water can be used for the next demand. Such an allocation process continues until no water is available.Ethics approvalAll authors accept all ethical approvals.Consent to participateAll authors consent to participate.Consent to publishAll authors consent to publish. More

Meng, L. et al. Proc. Natl Acad. Sci. USA 117, 4228 (2020).CAS 
Article 

Google Scholar 
Wohlfahrt, G., Tomelleri, E. & Hammerle, A. Nat. Ecol. Evol. 3, 1668–1674 (2019).Article 

Google Scholar 
Wortman, S. E. & Lovell, S. T. J. Environ. Qual. 42, 1283–1294 (2013).CAS 
Article 

Google Scholar 
Su, Y. et al. Agri. For. Meterol. 280, 107765 (2020).Article 

Google Scholar 
Smith, I. A., Dearborn, V. K. & Hutyra, L. R. PLoS ONE 14, e0215846 (2019).Article 

Google Scholar 
Richardson, A. D. et al. Nature 560, 368–371 (2018).CAS 
Article 

Google Scholar 
Meineke, E. K., Dunn, R. R. & Frank, S. D. Biol. Lett. 10, 20140586 (2014).Article 

Google Scholar 
Liu, J. et al. Tour. Manag. 70, 262–272 (2019).Article 

Google Scholar 
Li, X. et al. Remote Sens. Environ. 222, 267–274 (2019).CAS 
Article 

Google Scholar 
Wang, S. et al. Nat. Ecol. Evol. 3, 1076–1085 (2019).Article 

Google Scholar 
Feeley, K. J. et al. Nat. Clim. Change 10, 965–970 (2020).CAS 
Article 

Google Scholar 
Li, D. et al. Nat. Ecol. Evol. 3, 1661–1667 (2019).Article 

Google Scholar 
Li, X. et al. Earth Syst. Sci. Data 11, 881–894 (2019).Article 

Google Scholar 
Román, M. O. et al. Remote Sens. Environ. 210, 113–143 (2018).Article 

Google Scholar 
Li, X. et al. Remote Sens. Environ. 215, 74–84 (2018).Article 

Google Scholar  More

Kaplan, J. O., Krumhardt, K. M. & Zimmermann, N. Quat. Sci. Rev. 28, 3016–3034 (2009).Article 

Google Scholar 
Zheng, Z. et al. Proc. Natl Acad. Sci. USA 118, e2022210118 (2021).CAS 
Article 

Google Scholar 
Lewis, S. L. & Maslin, M. A. Nature 519, 171–180 (2015).CAS 
Article 

Google Scholar 
Steffen, W. et al. Anthr. Rev. 2, 81–98 (2015).
Google Scholar 
Ripple, W. J. et al. BioScience 70, 8–12 (2020).
Google Scholar 
Cohen, J. M., Lajeunesse, M. J. & Rohr, J. R. Nat. Clim. Change 8, 224–228 (2018).Article 

Google Scholar 
Vitasse, Y. et al. Biol. Rev. 96, 1816–1835 (2021).Article 

Google Scholar 
Aono, Y. & Kazui, K. Int. J. Climatol. 28, 905–914 (2008).Article 

Google Scholar 
Sparks, T. H. & Carey, P. D. J. Ecol. 83, 321–329 (1995).Article 

Google Scholar 
Ge, Q. et al. J. Geophys. Res. Biogeosci. 119, 301–311 (2014).Article 

Google Scholar 
IPCC Climate Change 2021: The Physical Science Basis (eds Masson-Delmotte, V. et al.) (Cambridge Univ. Press, 2021).Post, E., Steinman, B. A. & Mann, M. E. Sci. Rep. 8, 3927 (2018).Article 

Google Scholar 
Primack, R. B. & Miller-Rushing, A. J. Bioscience 62, 170–181 (2012).Article 

Google Scholar 
Kharouba, H. M. et al. Proc. Natl Acad. Sci. USA 115, 5211–5216 (2018).CAS 
Article 

Google Scholar 
Richardson, A. D. et al. Agric. For. Meteorol. 169, 156–173 (2013).Article 

Google Scholar 
Parker, D. E., Legg, T. P. & Folland, C. K. Int. J. Climatol. 12, 317–342 (1992).Article 

Google Scholar  More

Wang, L.C.H. Ecological, physiological, and biochemical aspects of torpor in mammals and birds in Advances in comparative and environmental physiology (ed. Wang, L.C.H.) 361–401 (Springer, 1989).Humphries, M.M., Speakman, J.R., & Thomas, D.W. Temperature, hibernation energetics, and the cave and continental distributions of little brown myotis in Functional and Evolutionary Ecology of Bats (ed. Zubaid, A., McCracken, G.F., & Kunz, T.H.) 23–37 (Oxford Press, 2005).Hudson, J.W. Torpidity in mammals in Comparative physiology of thermoregulation. (ed. Whittow, G.C., Hudson, J.W., & Deavers, D.R.) 97–165 (Academic Press, 1973).Geiser, F. & Ruf, T. Hibernation versus daily torpor in mammals and birds: Physiological variables and classification of torpor patterns. Phys. Zool. 68, 935–966 (1995).Article 

Google Scholar 
Davis, W. H. & Hitchcock, H. B. Biology and migration of the bat, Myotis lucifugus, New England. J. Mamm. 46, 296–313. https://doi.org/10.2307/1377850 (1965).Article 

Google Scholar 
Speakman, J. R. & Rowland, A. Preparing for inactivity: How insectivorous bats deposit a fat store for hibernation. Proc. Nutr. Soc. 58, 123–131. https://doi.org/10.1079/PNS19990017 (1999).CAS 
Article 
PubMed 

Google Scholar 
Boyles, J. G., Johnson, J. S., Blomberg, A. & Lilley, T. M. Optimal hibernation theory. Mammal Rev. 50, 91–100. https://doi.org/10.1111/mam.12181 (2020).Article 

Google Scholar 
Britzke, E. R., Sewell, P., Hohmann, M. G., Smith, R. & Scott, R. Use of temperature-sensitive transmitters to monitor the temperature profiles of hibernating bats affected with white-nose syndrome. Northeast. Nat. 17, 239–246. https://doi.org/10.1656/045.017.0207 (2010).Article 

Google Scholar 
Halsall, A. L., Boyles, J. G. & Whitaker, J. O. Jr. Body temperature patterns of big browns during winter in a building hibernaculum. J. Mamm. 93, 497–503. https://doi.org/10.1644/11-MAMM-A-262.1 (2012).Article 

Google Scholar 
Johnson, J. S., Lacki, M. J., Thomas, S. C. & Grider, J. F. Frequent arousals from winter torpor in Rafinesque’s big-eared bat (Corynorhinus rafinesquii). PLoS ONE 7, e49754. https://doi.org/10.1371/journal.pone.0049754 (2012).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Jonasson, K. A. & Willis, C. K. R. Hibernation energetics of free-ranging little brown bats. J. Exp. Bio. 215, 2141–2149. https://doi.org/10.1242/jeb.066514 (2012).Article 

Google Scholar 
Day, K. M. & Tomasi, T. E. Winter energetics of female Indiana bats Myotis sodalis. Physiol. Biochem. Zool. 87, 56–64. https://doi.org/10.1086/671563 (2014).Article 
PubMed 

Google Scholar 
Meierhofer, M. B. et al. Winter habitats of bats in Texas. PLoS ONE 14, e0220839. https://doi.org/10.1371/journal.pone.0220839 (2019).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Boyles, J. G. Benefits of knowing the costs of disturbance to hibernating bats. Wildl. Soc. Bull. 41, 388–392. https://doi.org/10.1002/wsb.755 (2017).Article 

Google Scholar 
Frick, W. F. et al. Pathogen dynamics during invasion and establishment of white-nose syndrome explain mechanisms of host persistence. Ecology 98, 624–631. https://doi.org/10.1002/ecy.1706 (2017).Article 
PubMed 

Google Scholar 
Cheng, T. L. et al. The scope and severity of white-nose syndrome on hibernating bats in North America. Conserv. Biol. 35, 1586–1597. https://doi.org/10.1111/cobi.13739 (2021).Article 
PubMed 
PubMed Central 

Google Scholar 
Cryan, P. M., Meteyer, C. U., Boyles, J. G. & Blehert, D. S. Wing pathology of white-nose syndrome in bats suggests life-threatening disruption of physiology. BMC Biol. 8, 1–8 (2010).Article 

Google Scholar 
Cryan, P. M. et al. Electrolyte depletion in white-nose syndrome bats. J. Wild. Dis. 49, 398–402 (2013).CAS 
Article 

Google Scholar 
Frick, W. F. et al. Disease alters macroecological patterns of North American bats. Glob. Ecol. Biogeogr. 24, 741–749. https://doi.org/10.1111/geb.12290 (2015).Article 

Google Scholar 
Bernard, R.F., Willcox, E.V., Parise, K.L., Foster, J.T., & McCracken, G.F. White-nose syndrome fungus, Pseudogymnoascus destructans, on bats captured emerging from caves during winter in the southeastern United States. BMC Zool. https://doi.org/10.1186/s40850-017-0021-2 (2017).Davy, C. M. et al. The other white-nose syndrome transcriptome: Tolerant and susceptible hosts respond differently to the pathogen Pseudogymnoascus destructans. Ecol. Evol. 7, 7161–7170. https://doi.org/10.1002/ece3.3234 (2015).Article 

Google Scholar 
Lilley, T. M. et al. Resistance is futile: RNA-sequencing reveals differing responses to bat fungal pathogen in Nearctic Myotis lucifugus and Palearctic Myotis myotis. Oecologia 191, 295–309. https://doi.org/10.1007/s00442-019-04499-6 (2019).ADS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Bernard, R. F. & McCracken, G. F. Winter behavior of bats and the progression of white-nose syndrome in the southeastern United States. Ecol. Evol. 7, 1487–1496. https://doi.org/10.1002/ece3.2772 (2017).Article 
PubMed 
PubMed Central 

Google Scholar 
Moosman, P.R., Warner, D.P., Hendren, R.H., & Hosler, M.J. Potential for monitoring eastern small-footed bats on talus slopes. Northeast. Nat. https://www.jstor.org/stable/26453719 (2015).Bernard, R. F. et al. Identifying research needs to inform white-nose syndrome management decisions. Cons. Sci. Prac. 2(e220), 2020. https://doi.org/10.1111/csp2.220 (2020).Article 

Google Scholar 
Reynolds, D. S., Shoemaker, K., Oettingen, S. V. & Najjar, S. High rates of winter activity and arousals in two New England bat species: implications for a reduced white-nose syndrome impact?. Northeast. Nat. 24, B188–B208. https://doi.org/10.1656/045.024.s720 (2017).Article 

Google Scholar 
Bernard, R. F., Willcox, E. V., Jackson, R. T., Brown, V. A. & McCracken, G. F. Feasting, not fasting: Winter diets of cave hibernating bats in the United States. Front. Zool. 18, 1–13. https://doi.org/10.1186/s12983-021-00434-9 (2021).Article 

Google Scholar 
Prendergast, B. J., Freeman, D. A., Zucker, I. & Nelson, R. J. Periodic arousal from hibernation is necessary for initiation of immune responses in ground squirrels. Am. J. Phys. 282, 1054–1062. https://doi.org/10.1152/ajpregu.00562.2001 (2002).Article 

Google Scholar 
Dobony, C. A. et al. Little brown myotis persist despite exposure to white-nose syndrome. J. Fish Wild. 2, 190–195. https://doi.org/10.3996/022011-JFWM-014 (2011).Article 

Google Scholar 
Rowley, J. J. & Alford, R. A. Hot bodies protect amphibians against chytrid infection in nature. Sci. Rep. 3, 1–4. https://doi.org/10.1038/srep01515 (2013).CAS 
Article 

Google Scholar 
Verant, M. L. et al. White-nose syndrome initiates a cascade of physiologic disturbances in the hibernating bat host. BMC Physiol. 14, 10 (2014).Article 

Google Scholar 
Verant, M. L. et al. Temperature-dependent growth of Geomyces destructans, the fungus that causes bat white-nose syndrome. PLoS ONE 7, e46280. https://doi.org/10.1371/journal.pone.0046280 (2012).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Brownlee-Bouboulis, S. A. & Reeder, D. M. White-nose syndrome-affected little brown Myotis (Myotis lucifugus) increase grooming and other active behaviors during arousals from hibernation. J. Wildl. Dis. 49, 850–859. https://doi.org/10.7589/2012-10-242 (2013).Article 
PubMed 

Google Scholar 
Campbell, J. Tennessee winter bat population and white-nose syndrome monitoring report for 2018–2019. TWRA Wildlife Technical Report 16-4. http://www.tnbwg.org/2019%20Annual%20Monitoring%20Report.pdf (2019).Langwig, K. E. et al. Sociality, density-dependence and microclimates determine the persistence of populations suffering from a novel fungal disease, white-nose syndrome. Eco. Let. 15, 1050–1057. https://doi.org/10.1111/j.1461-0248.2012.01829.x (2012).Article 

Google Scholar 
Langwig, K.E., et al. Drivers of variation in species impacts for a multi-host fungal disease of bats. Philos. Trans. R. Soc. B Biol. Sci. 371, 20150456. https://doi.org/10.1098/rstb.2015.0456 (2016).Aldridge, H.D.J.N., & Brigham, R.M. Load carrying and maneuverability in an insectivorous bat: A test of the 5% ‘rule’ of radio-telemetry. J. Mamm. 69, 379–382. https://doi.org/10.2307/1381393 (1988)Sikes, R. S. et al. Guidelines of the American Society of Mammalogists for the use of wild mammals in research and education. J. Mamm. 97, 663–688. https://doi.org/10.1093/jmammal/gyw078 (2016).Article 

Google Scholar 
Barclay, R. M. R. et al. Can external radiotransmitters be used to assess body temperature and torpor in bats?. J. Mammal. 77, 1102–1106. https://doi.org/10.2307/1382791 (1996).Article 

Google Scholar 
Turbill, C., & Geiser, F. Hibernation by tree-roosting bats. J. Comp. Physiol. B. 178, 597–605. https://doi.org/10.1007/s00360-007-0249-1 (2008)Park, K.J., Jones, G., & Ransome, R.D. Torpor, arousal and activity of hibernating greater horseshoe bats (Rhinolophus ferrumequinum). Funct. Ecol. 14, 580–588 (2000).Reeder, D. M. et al. Frequent arousal from hibernation linked to severity of infection and mortality in bats with white-nose syndrome. PLoS ONE 7, e38920. https://doi.org/10.1371/journal.pone.0038920 (2012).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
Sirajuddin, P. Vulnerability of tri-colored bats (Perimyotis subflavus) to white-nose syndrome in the southeastern United States. M.S thesis at https://www.proquest.com/dissertations-theses/vulnerability-tri-colored-bats-i-perimyotis/docview/2185672601/se-2?accountid=8361 (2018).R Development Core Team 3.6.1. A Language and Environment for Statistical Computing. http://www.r-project.org (2019).Bates, D., Maechler, M., Bolker, B., & Walker, S. lme4: Linear mixed-effects models using ‘Eigen’ and S4. R Packag. version 1.1-13. ftp://cran.r-project.org/pub/R/web/packages/lme4/lme4.pdf. (2017)Czenze, Z. J. & Willis, C. K. R. Warming up and shipping out: Arousal and emergence timing in hibernating little brown bats (Myotis lucifugus). J. Comp. Physiol. B-Biochem. Syst. Environ. Physiol. 185, 575–586. https://doi.org/10.1007/s00360-015-0900-1 (2015).Article 

Google Scholar 
Conover, W.J. Practical Nonparametric Statistics. (Wiley, 1999).Crawley, M.J. The R Book: Second Edition. (Wiley, 2013)Lenth, R., Singmann, H., Love, J., Buerkner, P., & Herve, M. Package ‘emmeans’. R Packag. version 1.7.0. https://cran.r-project.org/web/packages/emmeans/emmeans.pdf (2021)Best, T.L. & Jennings, J.B. Myotis leibii. Mammalian Species 547, 1–6. https://doi.org/10.2307/3504255 (1997)Frank, C.L., Herzog, C.,. Laske, J.P., & Cardino V. The role of skin temperature in the resistance of Myotis leibii to white-nose syndrome presented at the 49th Symposium of the North American Society for Bat Research, Kalamazoo, Michigan https://www.nasbr.org/resources/docs/meetings/year/2019/NASBR-2019-Abstracts-Draft-191001.pdf (2019)Johnson, J.S., Scafini, M.R., Sewall, B.J., & Turner, G.G. Hibernating bat species in Pennsylvania use colder winter habitats following the arrival of white-nose syndrome in Conservation and Ecology of Pennsylvania’s Bats (ed. Butchkoski, C.M., Reeder, D.M., Turner, G.G., & Whidden, H.P) 181–199 (The Pennsylvania Academy of Science, 2016)Haase, C.G., et al. Body mass and hibernation microclimate may predict bat susceptibility to white‐nose syndrome. Eco. Evo. 11, 506–515 https://doi.org/10.1002/ece3.7070 (2021)Veilleux, J.P. A Noteworthy Hibernation Record of Myotis leibii (Eastern Small-footed Bat) in Massachusetts. Northeast. Nat. 14, 501–502 https://doi.org/10.1656/1092-6194(2007)14[501:ANHROM]2.0.CO;2 (2007)Boyles, J. G., Dunbar, M. B. & Whitaker, J. O. Activity following arousal in winter in North American vespertilionid bats. Mamm. Rev. 36, 267–280. https://doi.org/10.1111/j.1365-2907.2006.00095.x (2006).Article 

Google Scholar 
Webb, P. I., Speakman, J. R. & Racey, P. A. How hot is a hibernaculum? A review of the temperatures at which bats hibernate. Can. J. Zool. 74, 761–765. https://doi.org/10.1139/z96-087 (1996).Article 

Google Scholar 
Jackson, R.T., Willcox, E.V., Zobel, J.M., & Bernard, R.F. Hibernation behavior of four bat species with differing susceptibility to white-nose syndrome. In review (2021).Fujita, M. S. & Kunz, T. H. Pipstrellus subflavus. Mamm. Spec. 228, 1–6. https://doi.org/10.2307/3504021 (1984).Article 

Google Scholar 
Bohn, S. J. et al. Evidence of ‘sickness behaviour’ in bats with white-nose syndrome. Behaviour 152, 981–1003. https://doi.org/10.1163/1568539X-00003384 (2016).Article 

Google Scholar 
Tuttle, M. D. Status, causes of decline, and management of endangered gray bats. J. Wild. Manag. 43, 1–17. https://doi.org/10.2307/3800631 (1979).Article 

Google Scholar 
Harvey, M.J., Altenbach, J.S., & Best, T.L. Bats of the Eastern United States. (JHU Press, 2011).Klüg-Baerwald, B. J., Lausen, C. L., Willis, C. K. & Brigham, R. M. Home is where you hang your bat: Winter roost selection by prairie-living big brown bats. J. Mamm. 98, 752–760. https://doi.org/10.1093/jmammal/gyx039 (2017).Article 

Google Scholar 
Sandel, J. K. et al. Use and selection of winter hibernacula by the eastern pipistrelle (Pipistrellus subflavus) in Texas. J. Mamm. 82, 173–178. https://doi.org/10.1644/1545-1542(2001)082%3c0173:UASOWH%3e2.0.CO;2 (2001).Article 

Google Scholar 
Hayes, J. P., Ober, H. K., Sherwin, R. E. Survey and monitoring of bats in Ecological and behavioral methods for the study of bats (ed. Kunz, T. H., Parsons, S.). 120–129 (JHU Press, 2009). More

Study areaThe study area covers the total land area of Thailand. Where it is useful, we divided Thailand into six regions (Fig. 2a), the names and boundaries of which are widely used, although they have no official administrative status. We focused on the elements of Thailand’s protected area system that were concerned principally with the in-situ conservation of biodiversity: existing and proposed National Parks, Wildlife Sanctuaries, Non-hunting Areas, and Forest Parks, covering 111, 201 km2 or 21.7% of the country’s land area37 (Fig. 1).Environmental dataA set of environmental variables that were expected to be directly or indirectly related to species distributions in Thailand was used to model suitable habitat in the present and future (Supplementary Material Table S1). These variables were chosen to encompass ecologically relevant variables and enable consistent comparison between species, regardless of species-specific preferences. GIS layers for the whole of the study area were compiled using a variety of data sources at 1-km2 resolution. For variables originally at higher than 1-km resolutions, we used the plus function in ArcMap to combine them with a mask of the study area to use the mask dimensions for all cells.The physical variables, altitude, slope, aspect, and soil pH are widely used in species distribution modeling. Slope and aspect have biologically significant impacts on both temperature and rainfall at these latitudes8 and are particularly important at the poleward margins of species ranges where species may be confined to one aspect. Slope also affects soil maturity and depth. Soil pH is a consistently measured soil variable that broadly correlates with fertility in tropical soils8. Additional soil variables, particularly soil phosphorus, have been shown to be important filters of plant species distributions in the tropics38, but they are not available for Thailand with a useful accuracy and spatial resolution. Altitudes were downloaded from the CGIAR-Consortium for Spatial Information, CGIAR-CSI version 4.1. Slope and aspect were generated by using surface tools in ArcGIS. Soil pH was extracted from ISRIC-World Soil Information version 2.0.Unlike the temperate zone, where tolerances of winter cold and requirements for summer warmth dominate plant and animal distributions, our understanding of how tropical climates filter species distributions is still weak38,39. In Thailand, as in most of the tropics, there are two major climatic gradients which correlate with changes in species composition: a rainfall gradient in the lowlands, along which total rainfall declines and the length of the dry season increases, and a gradient of steadily declining temperature with elevation7. There is no simple relationship between elevation, and thus temperature, and rainfall. An additional complication is that temperature seasonality may be significant in northern Thailand (north of c.18° N), where cooler winters reduce dry-season water stress and extreme low temperatures at high altitudes may exceed physiological tolerances. We therefore chose 8 bioclimatic variables (Supplementary Material Table S1) related to precipitation and temperature, and their seasonality, all of which have previously been used in species distribution modelling in this region9,40. These are available at a resolution of 30 arc sec (approximately 1 km at the equator) from WorldClim ver. 1.4 based on averages of 1970–1990. These variables are available from the same source (and downscaled using the same methods) for the future climate projections.Vegetation structure is an additional major influence on plant and animal distributions in the tropics, both in intact natural vegetation38,39 and when the original vegetation has been degraded or cleared8. Vegetation structure was represented through the inclusion of two continuous variables, percentage forest cover and tree density, as most of the modelled species are known to be sensitive to both the presence of forest and the degree of intactness of the tree cover9. Mean tree density per km2 was extracted from Crowther et al.41 version 2 and percentage coverage of forest per km2 was extracted from the European Space Agency (ESA) GlobCover Version 2.3.Note that the mechanistic basis of the correlations between all these variables and the current distributions of tropical plants and animals are rarely known. Temperature has a direct physiological impact on all organisms, and water supply may be seasonally limiting for plants and some amphibians, but indirect links through biotic interactions are expected to be more important in the tropics, including pest pressure on plants38 and food supply for animals39. Competition is probably also important in shaping local species assemblies. For future projections, we assumed that temperature and precipitation were changing, and that other variables (topography, soil, and vegetation) were stable, so our analysis represents the impacts of climate alone. For 2070, we used the same variables projected by three CMIP5 Earth System Models, CNRM-CM5, GFDL-CM3 and HadGEM2-ES, which have been previously used in Southeast Asia9,42 and in Thailand7. We used two Representative Concentration Pathways, RCP2.6 and RCP8.5, representing low and high greenhouse-gas concentration scenarios, respectively, and thus the potential range of radiative forcing by the end of the century43. RCP2.6 is consistent with meeting the Paris Agreement’s 2 °C global warming target.Species occurrence dataMany locality records for vertebrates were supplied by the Department of National Parks, Wildlife and Plant Conservation (DNP). Trained DNP staff walked along trails throughout the protected areas in Thailand during 2017–2018. They recorded 271,695 locations for 70 mammal species, 18 locations for 3 amphibian species, 318 locations for 18 reptile species, and 43,057 locations for 65 bird species44. We supplemented this with data downloaded from the Global Biodiversity Information Facility (GBIF, https://www.gbif.org/) for 1960–2019 for amphibians (2063 localities for 86 species)45, reptiles (1722 localities from 196 species)46, mammals (2508 localities from 191 species)47, and birds (1,559,222 localities from 884 species)48. More than 95% of the bird records from GBIF were identified as coming from eBird49, which is popular among birders in Thailand. For plants, we used occurrence data from the DNP’s forest resource inventory project from 221 plots, including 24,605 localities for 363 species, the DNP’s Forest Herbarium, including 227 localities for 141 species, and locations for 12 rare and endangered forest species collected from all over Thailand. We also downloaded data from the Botanical Information and Ecology Network (BIEN, https://bien.nceas.ucsb.edu/bien/), including 7209 localities for 1422 species.We removed suspect records (coordinate issues, name problems, etc.), duplicates from the same locality (i.e., more than one individual of the same species recorded in a cell), and species with  0.5 as adequate, but since only five SDMs out of the 1457 generated in this study had values lower than this (0.3–0.5), we retained all the models.Assessment of climate change impactsThe estimated current distribution for each species from Maxent was used as the baseline for comparison with projected distributions of suitable habitat for these species by 2070, under the two emission scenarios and three ESMs, and with and without unlimited dispersal into newly available habitat. We then assessed the impacts of climate change, both on the spatial distribution of individual species and on the pattern of species richness. To generate a species richness map, the binary habitat suitability maps for all species were stacked to produce a consolidated map, which showed the number of species for each 1 km grid cell, and then classified them into five classes (lowest, low, moderate, high, and highest), using the mean ± standard deviation as a break class40.Current and future maps were then compared for each species to calculate the change in species richness, and contingency tables were generated containing the numbers of cells (each of 1 km2) in each richness class. Suitable habitat areas were calculated for the current climate and projected for the future climate. For each species we estimated gained habitat as the areas that will become suitable for a species in future under that scenario, lost habitat as the areas currently predicted as suitable now but projected to become unsuitable under future climatic change, and stable habitat as the areas predicted as suitable now which will remain suitable into the future.We then assessed the vulnerability of each species by estimating the projected change in its range over the next 50 years and using a criteria-based approach, which combined the mean of the suitable habitat area (interpreted as equivalent to extent of occurrence) in the three models and a simplified version of the IUCN Red List criteria51. For 2070, we modified criterion A3(c) as follows; Extinct (Ex) species are projected to lose 100% of suitable habitat by 2070, Critically Endangered (CR) species are projected to lose over 80%, Endangered (EN) species are projected to lose 50–80%; Vulnerable (VU) species are projected to lose 30–50%, Near Threatened (NT) species are projected to lose  More

Collection of wolf hair samplesHair samples were collected by researchers from opportunistically found-dead wolves upon standard necropsy (all the Alpine and part of the Iberian samples) or in the field (all the Dinaric-Balkan and most of the Iberian samples), or from legally harvested wolves (only in the Scandinavian population). At the time of sample collection, wolves were legally harvested in Sweden, Slovenia, and Spain, and under total protection in Portugal and Italy. Hair samples were collected from four body regions, when possible: lumbar (n = 133), dorsal cervical (n = 66), tail (n = 33) and ventral thorax (n = 27) (Tables S1 and S2). The hair was cut as close as possible to the skin with scissors to avoid collecting hair follicles, but in some samples, hairs were pulled from the carcass. Samples were stored at room temperature in paper envelopes. Age, sex, date, and cause of death/capture, geographical location, body mass, and total length were obtained for most of the wolves.Age was estimated by the dental eruption and wear or cementum age analysis and classified as ‘juveniles’ ( 2 years)40, or ‘unknown’. Sex was assessed by inspection of genitalia. Causes of death were classified as ‘acute’, likely lasting minutes to hours (vehicle accident and legal or illegal shooting); ‘subacute’, likely lasting hours to days (drowning, poisoning, trapping and intraspecific aggression); ‘chronic’, likely lasting several weeks (infectious diseases—canine distemper, canine parvovirosis, leptospirosis; sarcoptic mange; or neoplastic diseases) or ‘unknown’. Total length was obtained by measuring with metric tape (1 mm precision) the distance from snout to the distal end of the last tail vertebrae. The body mass was measured with 100 g precision with scales.The detailed protocol for the handling of wolves live trapped in the scope of ecological and conservation studies (n = 7, all from the Iberian population) has been previously described5. Traps were monitored twice every day, in the early morning and late afternoon, hence the duration of restraint after capture was unknown for 8 wolves, potentially up to 12 h. Trap-alarms were deployed in the capture of 2 wolves, with 41 and 70 min intervals between activation of the alarm and administration of the drugs. Live trapping was conducted under permits issued by the nature conservation authorities of Portugal (Instituto de Conservação da Natureza e das Florestas: 338/2007/CAPT, 258/2008/CAPT, 286/2008/CAPT, 260/2009/CAPT, 332/2010/MANU, 333/2010/CAPT, 336/2010/MANU, 26/2012/MANU, and 72/2014/CAPT) and Spain (Dirección Xeral de Conservación da Naturaleza, Xunta de Galicia: E-0020/13-PNPE, 095/2013; Consejería de Medio Ambiente, Principado de Asturias: 31/08/2017-BOPA 05/09/17) and according to European Union directives on the protection of animals used for scientific purposes (Directive 2010/63/EU) and international wildlife standards41,42. The study was undertaken in compliance with the ARRIVE guidelines43.Cortisol extractionThe protocol for the extraction of cortisol from the hair was adapted from previously described procedures15,27. Forty mg of guard hairs were separated from the undercoat and placed in 15 ml falcon tubes. Hair follicles were cut whenever found in the sample. For each sample, the length of three intact hairs was recorded. The samples were washed twice with 40 µl of distilled water/mg hair and three times with the same amount of isopropanol. In each washing step, the samples and washing solution were vortexed, the supernatant discarded, and the hair dried using clean paper towels. After the final wash, samples were dried overnight at room temperature and 30 mg of hair cut into a 2 ml polypropylene screw cap plastic tube with five 4 mm steel beads added to each tube.The hair was ground to a fine powder in a FastPrep sample homogenizer (MP Biomedicals, USA) for four times 1 min at 6.0 m/s. 50 µl methanol/mg hair were added to each sample and sonicated for 30 min at 50 Hz at 50 °C. The samples were incubated for 18 h at 50 °C in an orbital shaker at 160 rpm, centrifuged for 15 min at 14,000g at 20 °C, and 1000 µl of supernatant was collected to a screw cap glass chromatography vial and dried at room temperature in a gentle stream of nitrogen gas. Due to restrictions on laboratory use during the SARS-Cov-2 pandemic, some batches of samples were instead evaporated overnight on a suction hood. This unexpected change in the methanol evaporation protocol was recorded and accounted for in the statistical analysis.Cortisol quantificationA commercial competitive ELISA kit (Cortisol free in Saliva ELISA, Demeditec, Germany) was used to quantify the concentration of cortisol, following the manufacturer’s instructions. The kit plate wells are provided coated with polyclonal rabbit antibody against cortisol, and cortisol-horseradish peroxidase was used as conjugate. According to the manufacturer, the cross-reactivity of the test to selected steroids is low (Table S3), the intra-assay variation is 3.8–5.8% and the inter-assay variation is 6.2–6.4%. Samples, standards, and controls were tested in duplicate.The 4-parameter standard curve was calculated from the log-transformed cortisol concentration of the standard solutions and their measured OD45044. Standard curves were estimated using the software GraphPad Prism 6.04 (GraphPad Software, La Jolla, California USA), and yielded an average R2adjusted = 0.991 (range 0.968–0.999). The cortisol concentration of the reconstituted samples was estimated from the standard curve and converted to cortisol concentration as picograms (pg) of cortisol/mg of guard hair.Intra and inter-assay coefficients of variation were estimated for six ELISA assays of 37–40 samples each. The low and high controls included in the kit were used to estimate the inter-assay coefficient of variation and the duplicate runs of each sample were used to estimate the intra-assay coefficient of variation. Linearity was assessed by two-fold dilutions (1:1, 1:2, 1:4 and 1:8) of 4 extracted samples, comparing the expected and observed concentrations. Recovery was assessed by spiking 6 ground hair samples with known concentrations of cortisol (50, 25, 12.5, 6.25 pg/mg, and no spiking), comparing the expected and observed concentrations.The intra-assay coefficient of variation of the ELISA assays ranged from 6.50 to 9.97% (average 7.66%). The inter-assay coefficient of variation was 11.54% for the low concentration controls and 9.08% for the high concentration controls (average 10.31%). Assay linearity was 91% for the 1:2 dilution, 103% for 1:4, and 117% for 1:8 (average 103%). The recovery of cortisol averaged 94%, being 73% for the 50 pg/mg spiked samples, 74% for 25 pg/mg, 95% for 12.5 pg/mg, and 113% for 6.25 pg/mg.Determinants of hair cortisol concentrationThe potential determinants of HCC investigated included wolf intrinsic variables: sex, age, body condition, body structural size, month of death/capture, and wolf population. The scaled mass index was selected as a measure of body condition45 and estimated from the log-transformed body weight (g) and total length (mm). Log-transformed total length was used as an indicator of body structural size46. Samples were assigned to the Iberian, Alpine, Dinaric-Balkan, or Scandinavian wolf populations16 from the geographical location of the death or live-trapping sites (Fig. 1).The relationship between HCC and additional variables related to the sampling procedure or to the work conducted in the laboratory (length of hair used for cortisol extraction, sample storage time, body region, cause of death/capture, and methanol evaporation protocol), herein referred to as methodological variables, was also investigated as potential confounding variables. Sample storage time was the period in months between death/capture and cortisol extraction. In those samples for which only the year of death was available, 30 June was assigned as the date of death, solely to estimate storage time. All continuous variables were standardized to their z-scores.Statistical analysisFirst, the effect of body region was investigated by a linear mixed model with HCC as the dependent variable, and the independent variables body region, as a categorical fixed effect, and individual wolf, as a random effect. The lumbar region was set as the reference class as it was the most represented in our sample (Table S1). Data from 27 wolves for which samples were available from all 4 body regions were used in this analysis. Four outliers in the dataset violated the assumption of normality in the residuals of the model comparing HCC across body regions (Fig. S1A) and were excluded from this model’s dataset (Fig. S1B).Second, the effect of intrinsic and methodological variables on HCC from the lumbar body region was investigated by another linear mixed model with sex, age, body condition, body structural size (standardized log-transformed total length), cause of death/capture, wolf population, hair length, sample storage time, and methanol evaporation protocol as fixed effect independent variables. The month of death/capture was included as a random effect. Reference classes for the categorical variables were set as female, adult, acute death, Iberian population, and methanol evaporation by nitrogen gas stream. Two outliers in the dataset violated the assumption of normality in the residuals of the model (Full model, Table S4) and were excluded from this analysis (Fig. S1C,D).The goal of this analysis was to assess the relationship between HCC and wolf intrinsic variables, controlling for the potential confounding effect of the methodological variables. Starting from the full model (Table S4), models including all possible combinations of variables were ranked by their AICc using the package “MuMIn”47 in R 3.6.148. The most supported model was selected for inference and models with ΔAICc  More