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

1.
Russell JB. Discussion on diphtheria. BMJ. 1891;2:631–40.
Article  Google Scholar 
2.
Joseph FH. Notes on some pathogenic bacteria as found in the Transvaal, and the variations from their European prototype. Rep. S Afr Assoc Adv Sci. 1904;2:237–42.
Google Scholar 

3.
Eyre JW, Washbourn JW. Varities and virulence of the pneumococcus. Lancet. 1899;153:19–22.
Article  Google Scholar 

4.
Gladstone RA, Lo SW, Lees JA, Croucher NJ, van Tonder AJ, Corander J, et al. International genomic definition of pneumococcal lineages, to contextualise disease, antibiotic resistance and vaccine impact. EBioMedicine. 2019;43:338–46.
PubMed  PubMed Central  Article  Google Scholar 

5.
Colijn C, Corander J, Croucher NJ. Designing ecologically optimized pneumococcal vaccines using population genomics. Nat Microbiol. 2020;5:473–85.
CAS  PubMed  Article  PubMed Central  Google Scholar 

6.
Weinberger DM, Malley R, Lipsitch M. Serotype replacement in disease after pneumococcal vaccination. Lancet. 2011;378:1962–73.
PubMed  PubMed Central  Article  Google Scholar 

7.
Lefevre JC, Faucon G, Sicard AM, Gasc AM. DNA fingerprinting of Streptococcus pneumoniae strains by pulsed-field gel electrophoresis. J Clin Microbiol. 1993;31:2724 LP–8.
Article  Google Scholar 

8.
Selander RK, Levin BR. Genetic diversity and structure in Escherichia coli populations. Science. 1980;210:545–7.
CAS  PubMed  Article  PubMed Central  Google Scholar 

9.
Smith JM, Smith NH, O’Rourke M, Spratt BG. How clonal are bacteria? Proc Natl Acad Sci USA. 1993;90:4384–8.
CAS  PubMed  Article  PubMed Central  Google Scholar 

10.
Maiden MCJ, Bygraves JA, Feil E, Morelli G, Russell JE, Urwin R, et al. Multilocus sequence typing: a portable approach to the identification of clones within populations of pathogenic microorganisms. Proc Natl Acad Sci USA. 1998;95:3140–5.
CAS  PubMed  Article  PubMed Central  Google Scholar 

11.
Smith JM, Feil EJ, Smith NH. Population structure and evolutionary dynamics of pathogenic bacteria. BioEssays. 2000;22:1115–22.
CAS  PubMed  Article  PubMed Central  Google Scholar 

12.
Croucher NJ, Coupland PG, Stevenson AE, Callendrello A, Bentley SD, Hanage WP. Diversification of bacterial genome content through distinct mechanisms over different timescales. Nat Commun. 2014;5:5471.
PubMed  PubMed Central  Article  Google Scholar 

13.
Lees JA, Harris SR, Tonkin-Hill G, Gladstone RA, Lo S, Weiser JN, et al. Fast and flexible bacterial genomic epidemiology with PopPUNK. Genome Res. 2019;29:304–16.
CAS  PubMed  PubMed Central  Article  Google Scholar 

14.
Alikhan NF, Zhou Z, Sergeant MJ, Achtman M. A genomic overview of the population structure of Salmonella. PLoS Genet. 2018;14:e1007261.
PubMed  PubMed Central  Article  CAS  Google Scholar 

15.
Hanage WP, Spratt BG, Turner KME, Fraser C. Modelling bacterial speciation. Philos Trans R Soc B Biol Sci. 2006;361:2039–44.
Article  Google Scholar 

16.
Ford Doolittle W, Papke RT. Genomics and the bacterial species problem. Genome Biol. 2006;7:116.
PubMed  Article  CAS  PubMed Central  Google Scholar 

17.
Shapiro BJ, David LA, Friedman J, Alm EJ. Looking for Darwin’s footprints in the microbial world. Trends Microbiol. 2009;17:196–204.
CAS  PubMed  Article  PubMed Central  Google Scholar 

18.
Abu-Raddad LJ, Ferguson NM. The impact of cross-immunity, mutation and stochastic extinction on pathogen diversity. Proc R Soc B Biol Sci. 2004;271:2431–8.
Article  Google Scholar 

19.
Fraser C, Hanage WP, Spratt BG. Neutral microepidemic evolution of bacterial pathogens. Proc Natl Acad Sci USA. 2005;102:1968–73.
CAS  PubMed  Article  PubMed Central  Google Scholar 

20.
Corander J, Fraser C, Gutmann MU, Arnold B, Hanage WP, Bentley SD, et al. Frequency-dependent selection in vaccine-associated pneumococcal population dynamics. Nat Ecol Evol. 2017;1:1950–60.
PubMed  PubMed Central  Article  Google Scholar 

21.
Buckee CO, Jolley KA, Recker M, Penman B, Kriz P, Gupta S, et al. Role of selection in the emergence of lineages and the evolution of virulence in Neisseria meningitidis. Proc Natl Acad Sci USA. 2008;105:15082–7.
CAS  PubMed  Article  PubMed Central  Google Scholar 

22.
Croucher NJ, Finkelstein JA, Pelton SI, Mitchell PK, Lee GM, Parkhill J, et al. Population genomics of post-vaccine changes in pneumococcal epidemiology. Nat Genet. 2013;45:656–63.
CAS  PubMed  PubMed Central  Article  Google Scholar 

23.
Croucher NJ, Hanage WP, Harris SR, McGee L, van der Linden M, de Lencastre H, et al. Variable recombination dynamics during the emergence, transmission and ‘disarming’ of a multidrug-resistant pneumococcal clone. BMC Biol. 2014;12:49.
PubMed  PubMed Central  Article  CAS  Google Scholar 

24.
Croucher NJ, Chewapreecha C, Hanage WP, Harris SR, McGee L, van der Linden M, et al. Evidence for soft selective sweeps in the evolution of pneumococcal multidrug resistance and vaccine escape. Genome Biol Evol. 2014;6:1589–602.
PubMed  PubMed Central  Article  CAS  Google Scholar 

25.
Vos M. A species concept for bacteria based on adaptive divergence. Trends Microbiol. 2011;19:1–7.
CAS  PubMed  Article  PubMed Central  Google Scholar 

26.
Whitaker RJ. Allopatric origins of microbial species. Philos Trans R Soc B Biol Sci. 2006;361:1975–84.
Article  Google Scholar 

27.
Fraser C, Alm EJ, Polz MF, Spratt BG, Hanage WP. The bacterial species challenge: Making sense of genetic and ecological diversity. Science. 2009;323:741–6.
CAS  PubMed  Article  PubMed Central  Google Scholar 

28.
Hanage WP, Fraser C, Spratt BG. Fuzzy species among recombinogenic bacteria. BMC Biol. 2005;3:6.
PubMed  PubMed Central  Article  CAS  Google Scholar 

29.
Fraser C, Hanage WP, Spratt BG. Recombination and the nature of bacterial speciation. Science. 2007;315:476–80.
CAS  PubMed  PubMed Central  Article  Google Scholar 

30.
Sheppard SK, McCarthy ND, Falush D, Maiden MCJ. Convergence of Campylobacter species: implications for bacterial evolution. Science. 2008;320:237–9.
CAS  PubMed  Article  PubMed Central  Google Scholar 

31.
Laible G, Spratt BG, Hakenbeck R. Interspecies recombinational events during the evolution of altered PBP 2x genes in penicillin‐resistant clinical isolates of Streptococcus pneumoniae. Mol Microbiol. 1991;5:1993–2002.
CAS  PubMed  Article  PubMed Central  Google Scholar 

32.
Dowson CG, Coffey TJ, Kell C, Whiley RA. Evolution of penicillin resistance in Streptococcus pneumoniae; the role of Streptococcus mitis in the formation of a low affinity PBP2B in S. pneumoniae. Mol Microbiol. 1993;9:635–43.
CAS  PubMed  Article  PubMed Central  Google Scholar 

33.
Majewski J, Zawadzki P, Pickerill P, Cohan FM, Dowson CG. Barriers to genetic exchange between bacterial species: Streptococcus pneumoniae transformation. J Bacteriol. 2000;182:1016–23.
CAS  PubMed  PubMed Central  Article  Google Scholar 

34.
Bobay L-M, Ochman H. Biological Species Are Universal across Life’s Domains. Genome Biol Evol. 2017;9:491–501.
PubMed Central  Article  Google Scholar 

35.
Croucher NJ, Harris SR, Barquist L, Parkhill J, Bentley SD. A high-resolution view of genome-wide pneumococcal transformation. PLoS Pathog. 2012;8:e1002745.
CAS  PubMed  PubMed Central  Article  Google Scholar 

36.
Cohan F. Bacterial species and speciation. Syst Biol. 2001;50:513–24.
CAS  PubMed  Article  PubMed Central  Google Scholar 

37.
Nosil P, Funk DJ, Ortiz-Barrientos D. Divergent selection and heterogeneous genomic divergence. Mol Ecol. 2009;18:375–402.
PubMed  Article  PubMed Central  Google Scholar 

38.
Cohan FM, Perry EB. A systematics for discovering the fundamental units of bacterial diversity. Curr Biol. 2007;17:R373–R386.
CAS  PubMed  Article  PubMed Central  Google Scholar 

