Ecology

Clothianidin exposure in field experiments
The experimental bee yard was situated in a suburban area in Kirchhain (Germany) with sufficient natural pollen sources in late summer, where the honeybees were allowed to forage freely. The field experiment took place with 20 Apis mellifera carnica colonies in Kirchhainer mating nuclei (25.5 × 19.8 × 17 cm), established with sister queens (A. m. carnica breeding line, mother: 17:27:20:11) mated on an island mating station (Norderney) and 180 g of worker bees each. We tried to reduce the variance between colonies by using young sister queens from an island mating station. Furthermore, we used small hives to minimize food competition so that we were able to set up the colonies within a rather small area of 20 by 20 m. In detail, the worker bees originated from four healthy colonies of the institute, located at an apiary of the institute. In order to get a homogenous composition, worker bees of brood combs and honeycombs of two colonies were shaken into a box and mixed thoroughly, before they were uniformly distributed into 10 mating nuclei. This procedure was repeated to fill 20 mating nuclei. Subsequently, the mating boxes were placed in a cool room (12 °C) for three day, to acclimatize, before they were transported to the island mating station Norderney. The queens started to lay eggs at June 25, 2014. The clothianidin exposure started at July 28. Thus, the colonies containing freshly mated, egg-laying queens were allowed to establish and build up an intact brood nest for four weeks. Two days before clothianidin exposure (sampling day 0, S0; SFig. 1), the colony strength was assessed and the treatment groups were randomly assigned to the hives such that differences between treatment groups were minimized with respect to the strength of their colonies. Every week, from July 30 to September 10, each colony received 400 mL of Apiinvert (Südzucker AG, Mannheim, Germany) sugar syrup (39% w/v fructose, 31% w/v saccharose and 30% w/v glucose) spiked with clothianidin (1, 10 or 100 µg/L). Control colonies received sugar syrup containing the same concentration of the solvent (water) as the clothianidin-treated groups. The syrup was fed in zip-look bags placed inside the food chamber containing a climbing aid. After 1 week, the leftovers were removed and weighed to record food consumption. For all four experimental groups, the analyzed clothinidin levels were close to the target concentrations (Suppl. Table 1). Environmental data were recorded during the study period using a USB data logger (EL-USB-2, Lascar Electronics Ltd., temperature accuracy ± 0.5 °C, relative humidity accuracy ± 3%) located under the colonies.
Sampling
During weeks 3 and 7, random samples of worker bees, larvae and worker jelly were taken from each colony. In detail, 20 randomly chosen worker bees located on a brood comb and five larvae (larval stage: day 7 or 8 after egg laying) were immediately frozen for chemical analysis. To collect worker jelly, extra thick blotting paper (Protean, Xi size; Bio-Rad, Hercules, CA, USA) was cut into strips, cleaned in pure ethanol and acetone (Carl Roth, Karlsruhe, Germany) and dried in a heating cabinet at 80 °C. Each strip was inserted into five brood cells containing a small larva (developmental day 4–5) to suck up the worker jelly and the strips were immediately frozen.
To document brood development, each side of each comb was photographed within an empty hive box transformed into a photo box containing a digital camera (Canon PowerShot A1000 IS, Tokyo, Japan) and a ring-flash (Aputure Amaran AHL-C60 LED, Shenzhen, China). Brood documentation took place every week. The colonies were sampled starting with the control and then from the lowest to highest concentration of clothianidin to minimize the risk of carryover.
HPG size measurements
To obtain worker bees of a defined age, single frames of late-stage capped brood (Binder, Tuttlingen, Germany) were brought to the laboratory and incubated in the dark at 32 °C, with humidity provided by open water jars. The frames with worker brood were collected from two full-sized colonies, which were regularly inspected for symptoms of disease and tested for Chronic bee paralysis virus, Deformed wing virus, Acute bee paralysis virus, and Sac brood virus65. Newly-emerged bees (≤ 24 h) were collected, color marked, and transferred to the experimental colonies on July 28 (15 marked bees per colony). During the second week of exposure, the marked bees (Suppl. Table 2) were removed from the colonies after 12 days in the hive and immobilized on ice. The HPGs were dissected in ice-cold phosphate-buffered saline (PBS, pH 7.4). The specimens were fixed in formaldehyde (4% in PBS, Carl Roth), rinsed three times in PBS, and mounted in Aquapolymount (Polysciences, Eppelheim, Germany). Three pictures of each gland (only one gland per bee) were photographed at 400x magnification using a phase contrast/fluorescence microscope (Leica DMIL, Leica camera DFC 420C) and LAS v4.4 image-capturing software (Leica Microsystems, Wetzlar, Germany). To measure the size of each gland, the diameters of 30 acini per bee were measured using ImageJ v1.49o (http://rsb.info.nih.gov/ij/index.html).
High-performance thin-layer chromatography
HPTLC was chosen as sensitive analytical method in order to detect qualitative and semi-quantitative differences in the composition of brood food and larvae of dosed hives. Individual larvae (one per colony) differing in their weights (200–500 mg/larva) were solely macerated and extracted in 1 mL n-hexane ( > 99% pure, Rotisolv, Carl Roth) in an ultrasonic bath for 1 min and then vortexed for 1 min. For the analysis of worker jelly (from three cells per colony), adsorptive filter strips (Sugi strips, Kettenbach, Eschenburg, Germany) were cut in half and one part was dipped in brood combs until maximal absorption of the material, whereas the other strip was used as background control strip, and both strips were extracted as above. The supernatants of larvae and worker jelly were transferred to a fresh vial 1000 μL isopropylacetate/methanol (3/2, v/v). After maceration for 1 min, the mix was vortexed for 1 min and the supernatant was transferred to another vial. Between extractions, the samples were cooled on ice and stored at –20 °C.
The isopropylacetate/methanol extracts (20 µL/band for worker jelly and 7 µL/band for larvae) were sprayed onto the silica gel 60 F254 HPTLC plate (Merck, Darmstadt, Germany) using an Automatic TLC Sampler 4. The plate of worker jelly was developed with a mobile phase consisting of chloroform/methanol/water/ammonia (30/17/2/1, v/v/v/v, all Carl Roth66 and the plate of larvae with an 8-step gradient development based on methanol, chloroform, toluene and n-hexane66,67. After drying in a stream of cold air for 2 min, the plate images were documented at UV 366 nm using the TLC Visualizer. For derivatisation, the chromatogram was dipped into the primuline solution (100 mg primuline in 200 mL acetone/water, 4/1 v/v, Sigma-Aldrich, Steinheim, Germany) at an immersion speed of 2.5 cm/s and an immersion time of 1 s, dried and derivatised as before. The data were processed using winCATS version 1.4.2.8121 (all instrumentation from CAMAG).
The bacterium Aliivibrio fischeri (NRRL-B11177, strain 7151), obtained from the German Collection of Microorganisms and Cell Cultures (DSMZ, Leibniz Institute, Berlin, Germany), was used to assess a non-targeted, broad range of effective substances within the worker jelly. HPTLC plates were developed as described above, neutralized and immersed in a bacterial suspension, prepared according to DIN EN ISO 11,348–1 845 at an immersion speed of 3 cm/s and an immersion time of 2 s. The bioluminescence of the wet bioautogram was recorded in an interval of 3 min over 30 min using the BioLuminizer (CAMAG).
Brood assessment
The free extension of the open source program ImageJ, Fiji (http://fiji.sc/Fiji) was used to count all cells with eggs, larvae, or sealed brood, for every colony on every sampling day. The photos of a single comb from different sampling days were aligned in the program and all brood cells were counted. Because the colonies had no wax foundations, some cells on the edges of these naturally-built combs were at an unfavorable angle. Therefore, most but presumably not all cells with brood were visible in the photos. This uncertainty was similar in all hives. To estimate the brood survival, we first tracked the development of individual eggs, but found a high mortality rate even in control colonies. Therefore, we tracked the development of individual larvae (day 4–5 after egg laying) over 4 weeks (= sampling weeks). Brood survival was estimated for two periods during the experiment (weeks 1–4 and 3–7).
Modeling of demographic compensation
The aim of the model was to estimate the number of new larvae that must hatch each week in order to maintain a stable number of larvae up to 7 weeks of age given the survivorship schedule of larvae week by week. Let h denote the number of eggs that hatch into larvae each week, and let the probability that any individual larva dies in each successive week be mi, where i takes values in the set {1, 2, …, 7} to indicate each of seven successive weeks, after which we assume that surviving larvae pupate. We used a demographic matrix model68 to describe the state of the population of larvae each week as follows. Let lx denote the number of larvae in the colony that are aged x weeks post-hatching and let mx denote their per capita weekly mortality rate. The population of larvae is distributed into seven age classes and we also assign a class to queens, which give rise to larvae by producing eggs. Using the Lefcovitch matrix approach, the larval population of a colony can thus be viewed as a state vector nt whose week-by-week change is the product of a matrix A and the population state vector nt, as shown in Eq. (1):

$$An_{t} = left[ {begin{array}{*{20}c} 0 & 0 & 0 & ldots & 0 & {varvec{h}} \ {(1 – {varvec{m}}_{1} )} & 0 & 0 & ldots & 0 & 0 \ 0 & {left( {1 – {varvec{m}}_{2} } right)} & 0 & ldots & 0 & 0 \ vdots & vdots & vdots & ddots & vdots & vdots \ 0 & 0 & 0 & ldots & {left( {1 – {varvec{m}}_{7} } right)} & 0 \ 0 & 0 & 0 & ldots & 0 & 1 \ end{array} } right]left[ {begin{array}{*{20}c} {{varvec{l}}_{1} } \ {{varvec{l}}_{2} } \ {{varvec{l}}_{3} } \ vdots \ {{varvec{l}}_{7} } \ {varvec{Q}} \ end{array} } right]$$
(1)

The lowermost element of nt is the number of queens (Q) which is here set to Q = 1 for all models, and the total number of larvae in the colony at any time is ({varvec{L}}_{{varvec{t}}} = mathop sum nolimits_{{varvec{x}}} {varvec{l}}_{{varvec{x}}}). Given the observed mortality schedule (mx), we used the model to solve for the value of h that produces a stable value for Lt, as shown in Eq. (2):

$$An_{t} = n_{t}$$
(2)

We determined the mortality schedule by observing the survivorship of a larval cohort. For example, if a cohort of St larvae was originally marked at time t and, of these, St+1 survived until the following week, the per capita weekly mortality rate would be estimated as shown in Eq. (3):

$${varvec{m}}_{{varvec{t}}} = 1 – frac{{{varvec{S}}_{{{varvec{t}} + 1}} }}{{{varvec{S}}_{{varvec{t}}} }}$$
(3)

BEEHAVE simulations
To predict the impact of clothianidin on colony development in standard hives we used BEEHAVE, a honeybee model that simulates colony dynamics and agent-based foraging in realistic landscapes44(http://beehave-model.net/). Although BEEHAVE does not explicitly allow the incorporation of pesticides, the effect of pesticides on behavior and mortality can nevertheless be addressed59,60. BEEHAVE simulates the development of a single honeybee colony, starting with 10,000 foragers on January 1. Colony dynamics are based on a daily egg laying rate, with the developmental stages eggs, larvae, pupae, and adults (in drones) or in-hive workers and foragers (in workers). The brood needs to be tended by in-hive bees, and the larvae additionally need to be fed with nectar and pollen. Foragers can scout for new food sources or collect nectar and pollen from sources already known. Successful foragers can recruit nestmates to the food source. Mortality rates depend on the developmental stage and the time spent on foraging. The colony dies if it either runs out of honey or if the colony size falls below 4000 bees at the end of the year. Swarming may take place when the brood nest grows to more than 17,000 bees before July 18. Under default conditions, two food sources are present at distances of 500 or 1,500 m. Daily foraging conditions are based on weather data from Rothamsted, UK.
We ran two sets of simulations: (A) We first set up BEEHAVE to mimic our experimental nucleus colonies. We then determined how the protein content of the jelly produced by nurses (ProteinFactorNurses) had to be modified to replicate the larval mortality we observed in our empirical data. (B) We then set up BEEHAVE under default conditions but modified ProteinFactorNurses according to the results from the previous simulations to assess the impact of clothianidin exposure under more realistic conditions.
To mimic the experimental nucleus colonies (simulation A), the maximum honey store (MAX_HONEY_STORE_kg) was reduced to 0.77 kg and the maximum size of the brood nest (MAX_BROODCELLS) was reduced to 250. Furthermore, HoneyIdeal was set to ‘true’, so that even though the honey store was small it was filled every day, reflecting the feeding of the experimental bees. In contrast, PollenIdeal was set to ‘false’, because the experimental colonies still had to forage for pollen. On day 209 (July 28), we set the number of pupae to 100 and the number of workers to 650, similar to the experimental colony sizes. During the exposure to clothianidin between days 211 (July 30) and 253 (September 10), the protein content of the jelly fed to the larvae (ProteinFactorNurses) was modified by the new variable ProteinNursesModifier_Exposed. We tested for ProteinNursesModifier_Exposed values from 0.6 to 1 in steps of 0.01. The main output of the simulation was the survival of the brood, calculated from the brood cohort sizes aged 19 days divided by the sizes of these cohorts when they were 3 days old. We calculated the mean brood survival over 30 replicates, using the last 10 cohorts only (i.e. those reaching the age of 19 days between September 1 and 10). Those parameter values for ProteinNursesModifier_Exposed resulting in brood survival most similar to the experimental brood survival were then chosen to represent the clothianidin concentrations of 100, 10 and 1 µg/L.
To assess the impact of clothianidin on standard colonies under more realistic conditions (simulation B), we ran BEEHAVE under default settings but reduced the protein content of the jelly (ProteinFactorNurses) during times of exposure. We assumed that colonies would be exposed when rapeseed plants are flowering, defined in the model as the period between days 95 (April 5) and 130 (May 10). ProteinNursesModifier_Exposed values representing the tested concentrations of clothianidin were derived from the previous set of simulations, and for the control we set ProteinNursesModifier_Exposed to 1 (i.e. no effect of clothianidin). Swarming was either prevented or allowed, in which case the simulation followed the colony remaining in the hive.
Statistical methods
Statistical analysis was carried out using R v3.4.269, including the add-on packages lme470 for linear mixed-effects models, pbkrtest71 for testing fixed effects in mixed-effects models, parallel69 to increase computational power, RLRsim72 for testing random effects in mixed-effects model, multcomp73 for multiple comparisons, and lattice74 for various graphical displays. We used linear mixed-effects models for one- and two-factorial analysis of variance (ANOVA) or regressions as indicated and where necessary. In those models, the colonies were modeled as random effects to reflect the (longitudinal) grouping structure in the data. For the analysis of larval survival, we used a logarithmic transformation of the proportion of surviving larvae. When testing fixed effects in mixed-effects models, we used the Kenward-Roger method (and double-checked the results by comparing them with parametric bootstrap values). Where sufficient, we simplified the analysis using linear (fixed-effects) models. Model diagnostics were performed for all fitted models using qualitative tools such as normal q-q-plots for residuals and plots of residuals versus fitted values to assess the validity of model assumptions like homoscedastic normality. Dunnett’s test (or customized contrasts as appropriate) was used for multiple comparisons in post hoc analysis, and P values were appropriately adjusted for multiple testing within well-defined test families using the single-step or Westfall’s method. Model and analysis details, model diagnostic graphs, and further information are available in the supplemental statistical report. More

A total of 779 persons opened the survey, 486 persons started responding to it, and 308 persons from 40 countries submitted their responses. Seven of these 308 respondents had never heard about biases, and the remaining 301 persons (i.e. 38.6% of those who opened the survey) answered our questions about their attitude to biases.
Nearly all (98%) scientists who responded to our survey were aware of the importance of biases in science. Among these, 33% reported ‘very well’ for their awareness on this topic, 52% classified their awareness as ‘well’ and 13% as ‘poor’. Most of respondents learned about biases from their university courses (36%), from contacts with colleagues (22%) or from scientific literature (20%). Among the different kinds of biases, the best known was observer/observation bias (82%), followed by publication bias (71%) and selection bias (70%); confirmation, reporting/presentation, researcher, measurement, geographic and funding biases were known to 50–60% of respondents (see Appendix S1). Among the seven suggested definitions of biases (see Appendix S1), ‘the tendency to search for, interpret, and publish information in a way that confirms one’s pre-existing beliefs or hypotheses’ corresponded to the understanding of biases by 80% of the respondents, and ‘preferential publication of statistically significant results’ was ranked the second, being selected by 61% of respondents.
In the opinion of our respondents, the stages of scientific research differed in their susceptibility to biases (χ220 = 90.0, P  planning/designing the study  > publishing the outcomes  > reporting the outcomes  > analysing the results  > implementing the study. Similarly, the types of publications were considered prone to biases to different extents (χ216 = 86.8, P  studies based on observational data  > studies based on modelling  > studies based on experiments  > meta-analyses.
Most scientists thought that the severity of the impact of biases on science in general and on their particular research field was medium or high, and a few considered that it was negligible. At the same time, our respondents estimated the impact of biases on their own studies as high almost three times less frequently and as negligible seven times more frequently when compared with their estimates of the impact of biases on other studies within their own research field (Fig. 1).
Figure 1

Impact of biases on science in general, on one’s own research field and on one’s own studies, as estimated by 301 respondents. Bars marked with different letters differ from each other at P  More

A unified framework model of the H/P ratio
Based on the Lotka–Volterra equations24,25, the biomass dynamics of primary producers (P) and herbivores (H) are described as follows:

$$begin{array}{l}dP/dt = gleft( P right)P – xP – fleft( P right)PH,\ dH/dt = kfleft( P right)PH – mHend{array}$$
(1)

where g(P) is biomass-specific primary production (gC gC−1 d−1) and may be a function of P (gC m−2) owing to density-dependent growth, f(P) is per capita grazing rate of herbivores (m2 gC−1 d−1) and also may be function of P depending on the functional response, x is biomass-specific loss rate of primary producers other than due to grazing loss (gC gC−1 d−1), k is the conversion efficiency of herbivores as a fraction of ingested food converted into herbivore biomass (dimensionless: 0~1), and m is per capita mortality rate of herbivores owing to predation and other factors (gC gC−1 d−1). A list of model variables is listed in Supplementary Table 1. If we assume both g(P) and f(P) are constants, Eq. (1) is basically the Lotka–Volterra model, whereas it is an expansion of the Rosenzweig–MacArthur model if we assume logistic growth for g(P) and Michaelis–Menten (Holling type II) functional responses for f(P)24. At the equilibrium state, i.e., dP/dt = 0 and dH/dt = 0, the abundance of producers (P*) and consumers (H*) can be represented as:

$$H ast /P ast = left{ {left[ {gleft( {P ast } right) – x} right]} right./left. {left. {fleft( {P ast } right)} right]} right}/left{ {{mathrm{m/}}left[ {kfleft( {P ast } right)} right]} right}$$
(2)

Thus, the relationship between H and P is not affected by the types of the functional response in herbivores (f(P)). If we set g(P) as the biomass-specific primary production rate at equilibrium, as in simple Lotka–Volterra equations (i.e., g(P*) = g), then the H/P ratio can be expressed with log transformation as:

$${mathrm{log}}left( {H ast /P ast } right) = {mathrm{log}}left( k right) + {mathrm{log}}(g – x) – {mathrm{log}}left( m right)$$
(3)

At equilibrium, (g − x)P* is the amount of primary production that herbivores consume per unit of time (f(P*)P*H*). Thus, if this amount is divided by primary production per unit of time (P*g), it corresponds to the fraction of primary production that herbivores consume (0~1). We define it as β (= 1 − x/g). Large β values imply that producers are efficiently grazed at the equilibrium state. Thus, β is a gauge of inefficiency in the producers’ defensive traits. Using these parameters, the H*/P* ratio can be expressed as:

$${mathrm{log}}left( {H ast /P ast } right) = {mathrm{log}}left( k right) + {mathrm{log}}left( beta right) + {mathrm{log}}left( g right) – {mathrm{log}}left( m right)$$
(4)

This equation implies that the H*/P* biomass ratio on a log scale is affected additively by the specific primary production rate (log(g)), the grazeable fraction of primary production (log(β)), the conversion efficiency (log(k)), and the mortality rate of herbivores (log(m)). According to this equation, communities with relatively low carnivore abundance would have a correspondingly low value of m and will exhibit high herbivore biomass relative to producer biomass (H*/P*), whereas those with low primary production (with low value of g) owing to, for example, low light supply will have a low H*/P* ratio. An increase in defended producers such as armored plants or a decrease in edible producers will decrease β by increasing the loss rate x owing to the cost of defense, and will result in a decreased H*/P* ratio. Finally, when the nutritional value of producers decreases, the conversion efficiency of herbivores (k) should be low, which in turn decreases the H*/P*biomass ratio.
The model for a test with plankton communities
To apply Eq. (4) to a natural community, some modifications are necessary. Here, we consider a plankton community composed of algae and zooplankton. A theory of ecological stoichiometry suggests that the carbon content of primary producers relative to their nutrient content such as nitrogen or phosphorus is an important property affecting growth efficiency in herbivores4. Supporting the theory, a number of studies have shown that growth rate in terms of carbon accumulation relative to ingestion rate strongly depends on the carbon contents of the food relative to nutrients26,27,28,29. Thus, k can be expressed as:

$$k = q_1 times a_{{mathrm{nut}}}^{varepsilon 1}$$
(5)

where αnut is carbon content relative to nutrient content of primary producers and q1 is the conversion factor adjusting to biomass units. In this study, we applied a power function with coefficient of ε1 as a first order approximation because effects of this factor on the H*/P* biomass ratio may not be proportionally related to plant nutrient content. For example, if ε1 is much smaller than zero, it means that negative effects of the carbon to phosphorus ratio of algal food on an herbivore’s k are more substantial when the carbon to phosphorus ratio is high compared with the case when the carbon to phosphorus ratio is low. However, if this factor does not affect the H*/P* biomass ratio, ε1 = 0 and k is constant.
As herbivorous zooplankton cannot efficiently graze on larger phytoplankton due to gape limitation30, the feeding efficiency of herbivores or the defense efficiency of the producers’ resistance traits, β, would be related to the fraction of edible algae in terms of size as follows:

$$beta = q_2 times a_{{mathrm{edi}}}^{varepsilon 2}$$
(6)

where αedi is a trait determining producer edibility, q2 is a factor for converting the traits to edible efficiency, and ε2 is how effective the trait is in defending against grazing. We expect ε2 = 0 if this factor does not matter in regulating the H*/P* biomass ratio but ε2  > 0 if it has a role. Similarly, g can be described as

$$g = q_3 times mu ^{varepsilon 3}$$
(7)

where μ is the specific growth rate of producers, q3 is a conversion factor, and ε3 is the effect of µ on growth rate. Again, we expect that ε3 ≠ 0 if g has a role in determining the H*/P* ratio. Finally, assuming a Holling type I functional response of carnivores, the mortality rate of herbivores, m, is expressed as:

$$m = q_4 times theta ^{varepsilon 4}$$
(8)

where θ is abundance of carnivores, q4 is specific predation rate, and ε4 is the effect size of carnivore abundance on m.
By inserting Eqs. (5–8) to Eq. (4), effects of factors on the H*/P* biomass ratio is formulated as:

$${!}{mathrm{log}}left( {H ast /P ast } right) {!}= varepsilon _1{mathrm{log}}(a_{{mathrm{nut}}}) + varepsilon _2{mathrm{log}}(a_{{mathrm{edi}}}) + varepsilon _3{mathrm{log}}(mu ) – varepsilon _4{mathrm{log}}(theta ) + gamma$$
(9)

where γ is log(q1) + log(q2) + log(q3) − log(q4). If differences in the H*/P* ratio among communities are regulated by growth rate (μ), edibility (αedi), and nutrient contents (αnut) of producers as well as by predation by carnivores (θ), we expect non-zero values for ε1–ε4. Thus, Eq. (9) can be used to evaluate the relative importance of the four hypothesized agents if all of them simultaneously affect the H*/P* ratio. Here we undertake this analysis using data for natural plankton communities in experimental ponds where primary production rate was manipulated with different abundance of carnivore fish.
Experimental test by plankton communities
The experiment was carried out at two ponds (pond ID 217 and 218) located at the Cornell University Experimental Ponds Facility in Ithaca, NY, USA during 4 June to 28 August 2016 (Fig. 1). Each pond has a 0.09 ha surface area (30 × 30 m) and is 1.5 m deep. To initiate the experiment, we equally divided each of the two ponds into four sections using vinyl-coated canvas curtains, and randomly assigned the four sections to either high-shade (64% shading), mid-shade (47% shading), low-shade (33% shading), or no-shade treatments (no shading). Shading in each treatment was made using opaque floating mats (6 m diameter; Solar-cell SunBlanket, Century Products, Inc., Georgia, USA)31. The floating mats were deployed silvered side up to reflect sunlight and blue side down to avoid pond heating. Sampling was performed biweekly for water chemistry and abundance of phytoplankton and zooplankton with measurements of vertical profiles of water temperature, dissolved oxygen (DO) concentration and photosynthetic active radiation (PAR).
Fig. 1: Ponds used in the experimental test.

Pond 217 (a) and Pond 218 (b) in the Cornell University Experimental Ponds Facility divided into four sections by vinyl-canvas curtains and partially shaded by floating mats to regulate primary production rate. Floating docks were placed at the center of the ponds for sampling.

Full size image

PAR in the water column was lower in the sections with larger shaded areas throughout the experiment (Fig. 2b). Water temperature varied from 18 to 25 °C during the experiment but showed no notable differences in mean values of the water columns or the vertical profiles among the four treatments of the two ponds regardless of the shading treatment (Supplementary Fig. 1a, Supplementary Fig. 2). In all treatments, pH values gradually decreased towards the end of experiment and were higher in pond 218 (Supplementary Fig. 1b). DO concentration varied among the treatments and between the ponds, but were within the range of 5–12 mg L−1 (Supplementary Fig. 1c).
Fig. 2: Relationships between zooplankton biomass and phytoplankton biomass, and between H/P mass ratio and photosynthetic active radiation.

Biplots of zooplankton biomass (μg C L−1) and phytoplankton biomass (mg C L−1) (a) and H/P mass ratio and photosynthetic active radiation (PAR, mol photon m−2 d−1) (b) in the water column during the experiment in no-shade (blue), low-shade (orange), mid-shade (red), and high-shade (gray) treatments in pond 217 (circles) and 218 (squares). In each panel, small symbols denote values at each sampling date, and large symbols denote the mean values among the sampling dates. Bars denote standard errors on the means (n = 7 sampling date in each section). Correlation coefficients (r) with p values between the mean values are inserted in each panel.

Full size image

Phytoplankton biomass (mg C L−1) correlated significantly with chlorophyll a (µg L−1) (r = 0.702, p  More

The concentration of signal molecule ComX increases by the square of bacterial density
In order to determine the dynamics of signal molecule (ComX) production, we have used experimental and mathematical modeling approaches. We quantified the ComX concentration over time in the spent medium of PS-216 (ΔcomP) with the biosensor strain BD2876 (for strain-description see Supplementary Table 1), which produces β-galactosidase in response to the exogenous addition of ComX34. The assay included proper controls and calibrations to assure the biosensor-derived ComX concentrations are accurate (for details see Materials and methods). We found that the ComX concentration correlated positively with population density of PS-216 (ΔcomP) and remained constant at 10 nM after entering the stationary phase (Fig. 2a). Importantly, the representation of ComX concentration versus cell density (OD650) (Fig. 2b) showed a non-linear trend between the two parameters. The experimental data were fitted by an allometric function:

$$SMleft( t right) = aNleft( t right)^b$$
(1)

Where SM(t) is a signal molecule (ComX) concentration in time, N(t) is bacterial cell density in time, expressed as optical density of the bacterial suspension (OD650). The fitting results for parameters a and b were (9.6 ± 0.6) nM a.u.−2.09 and 2.09 ± 0.10, respectively. The value of parameter a means that at OD650 = 1.0 a.u., which corresponds to the stationary growth phase in our experimental conditions and the bacterial density of 4 × 108 cells mL‒1, the ComX concentration is about 10 nM. In the early exponential growth phase the concentration was about 0.1 nM. Considering parameter b, the value obtained (2.09 ± 0.10) indicates that the ComX concentration increases by the square of bacterial density. This means that with increasing population density, the SM concentration (ComX) increases by the second power, while the amount of ComX per cell increases linearly. This relationship suggests that ComQXPA has an ultra-sensitive encoder module9, where signal molecule production is very sensitive to cell density. The same mathematical relationship can be obtained by assuming that SM production rate per cell corresponds to the product of a specific cell growth rate and a cell density (i.e., population growth rate, for details, see Supplementary Methods, Derivation of ComQXPA communication system model). The dependence of the SM production rate per cell on (a) cell density and (b) the specific cell growth rate can be seen as an alternative way to obtain the ultra-sensitivity of encoders, which is usually achieved by SM dependent positive feedback in many QS systems9. This makes ComX a true indicator of population density, which also encodes information about the cell growth rate.
Fig. 2: The accumulation of signal molecule (SM) during the growth of B. subtilis and fitness cost of SM production.

a The growth curve (OD650) of B. subtilis PS-216 ΔcomP (no signal receptor) producing SM (ComX) that is accumulating in the growth medium of fermenter working in the batch mode of n = 3 biologically independent replicates is presented; b The experimental data n = 3 biologically independent replicates where the data ≥ limit of detection of SM was fitted by Eq. (1); the error bars for SM concertation are standard errors calculated from 8 technical replicates for each biological replicate; c The comparison of growth curves of B. subtilis PS-216 with no signal molecule receptor (ΔcomP) and no signal molecule production and receptor (ΔcomQXP) of n = 3 biologically independent replicates; the OD650 at t = 0 h was corrected with respect to the measured CFU of the inoculum. The slopes of the fitted lines in c correspond to the growth rate divided by log 2; the exponential phase points in the most reliable OD650 region ( >0.1 a.u. and  1.75 h, P  1 values obtained for all fits indicate positive cooperativity (i.e., ultrasensitivity14) in the binding of the transcriptional factor ComA to the srfA promoter43,44. This agrees with the research showing that two molecules of the ComA homodimer cooperatively bind to the two promoter regions located upstream of the RNAP binding sites of srfA13,20,45. The inactivation of the second promoter region decreases the promoter activity of srfA by 100-fold (ref. 13), which underscores the importance of the second binding region, explains n ≥ 2 and the sharp increase in srfA promoter activity with ComX concentration. In addition, we show here that the critical concentration of ComX required to induce quantifiable response (designated here as lower limit response (LLR)) is 0.2–0.5 nM. These results, therefore, suggest that the response per cell depends cooperatively on the ComX concentration and that the cells respond to very low concentrations of ComX.
Fig. 4: Influence of oxygen on the presence of SM (ComX).

a Signal molecule deficient B. subtilis PS-216 (∆comQ, PsrfAA-yfp) was incubated in the presence of SM for 4 h and the maximum normalized response was determined from the activity of the srfA promoter. n = 4 biologically independent replicates were performed. Best, concatenated fit to the model in Eq. (2) is presented together with 95% confidence level. b the logistic fit to one of the three growth curves (n = 3, biologically independent replicates) measured as OD650 of the culture B. subtilis PS 216 (srfA-lacZ) producing signal molecule, SM that accumulated in the growth medium of batch system is shown. The response per cell data, obtained as the β-galactosidase activity of srfA promoter of B. subtilis PS-216 (srfA-lacZ) was fitted by Eq. (3), (R2  > 0.99). The time interval at which SM concentration is high enough to cause the measurable response, i.e., lower limit of response (LLR) as predicted from data in experiments in (a) is given as dashed window in (b). One of the five qualitatively and quantitatively similar experimental results is presented. Error bars represent SD of 8 technical replicates. For fits of additional replicates and data variability refer to Supplementary Tables 3 and 4.

Full size image

Fully functional ComQXPA communication system does not require a positive feedback loop–the validation of the ComQXPA communication system
The response curve in Fig. 4a is a function of the SM concentration only. In the more natural setting (i.e., during growth) the cells encounter growth-dependent changes in SM concentrations as well as changes in bacterial density and growth rate over time. We have therefore asked whether the response curve based on the modeling and results presented in Figs. 2b and  4a could fit the response data in the SM producing and responding strain exposed to changes in these three parameters.
We cultivated the SM producing and responding PS-216 strain carrying the response reporter (PsrfA-lacZ) in a large volume bioreactor system (Fig. 4b). This allowed sterile sampling of spent medium and cells (for response quantification) at several time points, without affecting growth conditions. Immediately after the inoculation of the fresh medium by overnight culture the β-galactosidase activity of PS-216 (PsrfA-lacZ) was high. We assumed that this was a consequence of the accumulation of the expressed PsrfA reporter (β-galactosidase, RM) during the overnight growth. As a consequence of the dilution of the intracellular β-galactosidase (RM) due to cell division, the activity of the β-galactosidase decreased sharply after 2 h incubation (Fig. 4b). Simultaneously, as predicted by (Eq. 1), the concentration of SM (ComX) in the medium was increased exponentially during growth, and soon reached a critical concentration to activate the srfA promoter. In particular, as elucidated by the fits of (Eq. 2) to the data in Fig. 4a, the lower limit of the response (LLR) is reached shortly before upturn of the cell response curve in Fig. 4b. At this point the culture is in exponential growth phase at the cell density of 3 to 8 × 107 cells mL−1. The steep slope of the response curve indicates that the rate by which the response molecule (RM) is synthesized now exceeds the dilution due to the cell division rate. From now on, the response per cell correlates approximately linearly with OD650, suggesting a strong coupling to cell growth. Taking these facts into account and considering that the response molecule (RM) concentration is sensitive to the concentration of the signal, SM, (Eq. 2) and that SM can be expressed in terms of cell density (Eq. 1), the concentration of a response molecule per cell, RM(t)/N(t), can be analytically described (see also Supplementary Methods, Derivation of ComQXPA communication system model) as:

$$frac{{RM(t)}}{{N(t)}} = frac{{RM0}}{{N(t)}} + frac{{RM1(t)}}{{N(t)}}$$
(3)

where RM0/N(t) is the response per cell of overnight culture, i.e., the overnight accumulated β-galactosidase. The second term, RM1(t)/N(t) accounts for the synthesis of the β-galactosidase after inoculation of a fresh medium and comprises the parameters describing the sensitivity of the response to a signal molecule, Wmax, Km, n (Eq. 2), the signal production, a (Eq. 1), cell density, N(t) and proportionality constant, k that gives the magnitude of the response per cell when the potential to respond to the signal is maximally fulfilled (i.e., at Wmax) and the specific growth rate is 1 h−1. The definition of RM1(t) is given in Supplementary Methods (eq S11). Note that for the derivation of Eq. (3) we assumed no degradation of SM occurs, as our experiments suggest SM was stable under the conditions studied (Supplementary Figs. 6a and 7, see also Supplementary Methods, Derivation of ComQXPA communication). All the parameters in (Eq. 3), except k in RM1(t) and RM0, were taken as constants obtained in the independent experiments by fits of (Eq. 1) and (Eq. 2). With k and RM0, as the only fitting parameters, we applied the mathematical model in (Eq. 3) (for details of the model equation see Supplementary Methods, Derivation of ComQXPA communication system model) to fit the experimental cell response data (Fig. 4b). The successful fit (see Supplementary Table 4 for details) indicates that the relationship assumed among cell density, cell growth, signal concentration and response in (Eq. 3) is valid and yields (760 ± 120) M.U. for k and (5.5 ± 1.5) M.U. a.u. for RM0. Again, we did not need to incorporate the SM feedback loop into our model, which is consistent with published results suggesting that this communication system lacks a feedback loop10,11,12,13.
The ComQXPA dependent signaling and response at the cellular level
So far, we have focused on the population averages, which is a traditional approach in studies on microbial communication systems17,46. We here report results on communication dynamics of B. subtilis at the single cell level using fluorescence-based molecular tools. This approach provides the means to track temporal changes in expression of genes involved in signal synthesis (signaling) and in response and thus provides the insight into a phenotypic heterogeneity within the population.
We used the double-labeled fluorescencent strain B. subtilis PS-216 (comQ-yfp, srfA-cfp), in which fluorescent reporters were fused to the comQ and srfA promoters. The two genes code for the ComX signal-processing protein and the communication-activated operon, respectively. Since comQ and comX share the same promoter and their genetic sequences often overlap15 expression level of comQ corresponds to the expression level of comX. The fluorescence of individual cells was observed under the microscope in different growth phases and quantitatively analyzed (Fig. 5).
Fig. 5: Single cell quantitative fluorescence microscopy of B. subtilis PS-216 (comQ-yfp, srfA-cfp).

The fluorescence microscopy images were taken periodically during incubation of B. subtilis PS-216 (comQ-yfp, srfA-cfp) in a batch fermenter by the YFP filter (a), CFP filter (b) or DIC (c). The example YFP and CFP images, taken after 3 h of incubation were pseudo-colored. The scale bar represents 10 µm. n = 3 biologically independent experiments were performed. d % of population of cells that are hyper-expressing comQ-yfp or srfA-cfp is depicted. Gene expression level determined by single cell fluorescence microscopy was measured as Na-fluoresceinate standard normalized mean fluorescence intensity per cells expressing comQ-yfp (e) and srfA-cfp (f); ON is the overnight culture. One of the three qualitatively similar cell distributions is shown in (g) and (h) for comQ-yfp and srfA-cfp, respectively; areas under the curves are the same in all time points. For additional replicate see Supplementary Fig. 8.

Full size image

The observed expression pattern for the signaling gene (comQ-yfp) follows lognormal distribution (Fig. 5g). A small number of outliers in the comQ expression (on average, 10x brighter than the majority) were easily detected in the qualitative image analysis (Fig. 5a). These hyperproducers were not present in the overnight culture and began to occur during exponential growth, after 1-hour incubation in fresh medium, (Fig. 5d). In general, hyperproducers accounted for about 0.1–1% of the population, and their frequency increased during the first 6 hours. These data suggest that bulk of the ComX is nevertheless produced by the majority of the population as expected for the true QS communication systems. The contribution of the srfA hyperproducers to the total surfactin production is even less pronounced since their occurrence did not exceed 0.1% of the population (Fig. 5b, d).
The most heterogeneous expression of the communication signaling gene (comQ-yfp) was observed in overnight culture, immediately after its transfer to the fresh medium (Fig. 5g, Supplementary Fig. 8a), but hyperproducers where not detectable at this time (Fig. 5d). Once the cells begun to divide, the distribution shifted to lower fluorescence intensities with a simultaneous decrease in heterogeneity, but from 3 to 4 h onwards single cell fluorescence gradually increased, along with an increase in population heterogeneity (Fig. 5e, g, Supplementary Fig. 8a). This suggests that the expression rate is now higher than the division rate (i.e., the production overpowers the dilution due to cell division). A similar pattern was observed in the communication response (srfA-cfp) (Fig. 5f, h, Supplementary Fig. 8b) with two major differences. The minimum level in comQ-Yfp fluorescence was reached 1 h later than srfA-Cfp fluorescence, which suggests the expression of ComX is in first hours low compared to the srfA expression. Nevertheless, the entire cell population, with the exception of hyperproducers, which represented only a fraction of the comQ/srfA expressing cells, followed unimodal lognormal distribution expression pattern. This suggests that the ComQXPA communication phenomenon in B. subtilis, at least under the conditions in our experiment, is not restricted to individuals and can be studied at the population level, i.e., the averages represent well the population.
The comQ-yfp and srfA-cfp expression co-localization analysis (Supplementary Table 5) reveals the correlation coefficient of about 0.5, which is significantly (P = 0.01) higher than in the overnight culture. The presence of the correlation suggests that on average, cells that produce the signal more intensively, also respond to signal more intensively, supporting the idea of self-sensing47. However, the correlation coefficient strength was only moderate, suggesting that sensing of the external signal (sensing-of others) still works as expected for a typical QS system.
The induced response in ComQXPA communication system is graded and almost switch-like
The perfect QS system does not produce a response until the threshold bacterial density is reached and then immediately switches to a full response. This minimum to maximum transition may be either a perfect switch or a graded induction. By combining the information from Figs. 2b and 4a in the form of eq S9 results in the normalized ComQXPA response curve as a function of bacterial density (Fig. 6a) that resembles a graded switch like induction (compare to Fig. 1a). The perfect switch like communication system is unrealistic, because it requires that all the cells are perfectly synchronized and immediately switch to maximal response, leaving no time for the adaptation to the signal stimuli. It is reasonable to expect that for the true quorum sensing system (QS) most of the response has to occur within the same generation of dividing cells (ngen  More

1.
USDAFS. Press Release: Survey Finds 18 Million Trees Died in California in 2018. https://www.fs.usda.gov/Internet/FSE_DOCUMENTS/FSEPRD609321.pdf (USDAFS, 2019).
2.
Griffin, D. & Anchukaitis, K. J. How unusual is the 2012-2014 California drought? Geophys. Res. Lett. 41, 9017–9023 (2014).
ADS  Article  Google Scholar 

3.
Robeson, S. M. Revisiting the recent California drought as an extreme value. Geophys. Res. Lett. 42, 6771–6779 (2015).
ADS  Article  Google Scholar 

4.
Asner, G. P. et al. Progressive forest canopy water loss during the 2012-2015 California drought. Proc. Natl Acad. Sci. USA 113, E249–E255 (2016).
CAS  PubMed  Article  PubMed Central  Google Scholar 

5.
Brodrick, P. G. & Asner, G. P. Remotely sensed predictors of conifer tree mortality during severe drought. Environ. Res. Lett. 12, 115013 (2017).
ADS  Article  Google Scholar 

6.
Fettig, C. J. in Managing Sierra Nevada Forests. PSW-GTR-237 Ch. 2 (USDA Forest Service, 2012).

7.
Kolb, T. E. et al. Observed and anticipated impacts of drought on forest insects and diseases in the United States. For. Ecol. Manag. 380, 321–334 (2016).
Article  Google Scholar 

8.
Waring, R. H. & Pitman, G. B. Modifying lodgepole pine stands to change susceptibility to mountain pine beetle attack. Ecology 66, 889–897 (1985).
Article  Google Scholar 

9.
Restaino, C. et al. Forest structure and climate mediate drought-induced tree mortality in forests of the Sierra Nevada, USA. Ecol. Appl. 0, e01902 (2019).
Article  Google Scholar 

10.
USDAFS. Press Release: Record 129 million dead trees in California. https://www.fs.usda.gov/Internet/FSE_DOCUMENTS/fseprd566303.pdf (USDAFS, 2017).

11.
Young, D. J. N. et al. Long-term climate and competition explain forest mortality patterns under extreme drought. Ecol. Lett. 20, 78–86 (2017).
PubMed  Article  PubMed Central  Google Scholar 

12.
Raffa, K. F. et al. Cross-scale drivers of natural disturbances prone to anthropogenic amplification: The dynamics of bark beetle eruptions. BioScience 58, 501–517 (2008).
Article  Google Scholar 

13.
Boone, C. K., Aukema, B. H., Bohlmann, J., Carroll, A. L. & Raffa, K. F. Efficacy of tree defense physiology varies with bark beetle population density: a basis for positive feedback in eruptive species. Can. J. Res. 41, 1174–1188 (2011).
Article  Google Scholar 

14.
Fettig, C. J., Mortenson, L. A., Bulaon, B. M. & Foulk, P. B. Tree mortality following drought in the central and southern Sierra Nevada, California, U.S. For. Ecol. Manag. 432, 164–178 (2019).
Article  Google Scholar 

15.
Stephenson, N. L., Das, A. J., Ampersee, N. J. & Bulaon, B. M. Which trees die during drought? The key role of insect host-tree selection. J. Ecol. 75, 2383–2401 (2019).
Article  Google Scholar 

16.
Senf, C., Campbell, E. M., Pflugmacher, D., Wulder, M. A. & Hostert, P. A multi-scale analysis of western spruce budworm outbreak dynamics. Landsc. Ecol. 32, 501–514 (2017).
Article  Google Scholar 

17.
Seidl, R. et al. Small beetle, large-scale drivers: how regional and landscape factors affect outbreaks of the European spruce bark beetle. J. Appl Ecol. 53, 530–540 (2016).
Article  Google Scholar 

18.
Fettig, C. J. in Insects and Diseases of Mediterranean Forest Systems (eds Lieutier, F. & Paine, T. D.) 499–528 (Springer International Publishing, 2016).

19.
Raffa, K. F. & Berryman, A. A. The role of host plant resistance in the colonization behavior and ecology of bark beetles (Coleoptera: Scolytidae). Ecol. Monogr. 53, 27–49 (1983).
Article  Google Scholar 

20.
Logan, J. A., White, P., Bentz, B. J. & Powell, J. A. Model analysis of spatial patterns in mountain pine beetle outbreaks. Theor. Popul. Biol. 53, 236–255 (1998).
CAS  PubMed  MATH  Article  PubMed Central  Google Scholar 

21.
Wallin, K. F. & Raffa, K. F. Feedback between individual host selection behavior and population dynamics in an eruptive herbivore. Ecol. Monogr. 74, 101–116 (2004).
Article  Google Scholar 

22.
Franceschi, V. R., Krokene, P., Christiansen, E. & Krekling, T. Anatomical and chemical defenses of conifer bark against bark beetles and other pests. N. Phytol. 167, 353–376 (2005).
CAS  Article  Google Scholar 

23.
Raffa, K. F., Grégoire, J.-C. & Staffan Lindgren, B. Natural History and Ecology of Bark Beetles 1–40 (Elsevier, 2015).

24.
Bentz, B. J. et al. Climate change and bark beetles of the western United States and Canada: direct and indirect effects. BioScience 60, 602–613 (2010).
Article  Google Scholar 

25.
DeRose, R. J. & Long, J. N. Drought-driven disturbance history characterizes a southern Rocky Mountain subalpine forest. Can. J. Res. 42, 1649–1660 (2012).
Article  Google Scholar 

26.
Hart, S. J., Veblen, T. T., Schneider, D. & Molotch, N. P. Summer and winter drought drive the initiation and spread of spruce beetle outbreak. Ecology 98, 2698–2707 (2017).
PubMed  Article  PubMed Central  Google Scholar 

27.
Netherer, S., Panassiti, B., Pennerstorfer, J. & Matthews, B. Acute drought Is an important driver of bark beetle infestation in Austrian Norway spruce stands. Front. For. Glob. Change 2, 39 (2019).

28.
Kaiser, K. E., McGlynn, B. L. & Emanuel, R. E. Ecohydrology of an outbreak: mountain pine beetle impacts trees in drier landscape positions first. Ecohydrology 6, 444–454 (2013).
Article  Google Scholar 

29.
Marini, L. et al. Climate drivers of bark beetle outbreak dynamics in Norway spruce forests. Ecography 40, 1426–1435 (2017).
Article  Google Scholar 

30.
Sambaraju, K. R., Carroll, A. L. & Aukema, B. H. Multiyear weather anomalies associated with range shifts by the mountain pine beetle preceding large epidemics. For. Ecol. Manag. 438, 86–95 (2019).
Article  Google Scholar 

31.
Hayes, C. J., Fettig, C. J. & Merrill, L. D. Evaluation of multiple funnel traps and stand characteristics for estimating western pine beetle-caused tree mortality. J. Econ. Entomol. 102, 2170–2182 (2009).
PubMed  Article  PubMed Central  Google Scholar 

32.
Thistle, H. W. et al. Surrogate pheromone plumes in three forest trunk spaces: composite statistics and case studies. For. Sci. 50, 610–625 (2004).

33.
Miller, J. M. & Keen, F. P. Biology and Control of the Western Pine Beetle: A Summary of The First Fifty Years of Research (US Department of Agriculture, 1960).

34.
Chubaty, A. M., Roitberg, B. D. & Li, C. A dynamic host selection model for mountain pine beetle, Dendroctonus ponderosae Hopkins. Ecol. Model. 220, 1241–1250 (2009).
Article  Google Scholar 

35.
Graf, M., Reid, M. L., Aukema, B. H. & Lindgren, B. S. Association of tree diameter with body size and lipid content of mountain pine beetles. Can. Entomol. 144, 467–477 (2012).
Article  Google Scholar 

36.
Geiszler, D. R. & Gara, R. I. in Theory and Practice of Mountain Pine Beetle Management in Lodgepole Pine Forests: Symposium Proceedings (eds Berryman, A. A., Amman, G. D. & Stark, R. W.) (1978).

37.
Klein, W. H., Parker, D. L. & Jensen, C. E. Attack, emergence, and stand depletion trends of the mountain pine beetle in a lodgepole pine stand during an outbreak. Environ. Entomol. 7, 732–737 (1978).
Article  Google Scholar 

38.
Mitchell, R. G. & Preisler, H. K. Analysis of spatial patterns of lodgepole pine attacked by outbreak populations of the mountain pine beetle. For. Sci. 37, 1390–1408 (1991).
Google Scholar 

39.
Preisler, H. K. Modelling spatial patterns of trees attacked by bark-beetles. Appl. Stat. 42, 501 (1993).
MATH  Article  Google Scholar 

40.
Jactel, H. & Brockerhoff, E. G. Tree diversity reduces herbivory by forest insects. Ecol. Lett. 10, 835–848 (2007).
PubMed  Article  Google Scholar 

41.
Faccoli, M. & Bernardinelli, I. Composition and elevation of spruce forests affect susceptibility to bark beetle attacks: implications for forest management. Forests 5, 88–102 (2014).
Article  Google Scholar 

42.
Berryman, A. A. in Bark Beetles in North American Conifers: A System for the Study of Evolutionary Biology 264–314 (University of Texas Press, 1982).

43.
Fettig, C. J. et al. The effectiveness of vegetation management practices for prevention and control of bark beetle infestations in coniferous forests of the western and southern United States. For. Ecol. Manag. 238, 24–53 (2007).
Article  Google Scholar 

44.
Moeck, H. A., Wood, D. L. & Lindahl, K. Q. Host selection behavior of bark beetles (Coleoptera: Scolytidae) attacking Pinus ponderosa, with special emphasis on the western pine beetle, Dendroctonus brevicomis. J. Chem. Ecol. 7, 49–83 (1981).
CAS  PubMed  Article  Google Scholar 

45.
Evenden, M. L., Whitehouse, C. M. & Sykes, J. Factors influencing flight capacity of the mountain pine beetle (Coleoptera: Curculionidae: Scolytinae). Environ. Entomol. 43, 187–196 (2014).
CAS  PubMed  Article  Google Scholar 

46.
Raffa, K. F. & Berryman, A. A. Accumulation of monoterpenes and associated volatiles following inoculation of grand fir with a fungus transmitted by the fir engraver, Scolytus ventralis (Coleoptera: Scolytidae). Can. Entomol. 114, 797–810 (1982).
CAS  Article  Google Scholar 

47.
Anderegg, W. R. L. et al. Tree mortality from drought, insects, and their interactions in a changing climate. N. Phytol. 208, 674–683 (2015).
Article  Google Scholar 

48.
Kane, V. R. et al. Assessing fire effects on forest spatial structure using a fusion of Landsat and airborne LiDAR data in Yosemite National Park. Remote Sens. Environ. 151, 89–101 (2014).
ADS  Article  Google Scholar 

49.
Larson, A. J. & Churchill, D. Tree spatial patterns in fire-frequent forests of western North America, including mechanisms of pattern formation and implications for designing fuel reduction and restoration treatments. For. Ecol. Manag. 267, 74–92 (2012).
Article  Google Scholar 

50.
Morris, J. L. et al. Managing bark beetle impacts on ecosystems and society: Priority questions to motivate future research. J. Appl. Ecol. 54, 750–760 (2017).
Article  Google Scholar 

51.
Shiklomanov, A. N. et al. Enhancing global change experiments through integration of remote-sensing techniques. Front. Ecol. Environ. 17, 215–224 (2019).

52.
Jeronimo, S. M. A. et al. Forest structure and pattern vary by climate and landform across active-fire landscapes in the montane Sierra Nevada. For. Ecol. Manag. 437, 70–86 (2019).
Article  Google Scholar 

53.
Roussel, J.-R., Auty, D., De Boissieu, F. & Meador, A. S. lidR: Airborne LiDAR Data Manipulation and Visualization for Forestry Applications (2019).

54.
McDowell, N. et al. Mechanisms of plant survival and mortality during drought: why do some plants survive while others succumb to drought? N. Phytol. 178, 719–739 (2008).
Article  Google Scholar 

55.
Seybold, S. J. et al. Management of western North American bark beetles with semiochemicals. Annu. Rev. Entomol. 63, 407–432 (2018).
CAS  PubMed  Article  PubMed Central  Google Scholar 

56.
Fettig, C. J., McKelvey, S. R. & Huber, D. P. W. Nonhost angiosperm volatiles and verbenone disrupt response of western pine beetle, Dendroctonus brevicomis (Coleoptera: Scolytidae), to attractant-baited traps. J. Econ. Entomol. 98, 2041–2048 (2005).
CAS  PubMed  Article  PubMed Central  Google Scholar 

57.
Fettig, C. J., Dabney, C. P., McKelvey, S. R. & Huber, D. P. W. Nonhost angiosperm volatiles and verbenone protect individual ponderosa pines from attack by western pine beetle and red turpentine beetle (Coleoptera: Curculionidae, Scolytinae). West J. Appl. 23, 40–45 (2008).
Google Scholar 

58.
Fettig, C. J. et al. Efficacy of ‘Verbenone Plus’ for protecting ponderosa pine trees and stands from Dendroctonus brevicomis (Coleoptera: Curculionidae) attack in British Columbia and California. J. Econ. Entomol. 105, 1668–1680 (2012).
PubMed  Article  PubMed Central  Google Scholar 

59.
Oliver, W. W. Is self-thinning in ponderosa pine ruled by Dendroctonus bark beetles? In Forest Health Through Silviculture: Proceedings of the 1995 National Silviculture Workshop 6 (1995).

60.
Fettig, C. & McKelvey, S. Resiliency of an interior ponderosa pine forest to bark beetle infestations following fuel-reduction and forest-restoration treatments. Forests 5, 153–176 (2014).
Article  Google Scholar 

61.
Fettig, C. J. & Hilszczański, J. Bark Beetles 555–584. https://doi.org/10.1016/B978-0-12-417156-5.00014-9 (Elsevier, 2015).

62.
Chesson, P. Mechanisms of maintenance of species diversity. Annu. Rev. Ecol. Syst. 31, 343–366 (2000).
Article  Google Scholar 

63.
Fricker, G. A. et al. More than climate? Predictors of tree canopy height vary with scale in complex terrain, Sierra Nevada, CA (USA). For. Ecol. Manag. 434, 142–153 (2019).
Article  Google Scholar 

64.
Ma, S., Concilio, A., Oakley, B., North, M. & Chen, J. Spatial variability in microclimate in a mixed-conifer forest before and after thinning and burning treatments. For. Ecol. Manag. 259, 904–915 (2010).
Article  Google Scholar 

65.
Stovall, A. E. L., Shugart, H. & Yang, X. Tree height explains mortality risk during an intense drought. Nat. Commun. 10, 1–6 (2019).
Article  CAS  Google Scholar 

66.
Stephenson, N. L. & Das, A. J. Height-related changes in forest composition explain increasing tree mortality with height during an extreme drought. Nat. Commun. 11, 3402 (2020).
ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

67.
Stovall, A. E. L., Shugart, H. H. & Yang, X. Reply to ‘Height-related changes in forest composition explain increasing tree mortality with height during an extreme drought’. Nat. Commun. 11, 3401 (2020).
ADS  CAS  PubMed  PubMed Central  Article  Google Scholar 

68.
Person, H. L. Tree selection by the western pine beetle. J. For. 26, 564–578 (1928).
Google Scholar 

69.
Person, H. L. Theory in explanation of the selection of certain trees by the western pine beetle. J. For. 29, 696–699 (1931).
CAS  Google Scholar 

70.
Pile, L. S., Meyer, M. D., Rojas, R., Roe, O. & Smith, M. T. Drought impacts and compounding mortality on forest trees in the southern Sierra Nevada. Forests 10, 237 (2019).
Article  Google Scholar 

71.
Frey, J., Kovach, K., Stemmler, S. & Koch, B. UAV photogrammetry of forests as a vulnerable process. A sensitivity analysis for a structure from motion RGB-image pipeline. Remote Sens. 10, 912 (2018).
ADS  Article  Google Scholar 

72.
James, M. R. & Robson, S. Mitigating systematic error in topographic models derived from UAV and ground-based image networks. Earth Surf. Process. Landf. 39, 1413–1420 (2014).
ADS  Article  Google Scholar 

73.
Gray, P. C. et al. A convolutional neural network for detecting sea turtles in drone imagery. Methods Ecol. Evol. 10, 345–355 (2019).
Article  Google Scholar 

74.
Millar, C. I., Stephenson, N. L. & Stephens, S. L. Climate change and forests of the future: Managing in the face of uncertainty. Ecol. Appl. 17, 2145–2151 (2007).
PubMed  Article  PubMed Central  Google Scholar 

75.
Vose, J. M. et al. in Impacts, Risks, and Adaptation in the United States: The Fourth National Climate Assessment Vol II (eds Reidmiller, D. R., et al.) 232–267. https://nca2018.globalchange.gov/chapter/6/https://doi.org/10.7930/NCA4.2018.CH6 (2018).

76.
Bedard, W. D. et al. Western pine beetle: field response to its sex pheromone and a synergistic host terpene, myrcene. Science 164, 1284–1285 (1969).
ADS  CAS  PubMed  Article  PubMed Central  Google Scholar 

77.
Byers, J. A. & Wood, D. L. Interspecific inhibition of the response of the bark beetles, Dendroctonus brevicomis and Ips paraconfusus, to their pheromones in the field. J. Chem. Ecol. 6, 149–164 (1980).
CAS  Article  Google Scholar 

78.
Shepherd, W. P., Huber, D. P. W., Seybold, S. J. & Fettig, C. J. Antennal responses of the western pine beetle, Dendroctonus brevicomis (Coleoptera: Curculionidae), to stem volatiles of its primary host, Pinus ponderosa, and nine sympatric nonhost angiosperms and conifers. Chemoecology 17, 209–221 (2007).
CAS  Article  Google Scholar 

79.
DJI. Zenmuse X3 – Creativity Unleashed. DJI Official https://www.dji.com/zenmuse-x3/info (2015).

80.
Micasense. MicaSense. https://support.micasense.com/hc/en-us/articles/215261448-RedEdge-User-Manual-PDF-Download- (2015).

81.
DJI. DJI – The World Leader in Camera Drones/Quadcopters for Aerial Photography. DJI Official https://www.dji.com/matrice100/info (2015).

82.
Wyngaard, J. et al. Emergent challenges for science sUAS data management: fairness through community engagement and best practices development. Remote Sens. 11, 1797 (2019).
ADS  Article  Google Scholar 

83.
Rouse, W., Haas, R. H., Deering, W. & Schell, J. A. Monitoring the Vernal Advancement and Retrogradation (Green Wave Effect) of Natural Vegetation (Remote Sensing Center, Texas A&M Univ., 1973).

84.
DronesMadeEasy. Map Pilot for DJI on iOS. App Store https://itunes.apple.com/us/app/map-pilot-for-dji/id1014765000?mt=8 (2018).

85.
Farr, T. G. et al. The shuttle radar topography mission. Rev. Geophys. 45, RG2004 (2007).

86.
Zhang, W. et al. An easy-to-use airborne LiDAR data filtering method based on cloth simulation. Remote Sens. 8, 501 (2016).
ADS  Article  Google Scholar 

87.
Hijmans, R. J. et al. Raster: Geographic Data Analysis and Modeling (2019).

88.
Gitelson, A. & Merzlyak, M. N. Spectral reflectance changes associated with autumn senescence of Aesculus hippocastanum L. and Acer platanoides L. leaves. Spectral features and relation to chlorophyll estimation. J. Plant Physiol. 143, 286–292 (1994).
CAS  Article  Google Scholar 

89.
Coops, N. C., Johnson, M., Wulder, M. A. & White, J. C. Assessment of QuickBird high spatial resolution imagery to detect red attack damage due to mountain pine beetle infestation. Remote Sens. Environ. 103, 67–80 (2006).
ADS  Article  Google Scholar 

90.
Clevers, J. G. P. W. & Gitelson, A. A. Remote estimation of crop and grass chlorophyll and nitrogen content using red-edge bands on Sentinel-2 and -3. Int. J. Appl. Earth Observ. Geoinf. 23, 344–351 (2013).
ADS  Article  Google Scholar 

91.
Li, W., Guo, Q., Jakubowski, M. K. & Kelly, M. A new method for segmenting individual trees from the LiDAR point cloud. Photogramm. Eng. Remote Sens. 78, 75–84 (2012).
Article  Google Scholar 

92.
Jakubowski, M. K., Li, W., Guo, Q. & Kelly, M. Delineating individual trees from LiDAR data: a comparison of vector- and raster-based segmentation approaches. Remote Sens. 5, 4163–4186 (2013).
ADS  Article  Google Scholar 

93.
Shin, P., Sankey, T., Moore, M. & Thode, A. Evaluating unmanned aerial vehicle images for estimating forest canopy fuels in a ponderosa pine stand. Remote Sens. 10, 1266 (2018).
ADS  Article  Google Scholar 

94.
Roussel, J.-R. lidRplugins: Extra Functions and Algorithms for lidR Package (2019).

95.
Eysn, L. et al. A benchmark of LiDAR-based single tree detection methods using heterogeneous forest data from the alpine space. Forests 6, 1721–1747 (2015).
Article  Google Scholar 

96.
Vega, C. et al. PTrees: a point-based approach to forest tree extraction from LiDAR data. Int. J. Appl. Earth Observ. Geoinf. 33, 98–108 (2014).
ADS  Article  Google Scholar 

97.
Plowright, A. ForestTools: Analyzing Remotely Sensed Forest Data (2018).

98.
Pau, G., Fuchs, F., Sklyar, O., Boutros, M. & Huber, W. EBImage: An R package for image processing with applications to cellular phenotypes. Bioinformatics 26, 979–981 (2010).
CAS  PubMed  PubMed Central  Article  Google Scholar 

99.
Meyer, F. & Beucher, S. Morphological segmentation. J. Vis. Commun. Image Represent. 1, 21–46 (1990).
Article  Google Scholar 

100.
Hunziker, P. Velox: Fast Raster Manipulation and Extraction (2017).

101.
Kuhn, M. Building predictive models in R using the caret package. J. Stat. Softw. 28, 1–26 (2008).
Article  Google Scholar 

102.
Wang, Y. et al. Is field-measured tree height as reliable as believed – a comparison study of tree height estimates from field measurement, airborne laser scanning and terrestrial laser scanning in a boreal forest. ISPRS J. Photogramm. Remote Sens. 147, 132–145 (2019).
ADS  Article  Google Scholar 

103.
Weinstein, B. G., Marconi, S., Bohlman, S., Zare, A. & White, E. Individual tree-crown detection in RGB imagery using semi-supervised deep learning neural networks. Remote Sens. 11, 1309 (2019).
ADS  Article  Google Scholar 

104.
dos Santos, A. A. et al. Assessment of CNN-based methods for individual tree detection on images captured by RGB cameras attached to UAVs. Sensors (Basel) 19, 3595 (2019).

105.
Stephenson, N. Actual evapotranspiration and deficit: biologically meaningful correlates of vegetation distribution across spatial scales. J. Biogeogr. 25, 855–870 (1998).
Article  Google Scholar 

106.
Flint, L. E., Flint, A. L., Thorne, J. H. & Boynton, R. Fine-scale hydrologic modeling for regional landscape applications: The California Basin Characterization Model development and performance. Ecol. Process. 2, 25 (2013).
Article  Google Scholar 

107.
Millar, C. I. et al. Forest mortality in high-elevation whitebark pine (Pinus albicaulis) forests of eastern California, USA: influence of environmental context, bark beetles, climatic water deficit, and warming. Can. J. For. Res. 42, 749–765 (2012).
Article  Google Scholar 

108.
Baldwin, B. G. et al. Species richness and endemism in the native flora of California. Am. J. Bot. 104, 487–501 (2017).
PubMed  Article  PubMed Central  Google Scholar 

109.
Bürkner, P.-C. brms: an R package for bayesian multilevel models using Stan. J. Stat. Softw. 80, 1–28 (2017).
Article  Google Scholar 

110.
Hoffman, M. D. & Gelman, A. The no-U-turn sampler: adaptively setting path lengths in Hamiltonian Monte Carlo. J. Mach. Learn. Res. 15, 31 (2014).
MathSciNet  MATH  Google Scholar 

111.
Carpenter, B. et al. Stan: a probabilistic programming language. J. Stat. Softw. 76, 1–32 (2017).
Article  Google Scholar 

112.
Brooks, S. P. & Gelman, A. General methods for monitoring convergence of iterative simulations. J. Comput. Graph. Stat. 7, 434 (1998).
MathSciNet  Google Scholar 

113.
Gabry, J., Simpson, D., Vehtari, A., Betancourt, M. & Gelman, A. Visualization in Bayesian workflow. J. R. Stat. Soc. Ser. A 182, 389–402 (2019).
MathSciNet  Article  Google Scholar 

114.
Koontz, M. J., Latimer, A. M., Mortenson, L. A., Fettig, C. J. & North, M. P. Drone-derived data supporting “Cross-scale interaction of host tree size and climatic water deficit governs bark beetle-induced tree mortality”. https://doi.org/10.17605/OSF.IO/3CWF9 (2020).

115.
Baldwin, B. G. et al. Master spatial file for native California vascular plants used by Baldwin et al. (2017 Amer. J. Bot.), Dryad, Dataset, 2017. https://doi.org/10.6078/D16K5W.

116.
R Core Team. R: A Language and Environment for Statistical Computing (R Foundation for Statistical Computing, 2018).

117.
Koontz, M. J., Latimer, A. M., Mortenson, L. A., Fettig, C. J. & North, M. P. Local-structure-wpb-severity. https://doi.org/10.17605/OSF.IO/WPK5Z (2019). More

Study design
This study began in 2013, before P. destructans invaded Michigan and Wisconsin, and continued through 2020, after all Michigan and Wisconsin hibernacula were invaded by the pathogen. Throughout this period, we visited 22 hibernacula twice per winter during bat hibernation to quantify bat colony sizes, individual bat roosting temperatures, fungal loads (to which we added a constant 0.0001 before transforming to the log10 scale), and recapture probabilities (Fig. 2). We sampled bats during early hibernation in November, when more than 95% of swarm activity was expected to be finished51, and again during late hibernation in March, when less than 1% of spring emergence activity was expected to have begun51,52. For each sampling event, we counted all bats of all species within the site. We focused our analyses on the little brown bat (Myotis lucifugus). This species suffered severe declines due to WNS24,53, and was the only species abundant enough to provide sufficient sample sizes for our analyses. We divided bat counts by sections within the hibernaculum (i.e., “rooms”) that potentially varied in microclimate, and we used HOBO Pro v2 data loggers to continuously record temperatures every 3 h in these sections (Supplementary Figs. 2 and 7).
After counting all bats, we haphazardly sampled 20–25 individual little brown bats stratified across sections and roughly in proportion to the number of bats in each section (Fig. 2). For each sampled bat, we used a Fluke 62 MAX IR laser thermometer to quantify the temperature of the substrate directly adjacent to the bat (2 years is defined here as the post-invasion period) and that we had sampled during both the pre-invasion period (up to and including year 0, when WNS was first detected) and invasion period (years 1–2 of invasion). From pre-invasion to post-invasion, bat temperature distributions might have changed across these 12 hibernacula via two mechanisms: (1) bat preferences could have remained the same while populations declined in hibernacula and/or microsites with unfavorable temperatures, causing distribution shifts via selection; and/or (2) bat preferences could have changed towards avoidance of hibernacula and/or microsites with high temperature-mediated disease mortality, causing distribution shifts via selection or learning. In both cases, we would expect the average regional bat roosting temperatures to change over time. In contrast, we would only expect the average bat roosting temperature in any given site (or section) to change over time if the site provided a wide range of temperatures to select from, which was not always the case in the hibernacula we studied. Therefore, we quantified how the regional mean early hibernation roosting temperatures for all counted bats changed across all 12 of these hibernacula from pre-invasion to post-invasion (2013–2020, up to four years after the fungus was first detected) using a Gaussian regression, without including a random site effect. We visualized this change using density plots with 1 °C smoothing bin widths.
Though we stratified bat roosting temperature sampling by section, roughly in proportion to the total number of bats present in a section, we only sampled up to 25 little brown bats per site, and we could not sample exactly in proportion to the total bats present. Therefore, calculating the mean temperatures based on only sampled bats might be misleading. Instead, we used the observed temperatures from our sampled bats to estimate the roosting temperatures for bats that were counted but not sampled. In particular, we randomly generated an additional set of temperatures for all unmeasured, counted bats based on the mean and standard deviation of the roosting temperatures of sampled bats in the same section during the same survey. We then used both known (sampled bats) and estimated (counted bats) temperatures in our Gaussian regression to determine how average temperatures changed from pre-invasion to post-invasion. We also confirmed that qualitative results were the same if we only used the temperatures from sampled bats, without estimation. By considering how bats’ early hibernation temperature distributions changed throughout invasion, regardless of infection status, we were quantifying bat roosting preferences within and between hibernacula, rather than behavioral responses to infection within a season (i.e., “behavioral fever”).
Finally, we used the number of bats counted per hibernaculum during early hibernation (November) and late hibernation (March) in all hibernacula that still had bats in post-invasion years to calculate two demographic parameters: over summer population growth rates (reflecting immigration and recruitment to sites between March and November) and over winter mortality rates (population growth rates from November to March). We used a Gaussian linear model to determine whether the log10 population growth rates were correlated with average early hibernation roosting temperatures for sampled bats in each hibernaculum, the demographic season (over summer versus over winter), and the interaction between temperature and demographic season. A significant interaction term would indicate that temperature-correlated bat immigration and recruitment rates were decoupled from temperature-mediated population declines.
All analyses described above were performed in R version 3.5.161 with packages ‘ggplot2’62, ‘lme4’63, ‘loo’64, ‘pROC’65, and ‘R2jags56. Visual inspection of residuals and predictions plots confirmed acceptable model fits for all regression models.
Reporting summary
Further information on research design is available in the Nature Research Reporting Summary linked to this article. More

1.
Huang, S. D. et al. Autotrophic and heterotrophic soil respiration responds asymmetrically to drought in a subtropical forest in the Southeast China. Soil Biol. Biochem. 123, 242–249 (2018).
CAS  Article  Google Scholar 
2.
Peichl, M. et al. Management and climate effects on carbon dioxide and energy exchanges in a maritime grassland. Agric. Ecosyst. Environ. 158, 132–146 (2012).
Article  Google Scholar 

3.
Zhang, X. B. et al. How do environmental factors and different fertilizer strategies affect soil CO2 emission and carbon sequestration in the upland soils of southern China?. Appl. Soil Ecol. 72, 109–118 (2013).
Article  Google Scholar 

4.
Curiel, Y. J. et al. Microbial soil respiration and its dependency on carbon inputs, soil temperature and moisture. Glob. Change Biol. 13, 2018–2035 (2010).
ADS  Article  Google Scholar 

5.
Peri, P. L., Bahamonde, H. & Christiansen, R. Soil respiration in Patagonian semiarid grasslands under contrasting environmental and use conditions. J. Arid Environ. 119, 1–8 (2015).
ADS  Article  Google Scholar 

6.
Davidson, E. A. et al. Effects of soil water content on soil respiration in forests and cattle pastures of eastern Amazonia. Biogeochemistry 48, 53–69 (2000).
CAS  Article  Google Scholar 

7.
Buchmann, N. Biotic and abiotic factors controlling soil respiration rates in Piceaavies stands. Soil Biol Biochem. 32, 1625–1635 (2000).
CAS  Article  Google Scholar 

8.
Li, X. et al. Contribution of root respiration to total soil respiration in a semi-arid grassland on the Loess Plateau, China. Sci. Total Environ. 62, 1209–1217 (2018).
ADS  Google Scholar 

9.
Jia, B. & Zhou, G. Integrated diurnal soil respiration model during growing season of a typical temperate steppe: Effects of temperature, soil water content and biomass production. Soil Biol. Biochem. 41, 681–686 (2009).
CAS  Article  Google Scholar 

10.
Jin, Z., Qi, Y. C., Dong, Y. S. & Domroes, M. Seasonal patterns of soil respiration in three types of communities along grass-desert shrub transition in Inner Mongolia, China. Adv. Atmos. Sci. 26, 503–512 (2009).
CAS  Article  Google Scholar 

11.
Li, Y. C., Hou, C. C., Song, C. C. & Guo, Y. D. Seasonal changes in the contribution of root respiration to total soil respiration in a freshwater marsh in Sanjiang Plain, Northeast China. Environ. Earth Sci. 75, 1–10 (2016).
ADS  Article  Google Scholar 

12.
Zhao, C. Y., Zhao, Z. M., Yilihamu, Hong, Z. & Jun, L. Contribution of root and rhizosphere respiration of Haloxylonammodendron to seasonal variation of soil respiration in the central Asian desert. Quatern. Int. 244, 304–309 (2011).
Article  Google Scholar 

13.
Wang, X. et al. Soil respiration under climate warming: Differential response of heterotrophic and autotrophic respiration. Glob. Change Biol. 20, 3229–3237 (2014).
ADS  Article  Google Scholar 

14.
Sun, Z., Han, J. C. & Wang, H. Y. Soft rock for improving crop yield in sandy soil of Mu Us sandy land, China. Arid Land Res. Manag. 33, 136–154 (2019).
CAS  Article  Google Scholar 

15.
Dong, O. G. et al. water infiltration of covering soils with Different textures and bulk densities in gravelmulched areas. Appl. Soil Ecol. 17, 14039–14052 (2019).
Google Scholar 

16.
Lei, N., Han, J. C., Mu, X. M., Sun, Z. H. & Wang, H. Y. Effects of improved materials on reclamation of soil properties and crop yield in hollow villages in China. J. Soil Sediment 19, 2374–2380 (2019).
CAS  Article  Google Scholar 

17.
Fu, G. et al. Response of ecosystem respiration to experimental warming and clipping at daily time scale in an alpine meadow of Tibet. J. Mt. Sci.-Engl. 10, 455–463 (2013).
Article  Google Scholar 

18.
Wu, J. P. et al. Response of soil respiration and ecosystem carbon budget to vegetation removal in Eucalyptus plantations with contrasting ages. Sci. Rep.-UK 4, 6262 (2014).
Article  Google Scholar 

19.
Han, T. F., Huang, W. J., Liu, J. X., Zhou, G. Y. & Xiao, Y. Different soil respiration responses to litter manipulation in three subtropical successional forests. Sci. Rep.-UK 5, 18166 (2015).
ADS  CAS  Article  Google Scholar 

20.
Lei, N. & Han, J. C. Effect of precipitation on respiration of different reconstructed soils. Sci. Rep.-UK 10, 7328 (2020).
ADS  CAS  Article  Google Scholar 

21.
Kosugi, Y. et al. Spatial and temporal variation in soil respiration in a Southeast Asian tropical rainforest. Agric. Forest Meteorol. 147, 35–47 (2007).
ADS  Article  Google Scholar 

22.
Hu, Z. D., Liu, S. R., Shi, Z. M., Liu, X. L. & He, F. Diel variations and seasonal dynamics of soil respirations in subalpine meadow in western Sichuan Province, China. Acta Ecol. Sin. 32, 6376–6386 (2012).
Article  Google Scholar 

23.
Shi, Z. et al. Seasonal variation and temperature sensitivity of soil V respiration under different plant communities along an elevation gradient in Wuyi Mountains of China. J. Appl. Ecol. 19, 2357–2363 (2008).
Google Scholar 

24.
Wu, C. S., Sha, L. Q., Gao, J. M. & Dong, L. Y. Analysis of soil respiration in a montane evergreen broad-leaved forest and an artificial tea garden in Ailao Mountains, Yunnan province. J. Nanjing For. Univ. (Nat. Sci.) 36, 64–68 (2012) (in Chinese).
Google Scholar 

25.
Zhao, J. X. et al. Soil respiration and its affecting factors of Pinusyunnanensis in the Middle Regions of Yunnan. J. Northwest For. Univ. 30, 8–13 (2015) (in Chinese).
Google Scholar 

26.
Rayment, M. B. & Jarvis, P. G. Temporal and spatial variation of soil CO2 efflux in a Canadian boreal forest. Soil Biol. Biochem. 32, 35–45 (2000).
CAS  Article  Google Scholar 

27.
Qi, Y. & Xu, M. Separating the effects of moisture and temperature on soil CO2 efflux in a coniferous forest in the Sierra Nevada Mountains. Plant Soil 237, 15–23 (2001).
CAS  Article  Google Scholar 

28.
Dilustro, J. J., Collins, B., Duncan, L. & Crawford, C. Moisture and soil texture effects on soil CO2 efflux components in southeastern mixed pine forests. Forest Ecol. Manag. 204, 87–97 (2005).
Article  Google Scholar 

29.
Gaumont-guay, D. et al. Influence of temperature and drought on seasonal and interannual variations of soil, bole and ecosystem respiration in a boreal aspen stand. Agric Forest Meteorol. 140, 203–219 (2006).
ADS  Article  Google Scholar 

30.
Yan, J. X., Hao, Z., Jing, X. K. & Li, H. J. Variations of soil respiration and its temperature sensitivity in an orchard in Jinci Region of Taiyuan City. J. Environ. Sci.-China 37, 3625–3633 (2016).
Google Scholar 

31.
Han, G. X., Zhou, G. S. & Xu, Z. Z. Research and prospects for soil respiration of farmland ecosystems in China. J. Plant Ecol.-UK 32, 719–733 (2008).
CAS  Google Scholar 

32.
Yuste, J. C., Janssens, I. A., Carrara, A., Meiresonne, L. & Ceulemans, R. Interactive effects of temperature and precipitation on soil respiration in a temperate maritime pine forest. Tree Physiol. 23, 1263–1270 (2003).
Article  Google Scholar 

33.
Conant, R. T., Dalla-bett, A. P., Klopatek, C. C. & Klopatek, J. Controls on soil respiration in semiarid soils. Soil Biol. Biochem. 36, 945–951 (2004).
CAS  Article  Google Scholar 

34.
Giardina, C. P. & Ryan, M. G. Evidence that decomposition rates of organic carbon in mineral soil do not vary with temperature. Nature 404, 858–861 (2000).
ADS  CAS  PubMed  Article  Google Scholar 

35.
Grace, J. & Rayment, M. Respiration in the balance. Nature 404, 819–820 (2000).
CAS  PubMed  Article  Google Scholar 

36.
Andrews, R. T. et al. Successful embolization of collaterals from the ovarian artery during uterine artery embolization for fibroids: A case report. J. Vasc. Interv. Radiol. 11, 607–610 (2000).
CAS  PubMed  Article  Google Scholar 

37.
Fu, G., Shen, Z. X., Zhang, X. Z. & Zhou, Y. T. Response of soil microbial biomass to short-term experimental warming in alpine meadow on the Tibetan Plateau. Appl. Soil Ecol. 61, 158–160 (2012).
Article  Google Scholar 

38.
Chang, S. X., Shi, Z. & Thomas, B. R. Soil respiration and its temperature sensitivity in agricultural and afforested poplar plantation systems in northern Alberta. Biol. Fertil. Soils 52, 629–641 (2016).
CAS  Article  Google Scholar 

39.
Uribe, C. et al. Effect of wildfires on soil respiration in three typical Mediterranean forest ecosystems in Madrid, Spain. Plant Soil 369, 403–420 (2013).
CAS  Article  Google Scholar 

40.
Bergaust, L., Mao, Y. J., Bakken, L. R. & Frostegard, A. Denitrification response patterns during the transition to anoxic respiration and posttranscriptional effects of suboptimal pH on nitrogen oxide reductase in Paracoccusdenitrificans. Appl. Environ. Microbiol. 76, 6387–6396 (2010).
CAS  PubMed  PubMed Central  Article  Google Scholar 

41.
Mu, Z. et al. Linking N2O emission to soil mineral N as estimated by CO2 emission and soil C/N ratio. Soil Biol. Biochem. 41, 2593–2597 (2009).
CAS  Article  Google Scholar 

42.
Fu, G. et al. Response of soil respiration to grazing in an alpine meadow at three elevations in Tibet. Sci. World J. 2014, 1–9 (2014).
Google Scholar 

43.
Yu, C. Q., Han, F. S. & Fu, G. Effects of 7 years experimental warming on soil bacterial and fungal community structure in the Northern Tibet alpine meadow at three elevations. Sci. Total Environ. 655, 814–822 (2019).
ADS  CAS  PubMed  Article  PubMed Central  Google Scholar  More

Chemicals and reagents
Sodium chloride (S3014), potassium chloride (P9541), sodium phosphate dibasic (S3264), potassium phosphate monobasic (P9791), potassium bicarbonate (12602), ammonium chloride (A9434), ethylene diamine tetra acetic acid di-sodium salt (E6635), bovine serum albumin (BSA) (5482), absolute ethanol (100983), methanol (34860) were purchased from Sigma–Aldrich, St. Louis, MO, United States. CD 34, Sca 1, and CBA Flex Set which contains IL-3 (558346), IL-6 (558301), TNFα (558299), IFN γ (562233), G-CSF (560152), GM-CSF (558347), IL-1α (560157), IL-1β (560232), Mouse/Rat Soluble Protein Master Buffer Kit (558266) were procured from BD Biosciences, United States. May-Grunwalds Stain (S039) and Giemsa Stain Solution (TCL083) were procured from Hi-Media India.
1 L 1X PBS (8 g of sodium chloride, 0.2 g of potassium chloride, 1.44 g of sodium phosphate dibasic, 0.25 g of potassium phosphate monobasic to 1 L at pH 7.4), 1 L 1X RBC lysis buffer (1 g of potassium bicarbonate, 8 g of ammonium chloride and 0.03 g of di-sodium EDTA), 1% BSA in PBS, 70% ethanol in PBS, and may-grunwald-giemsa stain mixed in 3:1 ratio were prepared in the laboratory.
Animals and groupings
Pathogen free Strain ‘A’ male mice weighing 25.0–30.0 gm were received from Institute of Nuclear Medicine and Allied Science (INMAS) experimental animal facility at the age of 9 weeks, which is equivalent to young adult humans. Six mice per cage were housed under 25 ± 3 °C temperature and relative humidity of 30–70% in 12 h light/dark cycle with standard food and water. The study design strictly adhered to the guidelines approved by Institutional Animal Care and Use Committee (IACUC) of our institute, Institute of Nuclear Medicine and Allied Sciences (INMAS), Defence Research and Development Organization (DRDO), Delhi, India (Institute Animals Ethics Committee number: INM/IAEC/2019/01).
After 1 week of acclimatization under ambient conditions, mice were grouped for the survival assay according to localized gamma radiation dose i.e. 7, 9, 10, 12, 15, 17 and 19 Gy. Eighteen animals were assigned in each group for this assay.
For all other parameters mice were divided into two groups with six animals in each. Group 1, consisted of untreated or sham irradiated animals. Group 2 was exposed to 15 Gy localized irradiation. Mice were sacrificed at different time points (1, 2, 4, 7, 10, 15, 20, 25, and 30 days) post treatment. Hematopoietic stem cell marker CD 34 and Sca1 were measured on 1, 4, 7, and 10th day post-irradiation.
Irradiation
A Cobalt-60 teletherapy (Model: Bhabhatron-II, Panacea Medical Technologies Pvt Ltd, India) machine was used to irradiate the hinds of the three mice at a time (Fig. 1). The cobalt-60 beam was calibrated following IAEA TRS-398 protocol, with measurement done in water phantom (30 cm × 30 cm × 30 cm), for a field size 10 cm × 10 cm at source to surface distance (SSD) of 80 cm. A Farmer type, 0.6 cc volume, ionization chamber was utilised and the same was recently calibrated for absorbed dose to water (NDW) at Secondary Standard Dosimetry Laboratory (SSDL) of India. The measured absorbed dose rate (water, 10 cm × 10 cm, 80 SSD, Dmax) was 1.39 Gy/min.
Figure 1

Representative image of animal exposure with γ-rays in a field size of 2 cm × 35 cm.

Full size image

The localized bone marrow irradiation was done using single Posteroanterior (PA) beam with 2 cm × 35 cm field size at 80 cm SSD opened through secondary collimators of the machine. The 5 mm thick acrylic re-strainers were used to immobilize each mouse which also provides the build-up thickness for the irradiation. The dose prescription point was the surface of the hinds i.e. 5 mm below the surface of the re-strainer. The whole irradiation platform was placed on 5 cm thick acrylic slab to provide sufficient back scatter and all the air gaps between the three mice were filled with tissue equivalent gel bolus or wet cotton to compensate for missing side scatter. The prescription doses were absorbed dose to water and not corrected for any tissue heterogeneity.
The absolute and relative dosimetry of collimated field (2 cm × 35 cm) were performed using pin point (0.04 cc volume) ionization chamber and EBT3 (ISP, Wayne, USA) Gafchromic dosimetry films in solid water phantom. The response of the pin point chamber and EBT3 films was calibrated against the SSDL calibrated ionization chamber in solid water phantom at the depth of 5 cm and field size 10 cm × 10 cm in the same beam quality (i.e. cobalt-60). The film scanning was performed using Epson 10000XL flatbed scanner as per the manufacturers scanning protocol. The beam profiles were obtained from the scanned films and the full width at half maximum (FWHM) in transverse direction was 1.78 cm for the opening of 2.0 cm. The dose uniformity was better than 4% in central 80% of the field. The absolute dose rate (dose to water, dmax, 2 cm × 35 cm, SSD 80 cm) measured at the centre of the field was 1 Gy/min. The overall uncertainty in dose calculations is 6.5% with enhancement or coverage factor (k = 2) at 95% confidence. The main body of the mice was covered with 15 mm lead to prevent any scatter radiation reaching at other body parts of the animal. The surface dose reaching at the shielded body parts of the mice was measured using 0.4 cc ionisation chamber. The measured mean surface dose was 0.024 ± 0.014% of the open field dose. This dose is negligible and considered to be acceptable.
Survival assay in mice
For survival assay, 18 mice from each experimental group i.e. 7, 9, 10, 12, 15, 17 and 19 Gy, were irradiated at a fixed dose rate of 1 Gy/min. Groups were inspected daily for their morbidity and mortality status for a total duration of 30 days. All the surviving mice were euthanized humanely using intra-peritoneal injection (i.p.) with dexmedetomidine (0.3 mg/kg) in combination with ketamine (190 mg/kg) at the completion of the observation period.
Validation of hematopoietic syndrome model
The model of hematopoietic syndrome was validated with localized radiation of 15 Gy as this was the maximum dose at which all the animals survived as observed in the survival results. Mice were euthanized using intra-peritoneal injection (i.p.) with dexmedetomidine (0.3 mg/kg) in combination with ketamine (190 mg/kg) at 1, 2, 4, 7, 10, 15, 20, 25, and 30th day time points. Further analysis such as organ weight, cell counting, and serum cytokines measurement assays were conducted on various organs such as femur bone, spleen, thymus, lungs, kidneys (L), kidneys (R), liver and heart and blood at above mentioned time points. Hematopoietic stem cell markers (CD 34 and Sca 1) were measured in bone marrow cells of mice at 1, 4, 7, and 10th day time points.
Body and organ weights
The control and 15 Gy locally irradiated group animals were sacrificed at 1, 4, 7, 10, 15, 20, 25, and 30th day time points. The organs (spleen, thymus, lungs, kidneys (L), kidneys (R), liver and heart) were dissected out and weighed individually. Relative organ weight was calculated as the ratio between organ weight and body weight.
Haematology
The blood samples were collected from untreated and 15 Gy localized radiation group at 1, 4, 7, 10, 15, 20, 25, and 30th day time points in EDTA vial by cardiac puncture. 20 µl of collected blood from each experimental animal was taken for haematology and analysed using Nihon Kohden Celltac α, Tokyo, Japan, a fully automatic 3 part haematology analyzer.
Bone marrow cell counting
Total bone marrow cells were flushed out at 1, 4, 7, 10, 15, 20, 25 and 30th day points in 1 ml PBS from both femurs of the untreated mice and 15 Gy localized radiation group. The PBS containing cells was centrifuged at 2000 rpm for 8 min and the supernatant was discarded. The cells were washed with 1 ml PBS and the cell pellet was re-suspended in 1 ml of PBS. The re-suspended sample was diluted as required. The cells were then counted using improved Neubauer chamber under Olympus BX-63 microscope.
Bone marrow smears study
Bone marrow cells collected at 1, 4, 7, 10, 15, 20, 25, and 30th day post-irradiation from different treatment groups were centrifuged at 2000 rpm for 8 min and re-suspended in 50 μl of Fetal Bovine Serum (FBS). A small drop of cells suspension was dropped on a clean microscopic slide and then smeared thinly over an area of 2–3 cm by pulling another glass slide held at a 45° angle. Cells were fixed with methanol and stained with may-grunwald-giemsa. Slides were then analysed using Olympus BX-63 microscope.
Splenocytes and thymocytes counts
The organs excised at 1, 4, 7, 10, 15, 20, 25, and 30th day time-points were cleaned in chilled PBS. Using sterile-chilled frosted slides both the thymus and spleen were minced into a cell suspension and which was further centrifuged at 2000 rpm for 8 min. After centrifugation the supernatant was discarded. The RBCs were lysed using RBC lysis buffer and the cells were washed with 1 ml PBS. The cells were then counted using improved Neubauer chamber under Olympus BX-63 microscope.
Measurement of hematopoietic stem cell marker CD 34 and Sca1
Bone marrow cells were isolated from both the femurs from various treatment groups at different time points. Cell were washed with chilled PBS and fixed with 70% chilled ethanol and stored overnight at – 20 °C. Following day, five million cells were blocked with 1% BSA-PBS and stained with CD 34 and Sca 1 cell-surface markers for 20 min at room temperature. 10,000 cells from each sample were analyzed using FACS LSR II (BD Biosciences, San Jose, CA) and results were visualized by the BD FACSDiva V7 software (BD Biosciences, San Jose, CA).
Cytokines expression changes
Mice blood samples were collected in BD serum separation tube at different time-points post treatment (1, 4, 7, 10, 15, 25, and 30 days) via. cardiac puncture from different treatment groups. Samples were incubated at room temperature for 2 h and after incubation centrifuged at 6000g for 15 min. Serum was separated and stored at − 80 °C until analysis. The Levels of IL-3, IL-6, TNFα, IFN γ, G-CSF, GM-CSF, IL-1α,and IL-1β cytokines in mice serum were determined by using respective BD Cytometric Bead Array (CBA) Flex Set (BD Biosciences, USA) on dual laser flow analyzer (LSR-II, Becton Dickinson Biosciences, USA) according to manufacturer’s instructions. Results were analysed using FCAP Array V3 software (BD Biosciences, San Jose, CA).
Statistical analysis
All values obtained in our results are represented as mean ± SD of three independent replicates. The 30 day survival data was plotted using Kaplan–Meier analysis. The difference between the experimental groups was evaluated by one-way analysis of variance, with Newman–Keuls multiple comparison test (V, 8.01; GraphPad Prism, San Diego, CA, USA). The comparisons were made among the untreated and 15 Gy locally irradiated groups for all experimental parameters except the survival study. A value of p  More