39.
Marttinen P, Hanage WP. Speciation trajectories in recombining bacterial species. PLoS Comput Biol. 2017;13:e1005640.
PubMed  PubMed Central  Article  CAS  Google Scholar 

40.
Gjini E, Valente C, Sá-Leão R, Gomes MGM. How direct competition shapes coexistence and vaccine effects in multi-strain pathogen systems. J Theor Biol. 2016;388:50–60.
PubMed  Article  PubMed Central  Google Scholar 

41.
Numminen E, Cheng L, Gyllenberg M, Corander J. Estimating the transmission dynamics of Streptococcus pneumoniae from strain prevalence data. Biometrics. 2013;69:748–57.
PubMed  Article  PubMed Central  Google Scholar 

42.
Majewski J, Cohan FM. Adapt globally, act locally: the effect of selective sweeps on bacterial sequence diversity. Genetics. 1999;152:1459–74.
CAS  PubMed  PubMed Central  Google Scholar 

43.
Wiedenbeck J, Cohan FM. Origins of bacterial diversity through horizontal genetic transfer and adaptation to new ecological niches. FEMS Microbiol Rev. 2011;35:957–76.
CAS  PubMed  Article  PubMed Central  Google Scholar 

44.
Barraclough TG, Balbi KJ, Ellis RJ. Evolving concepts of bacterial species. Evol Biol. 2012;39:148–57.
Article  Google Scholar 

45.
Chesson P. Mechanisms of maintenance of species diversity. Annu Rev Ecol Syst. 2000;31:343–66.
Article  Google Scholar 

46.
Neher RA, Shraiman BI. Competition between recombination and epistasis can cause a transition from allele to genotype selection. Proc Natl Acad Sci USA. 2009;106:6866–71.
CAS  PubMed  Article  PubMed Central  Google Scholar 

47.
Mostowy R, Croucher NJ, Andam CP, Corander J, Hanage WP, Marttinen P. Efficient inference of recent and ancestral recombination within bacterial populations. Mol Biol Evol. 2017;34:1167–82.
CAS  PubMed  PubMed Central  Article  Google Scholar 

48.
Croucher NJ, Harris SR, Fraser C, Quail MA, Burton J, van der Linden M, et al. Rapid pneumococcal evolution in response to clinical interventions. Science. 2011;331:430–4.
CAS  PubMed  PubMed Central  Article  Google Scholar 

49.
Romero P, Croucher NJ, Hiller NL, Hu FZ, Ehrlich GD, Bentley SD, et al. Comparative genomic analysis of ten Streptococcus pneumoniae temperate bacteriophages. J Bacteriol. 2009;191:4854–62.
CAS  PubMed  PubMed Central  Article  Google Scholar 

50.
Ferguson NM, Galvanl AP, Bush RM. Ecological and immunological determinants of influenza evolution. Nature. 2003;422:428–33.
CAS  PubMed  Article  PubMed Central  Google Scholar 

51.
Cobey S, Lipsitch M. Niche and neutral effects of acquired immunity permit coexistence of pneumococcal serotypes. Science. 2012;335:1376–80.
CAS  PubMed  PubMed Central  Article  Google Scholar 

52.
Cobey S. Pathogen evolution and the immunological niche. Ann N Y Acad Sci. 2014;1320:1–15.
CAS  PubMed  PubMed Central  Article  Google Scholar 

53.
Levin BR. Frequency-dependent selection in bacterial populations. Philos Trans R Soc Lond B Biol Sci. 1988;319:459–72.
CAS  PubMed  Article  PubMed Central  Google Scholar 

54.
Gomes MGM, Medley GF, Nokes DJ. On the determinants of population structure in antigenically diverse pathogens. Proc R Soc Lond Ser B Biol Sci. 2002;269:227–33.
Article  Google Scholar 

55.
Binsker U, Lees JA, Hammond AJ, Weiser JN. Immune exclusion by naturally acquired secretory IgA against pneumococcal pilus-1. J Clin Investig. 2020;130:927–41.
CAS  PubMed  Article  PubMed Central  Google Scholar 

56.
Gupta S, Maiden MCJ, Feavers IM, Nee S, May RM, Anderson RM. The maintenance of strain structure in populations of recombining infectious agents. Nat Med. 1996;2:437–42.
CAS  PubMed  Article  PubMed Central  Google Scholar 

57.
Gupta S, Ferguson N, Anderson R. Chaos, persistence, and evolution of strain structure in antigenically diverse infectious agents. Science. 1998;280:912–5.
CAS  PubMed  Article  PubMed Central  Google Scholar 

58.
Croucher NJ, Campo JJ, Le TQ, Liang X, Bentley SD, Hanage WP, et al. Diverse evolutionary patterns of pneumococcal antigens identified by pangenome-wide immunological screening. Proc Natl Acad Sci USA. 2017;114:E357–66.
CAS  PubMed  Article  PubMed Central  Google Scholar 

59.
Weinberger DM, Dagan R, Givon‐Lavi N, Regev‐Yochay G, Malley R, Lipsitch M. Epidemiologic evidence for serotype‐specific acquired immunity to pneumococcal carriage. J Infect Dis. 2008;197:1511–8.
PubMed  Article  PubMed Central  Google Scholar 

60.
Malley R, Lipsitch M, Bogaert D, Thompson CM, Hermans P, Claire Watkins A, et al. Serum antipneumococcal antibodies and pneumococcal colonization in adults with chronic obstructive pulmonary disease. J Infect Dis. 2007;196:928–35.
PubMed  Article  Google Scholar 

61.
Campo JJ, Le TQ, Pablo JV, Hung C, Teng AA, Tettelin H, et al. Panproteome-wide analysis of antibody responses to whole cell pneumococcal vaccination. Elife. 2018;7:e37015.
PubMed  PubMed Central  Article  Google Scholar 

62.
Buckee CO, Koelle K, Mustard MJ, Gupta S. The effects of host contact network structure on pathogen diversity and strain structure. Proc Natl Acad Sci USA. 2004;101:10839 LP–44.
Article  Google Scholar 

63.
Buckee CO, Recker M, Watkins ER, Gupta S. Role of stochastic processes in maintaining discrete strain structure in antigenically diverse pathogen populations. Proc Natl Acad Sci USA. 2011;108:15504 LP–9.
Article  Google Scholar 

64.
Watkins ER, Penman BS, Lourenço J, Buckee CO, Martin C, Maiden MCJ, et al. Vaccination drives changes in metabolic and virulence profiles of Streptococcus pneumoniae. PLoS Pathog. 2015;11:e1005034.
PubMed  PubMed Central  Article  CAS  Google Scholar 

65.
Skwark MJ, Croucher NJ, Puranen S, Chewapreecha C, Pesonen M, Xu YY, et al. Interacting networks of resistance, virulence and core machinery genes identified by genome-wide epistasis analysis. PLoS Genet. 2017;13:e1006508.
PubMed  PubMed Central  Article  CAS  Google Scholar 

66.
Pensar J, Puranen S, Arnold B, MacAlasdair N, Kuronen J, Tonkin-Hill G, et al. Genome-wide epistasis and co-selection study using mutual information. Nucleic Acids Res. 2019;47:e112.
CAS  PubMed  PubMed Central  Article  Google Scholar 

67.
Rocha EPC. Neutral theory, microbial practice: challenges in bacterial population genetics. Mol Biol Evol. 2018;35:1338–47.
CAS  PubMed  Article  PubMed Central  Google Scholar 

68.
Azarian T, Martinez PPP, Arnold BJ, Grant LR, Corander J, Fraser C, et al. Frequency-dependent selection can forecast evolution in Streptococcus pneumoniae. PLoS Biol. 2020;18:e3000878.
CAS  PubMed  PubMed Central  Article  Google Scholar 

69.
Johnston C, Campo N, Bergé MJ, Polard P, Claverys JP. Streptococcus pneumoniae, le transformiste. Trends Microbiol. 2014;22:113–9.
CAS  PubMed  Article  PubMed Central  Google Scholar 

70.
Marttinen P, Gutmann MU, Croucher NJ, Corander J, Hanage WP. Recombination produces coherent bacterial species clusters in both core and accessory genomes. Micro Genom. 2015;1:e000038.
Google Scholar 

71.
Croucher NJ, Mostowy R, Wymant C, Turner P, Bentley SD, Fraser C. Horizontal DNA Transfer mechanisms of bacteria as weapons of intragenomic conflict. PLoS Biol. 2016;14:e1002394.
PubMed  PubMed Central  Article  CAS  Google Scholar 

72.
Apagyi KJ, Fraser C, Croucher NJ. Transformation asymmetry and the evolution of the bacterial accessory genome. Mol Biol Evol. 2018;35:575–81.
CAS  PubMed  Article  PubMed Central  Google Scholar 

73.
Spearman C. The proof and measurement of association between two things. Am J Psychol. 1904;100:441–71.
Article  Google Scholar 

74.
Smirnov N. Table for estimating the goodness of fit of empirical distributions. Ann Math Stat. 1948;19:279–81.
Article  Google Scholar 

75.
Kelley T. Statistical method. New York: Macmillan; 1923.
Google Scholar 

76.
Simpson EH. Measurement of diversity. Nature. 1949;163:688.
Article  Google Scholar 

77.
Kimura M. On the probability of fixation of mutant genes in a population. Genetics. 1962;47:713–9.
CAS  PubMed  PubMed Central  Google Scholar 

78.
Hanage WP, Fraser C, Spratt BG. The impact of homologous recombination on the generation of diversity in bacteria. J Theor Biol. 2006;239:210–9.
CAS  PubMed  Article  PubMed Central  Google Scholar 

79.
Lawrence JG. Gene transfer in bacteria: speciation without species? Theor Popul Biol. 2002;61:449–60.
PubMed  Article  PubMed Central  Google Scholar 

80.
Pybus OG, Harvey PH. Testing macro-evolutionary models using incomplete molecular phylogenies. Proc R Soc B Biol Sci. 2000;267:2267–72.
CAS  Article  Google Scholar 

81.
Barraclough TG, Birky CW, Burt A. Diversification in sexual and asexual organisms. Evolution. 2003;57:2166–72.
PubMed  Article  PubMed Central  Google Scholar 

82.
Hudson RR. Gene genealogies and the coalescent process. In: Oxford surveys in evolutionary biology. Oxford: Oxford University Press; 1990. p. 1–44.

83.
Numminen E, Gutmann M, Shubin M, Marttinen P, Méric G, van Schaik W, et al. The impact of host metapopulation structure on the population genetics of colonizing bacteria. J Theor Biol. 2016;396:53–62.
PubMed  Article  PubMed Central  Google Scholar 

84.
Castillo-Chavez C, Hethcote HW, Andreasen V, Levin SA, Liu WM. Epidemiological models with age structure, proportionate mixing, and cross-immunity. J Math Biol. 1989;27:233–58.
CAS  PubMed  Article  PubMed Central  Google Scholar 

85.
McNally A, Kallonen T, Connor C, Abudahab K, Aanensen D, Horner C, et al. Signatures of negative frequency dependent selection in colonisation factors and the evolution of a multi-drug resistant lineage of Escherichia coli. MBio. 2019;10:e00644–19.
CAS  PubMed  PubMed Central  Article  Google Scholar 

86.
Kwun MJ, Oggioni MR, Bentley SD, Fraser C, Croucher NJ. Synergistic activity of mobile genetic element defences in Streptococcus pneumoniae. Genes. 2019;10:707.
CAS  PubMed Central  Article  Google Scholar 

87.
Iranzo J, Puigbo P, Lobkovsky AE, Wolf YI, Koonin EV. Inevitability of genetic parasites. Genome Biol Evol. 2016;8:2856–69. https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5631039/.
CAS  PubMed  PubMed Central  Article  Google Scholar 

88.
Budroni S, Siena E, Dunning Hotopp JC, Seib KL, Serruto D, Nofroni C, et al. Neisseria meningitidis is structured in clades associated with restriction modification systems that modulate homologous recombination. Proc Natl Acad Sci USA. 2011;108:4494–9.
CAS  PubMed  Article  PubMed Central  Google Scholar 

89.
Johnston C, Martin B, Granadel C, Polard P, Claverys JP. Programmed protection of foreign DNA from restriction allows pathogenicity island exchange during pneumococcal transformation. PLoS Pathog. 2013;9:e1003178.
CAS  PubMed  PubMed Central  Article  Google Scholar 

90.
Kwun MJ, Oggioni MR, De Ste Croix M, Bentley SD, Croucher NJ. Excision-reintegration at a pneumococcal phase-variable restriction-modification locus drives within- and between-strain epigenetic differentiation and inhibits gene acquisition. Nucleic Acids Res. 2018;46:11438–53.
CAS  PubMed  PubMed Central  Google Scholar  More

Experimental site
Seawater samples were collected from two stations (Fig. S1a, offshore OMZ station PS2 and coastal OMZ station PS3) in the ETNP in March and April 2018 on board R/V Sally Ride (Cruise ID: SR 1805). NO2− kinetics and O2 effects experiments (≤1 day incubation) were performed at the oxycline, the oxic/anoxic interface (top) of the ODZ and the core of the ODZ at each station (Fig. S1b and Table S1). Long incubations (6–7 days) were performed at the core of the ODZ at both stations to investigate nitrogen reaction rates at very low O2 concentrations.
Sampling and incubation experiments for NO2− kinetics and O2 effects
Twelve 30-L Niskin bottles on a rosette with a conductivity–temperature–depth (CTD) profiler were used to collect seawater while recording in situ O2 concentration, temperature, pressure, salinity, and chlorophyll fluorescence. O2 concentration at selected ODZ depths was measured by STOX sensor (on the CTD profiler) with detection limit of 10 nM [42]. NO2− concentrations were measured by standard spectrophotometric methods onboard. O2 and NO2− concentrations were used to select sampling depths. NO3− and NO2− concentrations in incubation samples were measured on a chemiluminescence NO/NOx Analyzer (Teledyne API, San Diego, CA, USA) in the lab.
Seawater was collected from Niskin bottles into 60 mL air-tight serum bottles after overflowing three times in order to minimize O2 contamination. Serum bottles were sealed with rubber septa and aluminum seals ensuring the absence of bubbles inside bottles. Septa were deoxygenated in anaerobic chambers with three cycles of vacuum/helium flushing over a period of 1 month prior to the cruise. A helium headspace was created for each sample collected from anoxic depths (the top and the core of the ODZ), and then samples were flushed with helium for at least 15 min to remove O2 that might have been introduced during sampling.
To determine a single NO2− oxidation rate, a set of five serum bottles amended with 15NO2− tracer was incubated on board at 12 °C in the dark. Incubations were terminated in time series (one bottle at day zero (T0), two at 0.5 day and two at 1 day) by adding 0.2 mL of saturated ZnCl2. The T0 bottles served as killed controls for both tracer contamination and abiotic reactions. The observed temporal changes in isotopic signals that occurred in the live samples over a few days in the time courses would not be detected if abiotic reactions were taking place during the >3 months that elapsed before all samples were measured on the mass spec.
In the lab, NO2− was removed from samples using sulfamic acid, and NO3− in serum bottles was converted into N2O using the denitrifier method [43, 44]. Both concentration and isotopic composition of N2O were measured on a mass spectrometer (Delta Vplus, Thermo Fisher Scientific, Waltham, MA, USA) for calculating nitrite oxidation rate from the linear regression of the five nitrate concentrations at the three time points as previously described [3]. NO2− kinetics of NO2− oxidation were determined by measuring NO2− oxidation rates under different 15NO2− tracer concentrations (0.5–13.8 μM). For responses of NO2− oxidation to O2, different amounts (0, 0.2, 0.5, 1, 2, and 5 mL) of O2 saturated seawater collected from the same Niskin bottle were added into serum bottles to achieve different O2 concentrations. O2 concentrations in serum bottles were monitored by optical oxygen sensors with a detection limit of 62.5 nM (PyroScience GmbH, Aachen, Germany).
Half-saturation constant
Half-saturation constant (Km) is the nitrite concentration at which nitrite oxidation rate (V) equals half of the potential maximum rate (Vm). The curve fitting tool in MATLAB_R2015a was used to fit the Michaelis–Menten model (Eq. (2)) to determine Km.

$${{V}} ,=, {{V}}_{mathrm{m}} ,times,left[ {{mathrm{NO}}_2^-} right]{mathrm{/}}left( {left[ {{mathrm{NO}}_2^-} right] ,+, {{K}}_{mathrm{m}}} right).$$
( 2)

Long incubation experiments
Longer (≥6 days) incubations were performed in 12-mL exetainers. Seawater from the core of the ODZ (250 m at station PS2, 160 m at station PS3) was sampled into 320 mL ground glass-stoppered bottles, which were immediately transferred into a N2 flushed glove bag. 15NO2− was added into these bottles to reach final concentrations of 7.24 and 8.01 µM at stations PS2 and PS3, respectively. 15NO2− labeled seawater was aliquoted into exetainers and capped within the glove bag. The septa had been stored under helium for at least 1 month. Exetainers were purged with helium for 5 min. Every 12 h, microbial activity in triplicate exetainers was terminated by adding 0.05 mL of 50% w/v ZnCl2. In the lab, N2 produced in exetainers was measured on a mass spectrometer (Europa Scientific 20–20, Crewe, UK) [45]. Denitrification and anammox rates were computed from linear regression of N2 produced at three time points. 15NO3− was measured (as described above) in the same exetainers, and NO2– oxidation rates were computed from linear regression of NO3− produced at three time points. O2 was monitored throughout the incubations in parallel exetainer vials using LUMOS sensors [12, 13]. Each O2 production or consumption rate was determined by linear regression of O2 concentrations at 32 time points. One sensing spot, glued inside an exetainer vial, allowed the optode to measure O2 concentrations every 5 min from the outside. Detection limit and resolution of LUMOS sensors was ≈1 nM.
DNA sampling, extraction, sequencing, metagenomics, and metatranscriptomics analysis
Particulate DNA samples were collected by filtration onto 0.22 µm Sterivex filters from the ODZ core samples (250 m at station PS2, 160 m at station PS3). DNA was extracted using the modified plant tissue protocol (All Prep DNA/RNA Mini Kit, Qiagen, Valencia, CA, USA), and subjected to paired-end sequencing on an Illumina MiSeq to generate over 10 million read pairs for each sample by the Genomics Core Facility of Lewis-Sigler Institute for Integrative Genomics at Princeton University. Quality control of raw reads was performed by BBDuk (DOE Joint Genome Institute, Walnut Creek, CA, USA), and assembled into contigs using metaSPAdes v3.12.0 [46] with specified options (-k 21, 33, 55, 77, 99, 127 -m 500). Metagenome-assembled genomes (MAGs) were constructed using MetaBAT v2.12.1 [47] with default (sensitive) mode and contigs longer than 1500 bp. The quality of MAGs was determined by checkM [48]. Taxonomy of MAGs was predicted by GTDB-tk [49]. ETNP PS2 MAG-11 from the ODZ core at station PS2 and ETNP PS3 MAG-54 from the ODZ core at station PS3 were identified as Nitrospina-like NOB. The taxonomy of these two NOB MAGs (ETNP PS2 MAG-11 and ETNP PS3 MAG-54) was further confirmed by comparing them with known OMZ NOB MAGs (MAG-1 and MAG-2) [8] using ANI. ANI between MAGs was assessed using enveomics [50]. MAG-2 [8], ETNP PS2 MAG-11 and ETNP PS3 MAG-54 belong to the same species (threshold for species: ANI ≥95%) based on their ANI values (Table S3). Considering the low completeness of MAG-11 and the high contamination of MAG-54 constructed from the two new metagenomes here (Table S4), we decided to use MAG-2 as the representative for this NOB species.
We estimated the relative abundance and the transcriptional activity of the two known OMZ NOB species in different oceanic regions (including the two ETNP stations in this study) by mapping metagenomes and metatranscriptomes from this and other studies to MAG-1 and MAG-2. The relative abundance of MAG-1 or MAG-2 at stations PS2 and PS3 in the ETNP was calculated as RPKG (the number of metagenomic reads obtained in this study mapped to a MAG per MAG length (kb) per genome equivalents) [15]. Genome equivalents were estimated using MicrobeCensus v1.1.1 [15]. The transcriptional activity of MAG-1 or MAG-2 in ETNP and ETSP OMZs was assessed by mapping published metatranscriptomic reads from the ETNP (PRJNA263621) [51] and the ETSP (SRA023632.1) [52] to MAG-1 and MAG-2. The relative abundance of RNA in Fig. S4 was calculated as the number of metatranscriptomic reads mapped to a MAG divided by the number of total reads. Mapping was performed by Bowtie2 [53] using “very-sensitive” mode, and only reads with a mapping quality above 20 were included as mapped reads.
In order to explore the possibility of the presence of other NOB in the ODZ core at stations PS2 and PS3, we also estimated the relative abundance of other (putative) NOB using their marker gene, nxrB (nitrite oxidoreductase). First, we downloaded previously identified (putative) nxrB sequences from the ETSP OMZ [8]. Then, we color coded the genes in Fig. 4 based on previously defined clusters in a phylogenetic tree (see Fig. 3 in [8]): the Nitrospina cluster (blue) contains nxrB grouped with cultured marine NOB, Nitrospina gracilis. Based on BLASTp search results, amino acid sequence identities between all the OMZ nxrB in this cluster and that of Nitrospina gracilis were 96.71%, 96.71%, and 96.87% for NODE_69234, NODE_114897, and NODE_102582, respectively. Only this cluster is associated with known NOB, and the nitrite oxidation capacity of all the other genes associated with nxrB needs to be confirmed. The anammox cluster (red) contains anammox nxrB sequences. The putative cluster (black) contains nxrB grouped with microbes in which nitrite oxidation capacity has not been proven. This putative cluster fell between known NOB and anammox, and was implied to represent unidentified NOB [54]. The last cluster is called unknown nxrB (gray) because neither their phylogeny nor function can be determined based on their distant relationship with known NOB or anammox nxrB. Finally, we mapped the ETNP metagenomic reads obtained in this study to each nxrB gene. Relative abundance of nxrB genes was also expressed as RPKG: relative abundance of nxrB gene = (number of mapped reads to a certain nxrB gene)/(length of this nxrB in kb)/(genome equivalents).
To explore the potential metabolisms of NOB in anoxic ODZ cores, we searched for chlorite dismutase and NO dismutase (nod) genes in the NOB MAGs. First, protein-coding sequences in the two new NOB MAGs obtained here (ETNP PS2 MAG-11 and ETNP PS3 MAG-54) were predicted by Prodigal v2.6.3 [55]. The protein-coding sequences were annotated by the best BLASTp hits against the nr protein database. DNA sequences encoding Cld were identified in both ETNP NOB MAGs. Predicted Cld amino acid sequences encoded by ETNP MAGs were too short (120 aa) to be compared to Cld of Candidatus Nitrospira defluvii, but they only had one mismatch with Cld amino acid sequences of MAG-2 from the ETSP OMZ based on MUSCLE alignment using MEGA 7 software (Fig. S5b). Thus, we looked for the arginine173, the Cld activity marker, in longer Cld sequences of MAG-1 and MAG-2 by aligning their Cld with the Cld of Candidatus Nitrospira defluvii using MUSCLE in MEGA 7 software. Since nod was not found in the ETNP MAGs via annotation and was not reported in MAG-1 and MAG-2 in the previous study [8], gene search (i.e., searching nod against NOB MAGs) was performed using Hidden Markov Models by HMMER3 [56]. Reference sequences of the search included a nod sequence retrieved from Candidatus Methylomirabilis oxyfera [32] genome (accession numbers FP565575.1), and three environmental nod sequences (accession numbers: KX364450.1, KX364454.1, and KU933965.1). Search queries were the two ETNP MAGs and two ETSP MAGs in Table S3.
Estimation of NO2− oxidation through disproportionation using an inverse isotope model
We simulated the distribution of NO2− oxidation rates via disproportionation using a 1-D inverse isotope model [37] in the anoxic ODZ core from two stations in the ETSP OMZ (Fig. S6), where the complete suite of isotope and rate data have been previously published (Table S7). Briefly, the net biochemical production or consumption rate of each nitrogen compound (R14Ammonium, R14Nitrite, R14Nitrate, R15Nitrite, and R15Nitrate) was balanced by vertical diffusion and advection at steady state, and then five equations were used to balance measured 14NH4+, 14NO2−, 15NO2−, 14NO3−, and 15NO3− concentrations (Eqs. (3–7)). The exclusion of horizontal processes was justified [37] since the model was not run all the way to the surface, and using a constant vertical advection term does not violate continuity of the model. F is the rate of each nitrogen cycle process represented in this model (Fig. 5a, b). We modified the previous model by replacing canonical NO2− oxidation with NO2− disproportionation and using the recently determined stoichiometry of nitrate production by anammox from 0.3 to 0.16 [24]. Since OMZ NOB genomes encode nitrite oxidoreductase (catalyzes NO2−→NO3−) and nitrite reductase (catalyzes NO2−→NO), we assumed that the fractionation factor of nitrate production from nitrite by disproportionation is the same as by canonical nitrite oxidation (αDis = αNxr) based on the close clustering between OMZ NOB nxr and aerobic NOB (Nitrospina gracilis) nxr [8], and the fractionation factor of N2 produced by disproportionation is the same as by denitrification (αDisN2 = αNir). Rates of NO2− disproportionation (FDis), NO3− reduction (FNar), denitrification (FNir), and anammox (FAmx) were solved from the equations by the nonnegative least squares optimization routine (lsqnonneg) in MATLAB_R2015a as described previously [37].

$$R_{14Ammonium} ,=, 0.11,times {,}^{14}F_{Nir} ,+, 0.07,times {,}^{14}F_{Nar} ,- {,}^{14}F_{Amx},$$
( 3)

$$R_{14Nitrite} = – ,^{14}F_{Nir} ,+, {,}^{14}F_{Nar} ,-, (1 ,+, c) ,times,{,}^{14}F_{Amx} \ -, left( {5/3} right) ,times,{,}^{14}F_{Dis},$$
( 4)

$$R_{14Nitrate} ,=, – ^{14}F_{Nar} ,+, c ,times {,}^{14}F_{Amx} ,+, {,}^{14}F_{Dis},$$
( 5)

$$R_{15Nitrite} = -!! {,}^{14}F_{Nir}/alpha_{Nir} ,times,left( {left[ {{,}^{15}{rm{NO}}_{2}^{-}} right]/left[ {{,}^{14}{rm{NO}}_{2}^{-}} right]} right) \ +, {,}^{14}F_{Nar}/alpha_{Nar},times,left( {left[ {{,}^{15}{rm{NO}}_{3}^{-}} right]/left[ {{,}^{14}{rm{NO}}_{3}^{-}} right]} right)\ -, {,}^{14}F_{Amx}/alpha_{Amx},times,left( {left[ {{,}^{15}{rm{NO}}_{2}^{-}} right]/left[ {{,}^{14}{rm{NO}}_{2}^{-}} right]} right) \ -, c times !{,}^{14}F_{Amx}/alpha_{NxrAmx},times,left( {left[ {{,}^{15}{rm{NO}}_{2}^{-}} right]/left[ {{,}^{14}{rm{NO}}_{2}^{-}} right]} right)\ -, {,}^{14}F_{Dis}/alpha_{Dis}timesleft( {left[ {{,}^{15}{rm{NO}}_{2}^{-}} right]/left[ {{,}^{14}{rm{NO}}_{2}^{-}} right]} right) \ -, left( {2/3} right) !times ! {,}^{14}F_{Dis}/alpha_{DisN2} !times left(! {left[ {{,}^{15}{rm{NO}}_{2}^{-}} right]/left[ !{{,}^{14}{rm{NO}}_{2}^{-}} right]} !right),$$
(6)

$$R_{15Nitrate} = -!! {,}^{14}F_{Nar}/alpha_{Nar} ,times,left( {left[ {{,}^{15}{rm{NO}}_{3}^{-}} right]/left[ {{,}^{14}{rm{NO}}_{3}^{-}} right]} right) \ +, c times !{,}^{14}F_{Amx}/alpha_{NxrAmx},times,left( {left[ {{,}^{15}{rm{NO}}_{2}^{-}} right]/left[ {{,}^{14}{rm{NO}}_{2}^{-}} right]} right)\ +, ^{14}F_{Dis}/alpha_{Dis},times,left( {left[ {{,}^{15}{rm{NO}}_{2}^{-}} right]/left[ {{,}^{14}{rm{NO}}_{2}^{-}} right]} right).$$
(7) More

1.
IPCC. Managing the risks of extreme events and disasters to advance climate change adaptation. A special report of working groups i and ii of the intergovernmental panel on climate change. Cambridge, UK, and New York, NY, USA: Cambridge University Press; 2012.
Google Scholar 
2.
Trenberth KE, Dai AG, van der Schrier G, Jones PD, Barichivich J, Briffa KR, et al. Global warming and changes in drought. Nat Clim Chang. 2014;4:17–22.
Article  Google Scholar 

3.
Kreyling J, Henry HAL. Vanishing winters in germany: soil frost dynamics and snow cover trends, and ecological implications. Clim Res. 2011;46:269–76.
Article  Google Scholar 

4.
Congreves KA, Wagner-Riddle C, Si BC, Clough TJ. Nitrous oxide emissions and biogeochemical responses to soil freezing-thawing and drying-wetting. Soil Biol Biochem. 2018;117:5–15.
CAS  Article  Google Scholar 

5.
Schimel J, Balser TC, Wallenstein M. Microbial stress-response physiology and its implications for ecosystem function. Ecology 2007;88:1386–94.
PubMed  Article  Google Scholar 

6.
Tecon R, Or D. Biophysical processes supporting the diversity of microbial life in soil. FEMS Microbiol Rev. 2017;41:599–623.
CAS  PubMed  PubMed Central  Article  Google Scholar 

7.
Koponen HT, Bååth E. Soil bacterial growth after a freezing/thawing event. Soil Biol Biochem. 2016;100:229–32.
CAS  Article  Google Scholar 

8.
Schostag M, Priemé A, Jacquiod S, Russel J, Ekelund F, Jacobsen CS. Bacterial and protozoan dynamics upon thawing and freezing of an active layer permafrost soil. ISME J. 2019;13:1345–59.
CAS  PubMed  PubMed Central  Article  Google Scholar 

9.
Barnard RL, Osborne CA, Firestone MK. Changing precipitation pattern alters soil microbial community response to wet-up under a mediterranean-type climate. ISME J. 2015;9:946–57.
CAS  PubMed  Article  Google Scholar 

10.
Drigo B, Nielsen UN, Jeffries TC, Curlevski NJA, Singh BK, Duursma RA, et al. Interactive effects of seasonal drought and elevated atmospheric carbon dioxide concentration on prokaryotic rhizosphere communities. Environ Microbiol. 2017;19:3175–85.
CAS  PubMed  Article  PubMed Central  Google Scholar 

11.
Priemé A, Christensen S. Natural perturbations, drying-wetting and freezing-thawing cycles, and the emission of nitrous oxide, carbon dioxide and methane from farmed organic soils. Soil Biol Biochem. 2001;33:2083–91.
Article  Google Scholar 

12.
Kim DG, Vargas R, Bond-Lamberty B, Turetsky MR. Effects of soil rewetting and thawing on soil gas fluxes: a review of current literature and suggestions for future research. Biogeosciences 2012;9:2459–83.
CAS  Article  Google Scholar 

13.
Meisner A, Leizeaga A, Rousk J, Bååth E. Partial drying accelerates bacterial growth recovery to rewetting. Soil Biol Biochem. 2017;112:269–76.
CAS  Article  Google Scholar 

14.
Blazewicz SJ, Schwartz E, Firestone MK. Growth and death of bacteria and fungi underlie rainfall-induced carbon dioxide pulses from seasonally dried soil. Ecology. 2014;95:1162–72.
PubMed  Article  PubMed Central  Google Scholar 

15.
Barnard RL, Osborne CA, Firestone MK. Responses of soil bacterial and fungal communities to extreme desiccation and rewetting. ISME J. 2013;7:2229–41.
CAS  PubMed  PubMed Central  Article  Google Scholar 

16.
Meisner A, Jacquiod S, Snoek BL, ten Hooven FC, van der Putten WH. Drought legacy effects on the composition of soil fungal and prokaryote communities. Front Microbiol. 2018;9:294.
PubMed  PubMed Central  Article  Google Scholar 

17.
de Vries FT, Griffiths RI, Bailey M, Craig H, Girlanda M, Gweon HS, et al. Soil bacterial networks are less stable under drought than fungal networks. Nat Commun. 2018;9:3033.
PubMed  PubMed Central  Article  CAS  Google Scholar 

18.
Bouskill NJ, Lim HC, Borglin S, Salve R, Wood TE, Silver WL, et al. Pre-exposure to drought increases the resistance of tropical forest soil bacterial communities to extended drought. ISME J. 2013;7:384–94.
CAS  PubMed  Article  Google Scholar 

19.
Ren J, Song C, Hou A, Song Y, Zhu X, Cagle GA. Shifts in soil bacterial and archaeal communities during freeze-thaw cycles in a seasonal frozen marsh, northeast china. Sci Total Environ. 2018;625:782–91.
CAS  PubMed  Article  Google Scholar 

20.
Mannisto MK, Tiirola M, Haggblom MM. Effect of freeze-thaw cycles on bacterial communities of arctic tundra soil. Micro Ecol. 2009;58:621–31.
Article  Google Scholar 

21.
Evans SE, Wallenstein MD. Climate change alters ecological strategies of soil bacteria. Ecol Lett. 2014;17:155–64.
PubMed  Article  PubMed Central  Google Scholar 

22.
Fuchslueger L, Bahn M, Hasibeder R, Kienzl S, Fritz K, Schmitt M, et al. Drought history affects grassland plant and microbial carbon turnover during and after a subsequent drought event. J Ecol. 2016;104:1453–65.
CAS  PubMed  PubMed Central  Article  Google Scholar 

23.
Hawkes CV, Waring BG, Rocca JD, Kivlin SN. Historical climate controls soil respiration responses to current soil moisture. Proc Natl Acad Sci. 2017;114:6322–7.
CAS  PubMed  Article  PubMed Central  Google Scholar 

24.
Preece C, Verbruggen E, Liu L, Weedon JT, Peñuelas J. Effects of past and current drought on the composition and diversity of soil microbial communities. Soil Biol Biochem. 2019;131:28–39.
CAS  Article  Google Scholar 

25.
de Nijs EA, Hicks LC, Leizeaga A, Tietema A, Rousk J. Soil microbial moisture dependences and responses to drying–rewetting: the legacy of 18 years drought. Glob Change Biol. 2019;25:1005–15.
Article  Google Scholar 

26.
Evans SE, Wallenstein MD. Soil microbial community response to drying and rewetting stress: does historical precipitation regime matter? Biogeochemistry 2012;109:101–16.
Article  Google Scholar 

27.
Butterly CR, Bunemann EK, McNeill AM, Baldock JA, Marschner P. Carbon pulses but not phosphorus pulses are related to decreases in microbial biomass during repeated drying and rewetting of soils. Soil Biol Biochem. 2009;41:1406–16.
CAS  Article  Google Scholar 

28.
Schimel JP, Clein JS. Microbial response to freeze-thaw cycles in tundra and taiga soils. Soil Biol Biochem. 1996;28:1061–6.
Article  Google Scholar 

29.
Jurburg SD, Nunes I, Brejnrod A, Jacquiod S, Priemé A, Sørensen SJ, et al. Legacy effects on the recovery of soil bacterial communities from extreme temperature perturbation. Front Microbiol. 2017;8:1832.
PubMed  PubMed Central  Article  Google Scholar 

30.
Calderón K, Philippot L, Bizouard F, Breuil M-C, Bru D, Spor A. Compounded disturbance chronology modulates the resilience of soil microbial communities and n-cycle related functions. Front Microbiol. 2018;9:2721.
PubMed  PubMed Central  Article  Google Scholar 

31.
Wilson SL, Grogan P, Walker VK. Prospecting for ice association: characterization of freeze–thaw selected enrichment cultures from latitudinally distant soils. Can J Microbiol. 2012;58:402–12.
CAS  PubMed  Article  Google Scholar 

32.
Evans SE, Wallenstein MD, Burke IC. Is bacterial moisture niche a good predictor of shifts in community composition under long-term drought? Ecology. 2014;95:110–22.
PubMed  Article  Google Scholar 

33.
Yergeau E, Kowalchuk GA. Responses of antarctic soil microbial communities and associated functions to temperature and freeze-thaw cycle frequency. Environ Microbiol. 2008;10:2223–35.
PubMed  Article  Google Scholar 

34.
Lennon JT, Aanderud ZT, Lehmkuhl BK, Schoolmaster DR. Mapping the niche space of soil microorganisms using taxonomy and traits. Ecology. 2012;93:1867–79.
PubMed  Article  Google Scholar 

35.
Mackey BM, Derrick CM. Conductance measurements of the lag phase of injured salmonella-typhimurium. J Appl Bacteriol. 1984;57:299–308.
CAS  PubMed  Article  Google Scholar 

36.
Nocker A, Fernandez PS, Montijn R, Schuren F. Effect of air drying on bacterial viability: a multiparameter viability assessment. J Microbiol Methods. 2012;90:86–95.
CAS  PubMed  Article  Google Scholar 

37.
Potts M. Desiccation tolerance of prokaryotes. Microbiological Rev. 1994;58:755–805.
CAS  Article  Google Scholar 

38.
Öquist MG, Sparrman T, Klemedtsson L, Drotz SH, Grip H, Schleucher J, et al. Water availability controls microbial temperature responses in frozen soil co2 production. Glob Change Biol. 2009;15:2715–22.
Article  Google Scholar 

39.
Manzoni S, Schimel JP, Porporato A. Responses of soil microbial communities to water-stress: results from a meta-analysis. Ecology. 2012;93:930–8.
PubMed  Article  Google Scholar 

40.
Iovieno P, Bååth E. Effect of drying and rewetting on bacterial growth rates in soil. FEMS Microbiol Ecol. 2008;65:400–7.
CAS  PubMed  Article  Google Scholar 

41.
Drotz SH, Sparrman T, Nilsson MB, Schleucher J, Öquist MG. Both catabolic and anabolic heterotrophic microbial activity proceed in frozen soils. Proc Natl Acad Sci USA. 2010;107:21046–51.
CAS  PubMed  Article  Google Scholar 

42.
Lennon JT, Jones SE. Microbial seed banks: the ecological and evolutionary implications of dormancy. Nat Rev Microbiol. 2011;9:119–30.
CAS  PubMed  Article  Google Scholar 

43.
Manzoni S, Schaeffer SM, Katul G, Porporato A, Schimel JP. A theoretical analysis of microbial eco-physiological and diffusion limitations to carbon cycling in drying soils. Soil Biol Biochem. 2014;73:69–83.
CAS  Article  Google Scholar 

44.
Larsen KS, Jonasson S, Michelsen A. Repeated freeze-thaw cycles and their effects on biological processes in two arctic ecosystem types. Appl Soil Ecol. 2002;21:187–95.
Article  Google Scholar 

45.
Williams MA, Xia K. Characterization of the water soluble soil organic pool following the rewetting of dry soil in a drought-prone tallgrass prairie. Soil Biol Biochem. 2009;41:21–8.
CAS  Article  Google Scholar 

46.
Placella SA, Brodie EL, Firestone MK. Rainfall-induced carbon dioxide pulses result from sequential resuscitation of phylogenetically clustered microbial groups. Proc Natl Acad Sci USA. 2012;109:10931–6.
CAS  PubMed  Article  Google Scholar 

47.
Fukami T. Historical contingency in community assembly: Integrating niches, species pools, and priority effects. Annu Rev Ecol Evol Syst. 2015;46:1–23.
Article  Google Scholar 

48.
Szekely AJ, Langenheder S. Dispersal timing and drought history influence the response of bacterioplankton to drying-rewetting stress. ISME J. 2017;11:1764–76.
PubMed  PubMed Central  Article  Google Scholar 

49.
Engelhardt IC, Welty A, Blazewicz SJ, Bru D, Rouard N, Breuil M-C, et al. Depth matters: effects of precipitation regime on soil microbial activity upon rewetting of a plant-soil system. ISME J. 2018;12:1061–71.
CAS  PubMed  PubMed Central  Article  Google Scholar 

50.
Guhr A, Borken W, Spohn M, Matzner E. Redistribution of soil water by a saprotrophic fungus enhances carbon mineralization. Proc Natl Acad Sci USA. 2015;112:14647–51.
CAS  PubMed  Article  PubMed Central  Google Scholar 

51.
Yergeau E, Hogues H, Whyte LG, Greer CW. The functional potential of high arctic permafrost revealed by metagenomic sequencing, qpcr and microarray analyses. ISME J. 2010;4:1206–14.
CAS  PubMed  Article  PubMed Central  Google Scholar 

52.
Gittel A, Barta J, Kohoutova I, Mikutta R, Owens S, Gilbert J, et al. Distinct microbial communities associated with buried soils in the siberian tundra. ISME J. 2014;8:841–53.
CAS  PubMed  Article  PubMed Central  Google Scholar 

53.
Hultman J, Waldrop MP, Mackelprang R, David MM, McFarland J, Blazewicz SJ, et al. Multi-omics of permafrost, active layer and thermokarst bog soil microbiomes. Nature 2015;521:208–12.
CAS  PubMed  Article  PubMed Central  Google Scholar 

54.
Feng X, Nielsen LL, Simpson MJ. Responses of soil organic matter and microorganisms to freeze–thaw cycles. Soil Biol Biochem. 2007;39:2027–37.
CAS  Article  Google Scholar 

55.
Carini P, Marsden PJ, Leff J, Morgan EE, Strickland MS, Fierer N. Relic DNA is abundant in soil and obscures estimates of soil microbial diversity. Nat Microbiol. 2017;2:6.
Article  CAS  Google Scholar 

56.
Lennon JT, Muscarella ME, Placella SA, Lehmkuhl BK. How, when, and where relic DNA affects microbial diversity. mBio 2018;9:e00637–18.
PubMed  PubMed Central  Google Scholar 

57.
Blazewicz SJ, Barnard RL, Daly RA, Firestone MK. Evaluating rrna as an indicator of microbial activity in environmental communities: limitations and uses. ISME J. 2013;7:2061–8.
CAS  PubMed  PubMed Central  Article  Google Scholar 

58.
Blagodatskaya E, Kuzyakov Y. Active microorganisms in soil: critical review of estimation criteria and approaches. Soil Biol Biochem. 2013;67:192–211.
CAS  Article  Google Scholar 

59.
van de Voorde TFJ, van der Putten WH, Martijn Bezemer T. Intra- and interspecific plant–soil interactions, soil legacies and priority effects during old-field succession. J Ecol. 2011;99:945–53.
Article  Google Scholar 

60.
Nunes I, Jurburg S, Jacquiod S, Brejnrod A, Falcão Salles J, Priemé A, et al. Soil bacteria show different tolerance ranges to an unprecedented disturbance. Biol Fertil Soils. 2018;54:189–202.
Article  Google Scholar 

61.
Christiansen CT, Haugwitz MS, Priemé A, Nielsen CS, Elberling B, Michelsen A, et al. Enhanced summer warming reduces fungal decomposer diversity and litter mass loss more strongly in dry than in wet tundra. Glob Change Biol. 2017;23:406–20.
Article  Google Scholar 

62.
Borg Dahl M, Brejnrod AD, Russel J, Sørensen SJ, Schnittler M. Different degrees of niche differentiation for bacteria, fungi, and myxomycetes within an elevational transect in the german alps. Micro Ecol. 2019;78:764–80.
CAS  Article  Google Scholar 

63.
RCoreTeam. R: a language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing; 2018.
Google Scholar 

64.
McMurdie PJ, Holmes S. Phyloseq: an r package for reproducible interactive analysis and graphics of microbiome census data. PLOS ONE. 2013;8:e61217.
CAS  PubMed  PubMed Central  Article  Google Scholar 

65.
Oksanen J, Blanchet FG, Friendly M, Kindt R, Legendre P, McGlinn D, et al. Vegan: community ecology package. R package version 2.5-5. https://CRAN.R-project.org/package=vegan. 2019.

66.
Kuznetsova A, Brockhoff PB, Christensen RHB. lmerTest Package: Tests in Linear Mixed Effects Models. J Stat Softw. 2017;82:1–26.
Article  Google Scholar 

67.
Bates D, Maechler M, Bolker B, Walker S. Fitting Linear Mixed-Effects Models Using lme4. J Stat Softw. 2015;67:1–48.
Article  Google Scholar 

68.
Warnes GR, Bolker B, Bonebakker L, Gentleman R, Liaw WHA, Lumley T, et al. Package ‘gplots’: various r programming tools for plotting data. R package version 3.0.1. https://CRAN.R-project.org/package=gplots. 2019.

69.
Wickham H. Ggplot2: elegant graphics for data analysis. New York: Springer-Verlag; 2016.
Google Scholar 

70.
Sharma S, Szele Z, Schilling R, Munch JC, Schloter M. Influence of freeze-thaw stress on the structure and function of microbial communities and denitrifying populations in soil. Appl Environ Microbiol. 2006;72:2148–54.
CAS  PubMed  PubMed Central  Article  Google Scholar 

71.
Tilston EL, Sparrman T, Öquist MG. Unfrozen water content moderates temperature dependence of sub-zero microbial respiration. Soil Biol Biochem. 2010;42:1396–407.
CAS  Article  Google Scholar 

72.
Ho A, Di Lonardo DP, Bodelier PLE. Revisiting life strategy concepts in environmental microbial ecology. FEMS Microbiol Ecol. 2017;93:14.
Google Scholar 

73.
Yuste JC, Penuelas J, Estiarte M, Garcia-Mas J, Mattana S, Ogaya R, et al. Drought-resistant fungi control soil organic matter decomposition and its response to temperature. Glob Change Biol. 2011;17:1475–86.
Article  Google Scholar 

74.
Kaisermann A, Maron PA, Beaumelle L, Lata JC. Fungal communities are more sensitive indicators to non-extreme soil moisture variations than bacterial communities. Appl Soil Ecol. 2015;86:158–64.
Article  Google Scholar 

75.
Clemmensen KE, Bahr A, Ovaskainen O, Dahlberg A, Ekblad A, Wallander H, et al. Roots and associated fungi drive long-term carbon sequestration in boreal forest. Science. 2013;339:1615–8.
CAS  PubMed  Article  Google Scholar 

76.
Bates ST, Clemente JC, Flores GE, Walters WA, Parfrey LW, Knight R, et al. Global biogeography of highly diverse protistan communities in soil. ISME J. 2013;7:652–9.
CAS  PubMed  Article  Google Scholar 

77.
Geisen S, Mitchell EAD, Adl S, Bonkowski M, Dunthorn M, Ekelund F, et al. Soil protists: a fertile frontier in soil biology research. FEMS Microbiol Rev. 2018;42:293–323.
CAS  PubMed  Article  Google Scholar 

78.
Geisen S, Bandow C, Römbke J, Bonkowski M. Soil water availability strongly alters the community composition of soil protists. Pedobiologia 2014;57:205–13.
Article  Google Scholar 

79.
Clarholm M. Protozoan grazing of bacteria in soil—impact and importance. Micro Ecol. 1981;7:343–50.
CAS  Article  Google Scholar 

80.
Harder CB, Rønn R, Brejnrod A, Bass D, Al-Soud WA, Ekelund F. Local diversity of heathland cercozoa explored by in-depth sequencing. ISME J. 2016;10:2488.
CAS  PubMed  PubMed Central  Article  Google Scholar 

81.
Ekelund F, Frederiksen HB, Rønn R. Population dynamics of active and total ciliate populations in arable soil amended with wheat. Appl Environ Microbiol. 2002;68:1096–101.
CAS  PubMed  PubMed Central  Article  Google Scholar 

82.
Müller H, Achilles-Day UEM, Day JG. Tolerance of the resting cysts of colpoda inflata (ciliophora, colpodea) and meseres corlissi (ciliophora, spirotrichea) to desiccation and freezing. Eur J Protistol. 2010;46:133–42.
PubMed  Article  Google Scholar 

83.
Stieglmeier M, Alves RJE, Schleper C. The phylum thaumarchaeota. In: Rosenberg E, DeLong EF, Lory S, Stackebrandt E, Thompson F, editors. The prokaryotes: other major lineages of bacteria and the archaea. Berlin, Heidelberg: Springer Berlin Heidelberg; 2014. p. 347–62.

84.
Shi Y, Adams JM, Ni Y, Yang T, Jing X, Chen L, et al. The biogeography of soil archaeal communities on the eastern tibetan plateau. Sci Rep. 2016;6:38893.
CAS  PubMed  PubMed Central  Article  Google Scholar 

85.
Pedrós-Alió C. The rare bacterial biosphere. Ann Rev Mar Sci. 2012;4:449–66.
PubMed  Article  Google Scholar 

86.
Lynch MDJ, Neufeld JD. Ecology and exploration of the rare biosphere. Nat Rev Microbiol. 2015;13:217–29.
CAS  PubMed  Article  Google Scholar 

87.
Chowdhury N, Nakatani AS, Setia R, Marschner P. Microbial activity and community composition in saline and non-saline soils exposed to multiple drying and rewetting events. Plant Soil. 2011;348:103–13.
CAS  Article  Google Scholar 

88.
Fierer N, Schimel JP. Effects of drying-rewetting frequency on soil carbon and nitrogen transformations. Soil Biol Biochem. 2002;34:777–87.
CAS  Article  Google Scholar 

89.
Schöler A, Jacquiod S, Vestergaard G, Schulz S, Schloter M. Analysis of soil microbial communities based on amplicon sequencing of marker genes. Biol Fertil Soils. 2017;53:485–9.
Article  CAS  Google Scholar 

90.
Lau JA, Lennon JT. Rapid responses of soil microorganisms improve plant fitness in novel environments. Proc Natl Acad Sci USA. 2012;109:14058–62.
CAS  PubMed  Article  PubMed Central  Google Scholar 

91.
Meisner A, De Deyn GB, de Boer W, van der Putten WH. Soil biotic legacy effects of extreme weather events influence plant invasiveness. Proc Natl Acad Sci USA. 2013;110:9835–8.
CAS  PubMed  Article  PubMed Central  Google Scholar 

92.
Solyanikova IP, Suzina NE, Egozarian NS, Polivtseva VN, Prisyazhnaya NV, El-Registan GI, et al. The response of soil arthrobacter agilis lush13 to changing conditions: Transition between vegetative and dormant state. J Environ Sci Health B 2017;52:745–51.
CAS  PubMed  Article  PubMed Central  Google Scholar 

93.
Levinson HS, Hyatt MT. Correlation of respiratory activity with phases of spore germination and growth in bacillus-megaterium as influenced by manganese and l-alanine. J Bacteriol. 1956;72:176–83.
CAS  PubMed  PubMed Central  Article  Google Scholar 

94.
Williamson KE, Fuhrmann JJ, Wommack KE, Radosevich M. Viruses in soil ecosystems: an unknown quantity within an unexplored territory. Annu Rev Virol. 2017;4:201–19.
CAS  PubMed  Article  PubMed Central  Google Scholar 

95.
van Kruistum H, Bodelier PLE, Ho A, Meima-Franke M, Veraart AJ. Resistance and recovery of methane-oxidizing communities depends on stress regime and history; a microcosm study. Front Microbiol. 2018;9:1714.
PubMed  PubMed Central  Article  Google Scholar 

96.
Bérard A, Ben Sassi M, Renault P, Gros R. Severe drought-induced community tolerance to heat wave. An experimental study on soil microbial processes. J Soils Sed. 2012;12:513–8.
Article  Google Scholar  More

Data selection and handling: death data
For mortality due to COVID-19, we used time series provided by the New York Times12. We selected the New York Times dataset because it is rigorously curated. We analyzed separately only counties that had records of 100 or more deaths by 23 May, 2020. The threshold of 100 was a balance between including more counties and obtaining reliable estimates of r(t). Preliminary simulations showed that time series with low numbers of deaths would bias r(t) estimates (Supplementary Fig. 2). However, we did not want to use the maximum daily number of deaths as a selection criterion, because this could lead to selection of counties based on data from a single day. It would also involve some circularity, because the information obtained, r(t), would be directly related to the criterion used to select datasets. Therefore, we used the threshold of 100 cumulative deaths. The District of Columbia was treated as a county. Also, because the New York Times dataset aggregated the five boroughs of New York City, we treated them as a single county. For counties with fewer than 100 deaths, we aggregated mortality to the state level to create a single time series. For thirteen states (AK, DE, HI, ID, ME, MT, ND, NH, SD, UT, VM, WV, and WY), the aggregated time series did not contain 100 or more deaths and were therefore not analyzed.
Data selection and handling: explanatory county-level variables
County-level variables were collected from several public data sources36,37,38,39,40,41,42. We selected socio-economic variables a priori in part to represent a broad set of population characteristics.
Time series analysis: time series model
We used a time-varying autoregressive model15,56 designed explicitly to estimate the rate of increase of a variable using nonlinear, state-dependent error terms16. We assume in our analyses that the susceptible proportion of the population represented by a time series is close to one, and therefore there is no decrease in the infection rate caused by a pool of individuals who were infected, recovered, and were then immune to further infection.
The model is

$$xleft( t right) = rleft( {t-1} right) + xleft( {t-1} right)$$
(1a)

$$rleft( t right) = rleft( {t-1} right) + omega_rleft( t right)$$
(1b)

$$Dleft( t right) = {mathrm{exp}}(xleft( t right) + phi left( t right))$$
(1c)

Here, x(t) is the unobserved, log-transformed value of daily deaths at time t, and D(t) is the observed count that depends on the observation uncertainty described by the random variable ϕ(t). Because a few of the datasets that we analyzed had zeros, we replaced zeros with 0.5 before log-transformation. The model assumes that the death count increases exponentially at rate r(t), where the latent state variable r(t) changes through time as a random walk with ωr(t) ~ N(0, σ2r). We assume that the count data follow a quasi-Poisson distribution. Thus, the expectation of counts at time t is exp(x(t)), and the variance is proportional to this expectation.
We fit the model using the extended Kalman filter to compute the maximum likelihood57,58. In addition to the parameters σ2r and σ2ϕ, we estimated the initial value of r(t) at the start of the time series, r0, and the initial value of x(t), x0. The estimation also requires terms for the variances in x0 and r0, which we assumed were zero and σ2r, respectively. In the validation using simulated data (Supplementary Methods: Simulation model), we found that the estimation process tended to absorb σ2r to zero too often. To eliminate this absorption to zero, we imposed a minimum of 0.02 on σ2r.
Time series analysis: parametric bootstrapping
To generate approximate confidence intervals for the time-varying estimates of r(t) (Eq. 1b), we used a parametric bootstrap designed to simulate datasets with the same characteristics as the real data that are then refit using the autoregressive model. We used bootstrapping to obtain confidence intervals, because an initial simulation study showed that standard methods, such as obtaining the variance of r(t) from the Kalman filter, were too conservative (the confidence intervals too narrow) when the number of counts was small. Furthermore, parametric bootstrapping can reveal bias and other features of a model, such as the lags we found during model fitting (Supplementary Fig. 1a, b).
Changes in r(t) consist of unbiased day-to-day variation and the biased deviations that lead to longer-term changes in r(t). The bootstrap treats the day-to-day variation as a random variable while preserving the biased deviations that generate longer-term changes in r(t). Specifically, the bootstrap was performed by calculating the differences between successive estimates of r(t), Δr(t) = r(t) – r(t-1), and then standardizing to remove the bias, Δrs(t) = Δr(t) – E[Δr(t)], where E[] denotes the expected value. The sequence Δrs(t) was fit using an autoregressive time-series model with time lag 1, AR(1), to preserve any shorter-term autocorrelation in the data. For the bootstrap, a new time series was simulated from this AR(1) model, Δρ(t), and then standardized, Δρs(t) = Δρ(t) – E[Δρ(t)]. The simulated time series for the spread rate was constructed as ρ(t) = r(t) + Δρs(t)/21/2, where dividing by 21/2 accounts for the fact that Δρs(t) was calculated from the difference between successive values of r(t). A new time series of count data, ξ(t), was then generated using equation 1 with the parameters from fitting the data. Finally, the statistical model was fit to the reconstructed ξ(t). In this refitting, we fixed the variance in r(t), σ2r, to the same value as estimated from the data. Therefore, the bootstrap confidence intervals are conditional of the estimate of σ2r.
Time series analysis: calculating R0
We derived estimates of R(t) directly from r(t) using the Dublin-Lotka equation21 from demography. This equation is derived from a convolution of the distribution of births under the assumption of exponential population growth. In our case, the “birth” of COVID-19 is the secondary infection of susceptible hosts leading to death, and the assumption of exponential population growth is equivalent to assuming that the initial rate of spread of the disease is exponential, as is the case in equation 1. Thus,

$$Rleft( t right) = 1/mathop {sum}nolimits_{_tau} {{mathrm{e}}^{ – r(t)}} tau p(tau)$$
(2)

where p(τ) is the distribution of the proportion of secondary infections caused by a primary infection that occurred τ days previously. We used the distribution of p(τ) from Li et al.59 that had an average serial interval of T0 = 7.5 days; smaller or larger values of T0, and greater or lesser variance in p(τ), will decrease or increase R(t) but will not change the pattern in R(t) through time. Note that the uncertainty in the distribution of serial times for COVID-19 is a major reason why we focus on estimating r0, rather than R0: the estimates of r0 are not contingent on time distributions that are poorly known. Computing R(t) from r(t) also does not depend on the mean or variance in time between secondary infection and death. We report values of R(t) at dates that are offset by 18 days, the average length of time between initial infection and death given by Zhou et al.60.
Time series analysis: Initial date of the time series
Many time series consisted of initial periods containing zeros that were uninformative. As the initial date for the time series, we chose the day on which the estimated daily death count exceeded 1. To estimate the daily death count, we fit a Generalized Additive Mixed Model (GAMM) to the death data while accounting for autocorrelation and greater measurement error at low counts using the R package mgcv61. We used this procedure, rather than using a threshold of the raw death count, because the raw death count will include variability due to sampling small numbers of deaths. Applying the GAMM to “smooth” over the variation in count data gives a well-justified method for standardizing the initial dates for each time series.
Time series analysis: validation
We performed extensive simulations to validate the time-series analysis approach (Supplementary Methods: Simulation model).
Regression analysis for r 0
We applied a Generalized Least Squares (GLS) regression model to explain the variation in estimates of r0 from the 160 county and county-aggregate time series:

$$r_0 = b_0 + b_1start.date + b_2logleft( {pop.size} right) + b_3pop.den^{0.25} + varepsilon$$
(3)

where start.date is the Julian date of the start of the time series, log(pop.size) and pop.den0.25 are the log-transformed population size and 0.25 power-transformed population density of the county or county-aggregate, respectively, and ε is a multivariate Gaussian random variable with covariance matrix σ2Σ. We used the transforms log(pop.size) and pop.den0.25 to account for nonlinear relationships with r0; these transforms give the highest maximum likelihood of the overall regression. The covariance matrix contains a spatial correlation matrix of the form C = uI + (1–u)S(g) where u is the nugget and S(g) contains elements exp(−dij/g), where dij is the distance between spatial locations and g is the range62. To incorporate differences in the precision of the estimates of r0 among time series, we weighted by the vector of their standard errors, s, so that Σ = diag(s) * C * diag(s), where * denotes matrix multiplication. With this weighting, the overall scaling term for the variance, σ2, will equal 1 if the residual variance of the regression model matches the square of the standard errors of the estimates of r0 from the time series. We fit the regression model with the function gls() in the R package nlme63.
To make predictions for new values of r0, we used the relationship

$$hat e_{mathrm{i}} = bar e + {mathbf{v}}_{mathbf{i}} ast ,{mathbf{V}}^{ – 1}(epsilon_i – bar e)$$
(4)

where ει is the GLS residual for data i, (hat e)i is the predicted residual, (bar e) is the mean of the GLS residuals, V is the covariance matrix for data other than i, and vi is a row vector containing the covariances between data i and the other data in the dataset64. This equation was used for three purposes. First, we used it to compute R2pred for the regression model by removing each data point, recomputing (hat e)i, and using these values to compute the predicted residual variance23. Second, we used it to obtain predicted values of r0, and subsequently R0, for the 160 counties and county-aggregates for which r0 was also estimated from time series. Third, we used equation (4) to obtain predicted values of r0, and hence predicted R0, for all other counties. We also calculated the variance of the estimates from64

$$hat v_{mathrm{i}} = sigma^2-{mathbf{v}}_{mathbf{i}} ast ,{mathbf{V}}^{ – 1} ast v_i^{mathbf{t}}$$
(5)

Predicted values of R0 were mapped using the R package usmap65.
Regression analysis for SARS-CoV-2 effects on r0
The GISAID metadata27 for SARS-CoV-2 contains the clade and state-level location for strains in the USA; strains G, GH, and GR contain the G614 mutation. For each state, we limited the SARS-CoV-2 genomes to those collected no more than 30 days following the onset of outbreak that we used as the starting point for the time series from which we estimated r0; from these genomes (totaling 5290 from all states), we calculated the proportion that had the G614 mutation. We limited the analyses to the 28 states that had five or more genome samples. For each state, we selected the estimates of r0 from the county or county-aggregate representing the greatest number of deaths. We fit these estimates of r0 with the weighted Least Squares (LS) model as in equation (3) with additional variables for strain. Figure 3 was constructed using the R packages usmap65 and scatterpie66.
Statistics and reproducibility
The statistics for this study are summarized in the preceding sections of the “Methods”. No experiments were conducted, so experimental reproducibility is not an issue. Nonetheless, we repeated analyses using alternative datasets giving county-level characteristics, and also an alternative dataset on SARS-CoV-2 strains (Supplementary Methods: Analysis of Nextstrain metadata of SARS-CoV-2 strains), and all of the conclusions were the same.
Reporting summary
Further information on research design is available in the Nature Research Reporting Summary linked to this article. More

As a wildlife ecologist, I study many species, but the platypus (Ornithorhynchus anatinus) is my speciality. I’ve worked with these duck-billed monotremes, which are indigenous to eastern Australia, for 13 years. So much remains unknown about them — as solitary animals, they are very hard to capture. The males have a venomous spur on each of their hind feet, so I always hold them by the tail. Luckily, I’ve never been spurred.
At the environmental consultancy firm Cesar Australia, we’ve focused on whether platypuses’ numbers are rising or falling. Their aquatic habitat has shrunk by 25% in the past 30 years — partly owing to reduced rainfall from climate change and partly to residential and industrial development, leaving less water, which is of poorer quality.
To monitor the population’s size and health, we do overnight trapping surveys — typically during the spring breeding season and in the autumn, when the babies emerge. For each survey, we set up nets along 5–10 kilometres of waterway and then drive round checking them through the night. It’s a long, sleepless process that usually yields just one or two catches. In the 15 minutes or so for which each animal is out of the water, we weigh and measure it. If it’s carrying an identification chip, we know we’ve caught it before. If not, we implant one in the animal. In this picture, I am releasing a platypus back into a stream outside Melbourne.
In the past five years, we’ve also been sampling environmental DNA for signatures of platypuses’ presence. It’s been a game-changer in terms of mapping their distribution. We can take a water sample, extract all the DNA and confidently identify whether platypus DNA is present. We’ve surveyed more than 2,000 sites spanning 500,000 square kilometres, enabling us to measure at scales never before possible.
We’ve also developed a citizen-science platform through which people can submit platypus sightings. This helps us still further. More