Ecology

We retrieved 684 Lokiarchaeota and 31 Heimdallarchaeota near-full-length 16S rRNA sequences from sequence libraries generated from sediment sampled at 27 m water depth in 5 cm intervals between 0 and 40 cm below seafloor (cm.b.s.f) in Aarhus Bay (Supplementary Information). The maximum relative read abundance of Lokiarchaeota was 1.6% at 15–20 cm.b.s.f. and 0.1% for Heimdallarchaeota at 10–15 cm.b.s.f. (Fig. 1). The sequences were grouped into 58 Loki- and 3 Heimdallarchaeota operational taxonomic units (OTUs) using a 98% sequence identity threshold and formed three distinct Lokiarchaeota clades and one monophyletic Heimdallarchaeota cluster (Fig. 1). The primer-free sequencing of RNA extracts enabled us to broadly sample the Asgard archaeal diversity in Aarhus Bay sediments and provided a solid database to design oligonucleotide probes for their visualization.Fig. 1: Phylogenetic analysis and depth distribution of Loki- and Heimdallarchaeota 16S rRNA sequences from Aarhus Bay sediments.A Maximum likelihood phylogeny of Loki- and Heimdallarchaeota operational taxonomic units (OTUs) and related sequences selected from the SILVA database (v. 132). Specificities of FISH probes and the number of sequences constituting each OTU are also depicted. TACK archaea were selected as outgroup. Bar: 0.1 substitutions per nucleotide position. B Heatmap and relative abundances of Loki- and Heimdallarchaeota sequences at different sediment depths.Full size imageBased on the newly retrieved full-length sequences, we designed four novel oligonucleotide probes specifically targeting Loki- and Heimdallarchaeota 16S rRNA with high coverage (Fig. 1, Supplementary Table 1). Probe LOK1183 targets almost all sequences in Lokiarchaeota Clade A, which contains 92% of the retrieved Lokiarchaeota sequences from Aarhus Bay sediments, while probe LOK1378 targets 85% of the sequences in all three Lokiarchaeota clades. Probe HEIM329 and HEIM529 each target >97% of the retrieved Heimdallarchaeota sequences. All designed probes cover >89% sequences in their target groups in the SILVA database (v. 132). The two Lokiarchaeota probes match 5 and 10 different non-target sequences in the SILVA database (v. 132), respectively, while the Heimdallarchaeota probes have no match outside their target group. The broad coverage and high specificity suggest that our probes can also be used to detect Loki- and Heimdallarchaeota in other habitats. Furthermore, designing two probes for each phylum enabled us to identify Lokiarchaeota clade A and Heimdallarchaeota cells via double hybridizations with two distinct dyes and thus confidently distinguish true- and false-positive signals (Supplementary Fig. 1). The general archaeal probe ARC915 also targets Lokiarchaeota and thereby provided yet another control for specific hybridization of the two Lokiarchaeota-specific probes, while the non-sense probe NON338 served as the negative control. We also designed competitor probes to minimize the theoretical false-positive hybridizations with the most frequent one and two mismatches [11] in the SILVA database (v. 132) and helper probes to facilitate probe binding [12]. This comprehensive experimental design with appropriate controls enabled reliable detection of low-abundant Loki- and Heimdallarchaeota cells in Aarhus Bay sediments.We used both confocal laser scanning microscopy (CLSM) and three-dimensional super-resolution structured illumination (SR-SIM) microscopy for detailed imaging of dual-labeled Loki- and Heimdallarchaeota signals. Loki- and Heimdallarchaeota cells featured coccoid shapes and often formed clusters (Fig. 2) (Supplementary Fig. 2). Based on SR-SIM imaging, Lokiarchaeota cells (n = 18) were 1.27 ± 0.24 µm in diameter and 1.43 ± 0.25 µm in length, while the width and the length of Heimdallarchaeota cells (n = 11) were 1.30 ± 0.20 µm and 1.37 ± 0.21 µm, respectively (Supplementary Table 2). In addition, we observed a few large ( >3 µm) ovoid and filamentous cells, resembling some of the Lokiarchaeota morphotypes reported from lake sediment [9]; however, we never detected these cell types in double hybridizations with two probes (Supplementary Fig. 1P–R), and therefore consider them false-positives.Fig. 2: Visualization of Loki- and Heimdallarchaeota cells in Aarhus Bay sediments by CARD-FISH.Probe names and the dyes are indicated for each panel. Representative cell morphotypes were imaged in a super-resolution structured illumination microscope (SR-SIM; panels (A), (B), (D), (E)) and confocal laser scanning microscope (CLSM; panels (C) and (F)). For SR-SIM images, single slices from the center of the focal plane are shown. For CLSM images, three-dimensional (3D) surface reconstructions are depicted. All z-stack images taken in CLSM are included in Supplementary Fig. 2. 360° rotation of 3D reconstructed images are also provided in Supplementary Video. Negative and positive controls are shown in Supplementary Fig. 1 together with large ovoid and filamentous false-positive signals. Images are representative of dual labeled Lokiarchaeota (n = 72) and Heimdallarchaeota (n = 70) cells in five individual experiments using two different sediment cores taken from the same sampling site. The scale bar is 1 µm.Full size imageThe DNA stain (4′,6-diamidino-2-phenylindole; DAPI) in the FISH-identified Loki- and Heimdallarchaeota cells was consistently confined to a single spherical central or lateral position (Fig. 2), corroborating the signal pattern suggested for some of the Asgard archaeal cells in lake sediments [9]. Using SR-SIM, we could image a clear gap, which separated the DNA from the ribosome-originated FISH signals with an average width of 0.18 ± 0.07 µm in Heimdallarchaeota and 0.16 ± 0.13 µm in Lokiarchaeota cells (Supplementary Table 2). The spatial separation of DNA and ribosomes in Loki- and Heimdallarchaeota cells represents an unusual observation since DAPI and FISH signals generally overlap partially or completely in prokaryotic cells [13]. Accordingly, SR-SIM imaging of benthic bacteria in Aarhus Bay sediments demonstrated the prevalence of this overlapping signal pattern (Supplementary Fig. 3). Also, the separated DNA signal observed in Loki- and Heimdallarchaeota cells appeared different from the condensed DNA formation previously described, for example, in Escherichia coli cells [14] and the Thaumarcheota Cenarcheum symbiosum [15] and Nitrosopumilus maritimus [16]. To corroborate this, we performed SR-SIM imaging of CARD-FISH-labeled E.coli and N. maritimus cells. Although their DNA was condensed in particular cellular locations, their FISH and DAPI signals always overlapped, indicating that their DNA and ribosomes are partially co-localized and not fully separated (Supplementary Fig. 4).To assess whether the gap between DAPI and FISH signals was indicative of an internal membrane, we tried various dyes to stain membranes of the CARD-FISH-labeled Asgard archaeal cells (Supplementary Information). However, none of these stainings was successful, not even for the outer cell membrane, most likely because cell membranes were disintegrated during the CARD-FISH protocol. We then used wheat germ agglutinin (WGA), a lectin primarily binding to N-acetyl-D-glucosamine but also other glycoconjugates and oligosaccharides [17] to at least be able to visualize the surfaces of Loki- and Heimdallarchaeota cells. WGA consistently decorated a cell surface that enclosed the proximal FISH and DAPI signals, suggesting that both signals originated from the same single cell (Supplementary Fig. 5). The WGA staining also demonstrated extracellular structures connected to some Heimdallarchaeota cells (Supplementary Fig. 5). These structures appear different than the membrane protrusions in the first cultured Lokiarchaeon “Ca. P. syntrophicum”, which has a considerably smaller cell size (550 nm in diameter) and does not possess the separated DNA and ribosome signals [5]. Our observations therefore indicate diverse cellular organizations and morphotypes within Asgard archaea superphylum.Our combined results suggest that genomic material is condensed and spatially distinct from the riboplasm within the detected Loki- and Heimdallarchaeota cells. Considering the anticipated role of Asgard archaea in eukaryogenesis, in particular the presence of ESPs potentially involved in dynamic cytoskeleton formation [18] and membrane remodeling [4, 19], the separation of DNA- and ribosome-derived signals might be indicative of cellular compartmentalization. Alternatively, the observed pattern could be the result of a membrane invagination to form a nucleoid region, similar to membrane organizations for example in Planctomycetes cells [20] or Atribacter laminatus [21].Our study demonstrates the first visualization of diverse Loki- and Heimdallarchaeota cells in the marine environment and provides a protocol for reliable in situ imaging of rare microorganisms in environmental samples. Future research should address the ultrastructure of Asgard archaeal cells using electron microscopy. This would help to elucidate the cell biology of Asgard archaea and provide insights into the emergence of subcellular complexity of the eukaryotic cell. More

The following experimental workflow was implemented to address the main questions raised in our study (Fig. 1).Fig. 1: The scheme of experimental pipeline used in this study to examine the impact of lytic phage infection on the P. aeruginosa population and the development of phage-resistance.Experiments were conducted as follows: culture preparation (1); biofilm formation (2); phage infection with single or cocktail preparations (3); incubation (4); biofilm and planktonic populations sampling (5); culture plating on TSA agar and isolation of discrete colonies (6); phage typing determination (7); to select isolates with unique patterns (8) for further phenotypic (9) and genome sequencing analyses (10).Full size imageThe P. aeruginosa PAO1 reference strain and four other clinical representatives were infected with distinct lytic phages in a single or different cocktail combination. Randomly picked colonies from the surviving cultures were then tested in terms of susceptibility to inoculated phages as well as to the others from the Pseudomonas phages panel (Table 1). We were interested in exploring the broadest clonal variability developed in phage infected Pseudomonas population. Therefore, the first phase of the study was focused on examining the phenotypic heterogeneity of PAO1 reference mutants (phage typing) within planktonic and biofilm populations. Since the consequences of introducing lytic phages into the bacterial population are difficult to predict, a representative pool of bacterial clones that have survived infection was sampled. A total of 780 P. aeruginosa PAO1 clones were typed with phages (planktonic (320), biofilm populations (400) and 60 control clones). No resistance to phages was observed among the control clones taken from untreated biofilm or plankton. Therefore, three biofilm and three planktonic representatives and the wild-type PAO1 were selected for further genetic and fitness analyses (Table S1). Finally, a pool of 95 isolates has been filtered, representing seventeen different phage susceptibility patterns (Tables S1, 2). This selection was based on the maximum variety of phage-type profiles, without accounting for the origin of the isolate (biofilm/plankton), as the infected planktonic bacteria turned out to be less diverse and all phage types were also present in the biofilm population.Since we did not aim to analyse the differences of planktonic versus sessile cells response to phage infection but rather look for maximum population heterogeneity, we decided to focus the investigation on the biofilm population for the other clinical strains during the second stage of this research. Accordingly, 880 (30 clones from every condition plus 10 control clones for each strain) isolated colonies from A5803, AA43, CHA, and PA biofilm populations were first subjected to phage typing. No phage resistance was observed among clones taken from phage-untreated samples compared to the wild-type strain. Ultimately, 35 phage-treated colonies, three controls, and the wild-type from each strain were selected for further investigation, resulting in a pool of 156 clones in total (39 × 4 strains) representing over twenty different phage susceptibility patterns (Table S1 and S3).Do phages always select for cross-resistance to other phages recognising the same bacterial receptor?The application of monovalent phage against reference PAO1 population generally led to the selection of cross-resistance against phages that recognise the same receptor as the applied one (Table S2). This was observed for 12/17 and 23/24 PAO1 clones isolated after LPS- and T4P-dependent phages treatment, respectively. Similar relation (15/20) was only observed for other clinical cultures infected with phiKZ phage (T4P-dependent) (Table S3). The resistance to both groups of phages was less frequent in monovalent infections (14.5% in PAO1 and 32.5% for other clinical strains) compared to polyvalent infections (61.1%; 33/54) and 51.6% (31/60) for PAO1 and clinical strains, respectively. The use of a cocktail of two phages recognising LPS selected for PAO1 clones resistant only to LPS-dependent phages. In contrast, LPS-dependent phages application was mostly accompanied by the emergence of resistance to phages recognising alternative receptors in clinical strains (28/60 cases).The introduction of a particular phage into the population did not guarantee the isolation of clones resistant to this phage. This event was recorded in the case of single phages, as well as for polyvalent combinations (23 PAO1 mutants). However, the cross-resistance to other phages recognising the same or both receptors did also occur. Interestingly, LUZ7 and KTN6 phages could still infect surviving clinical populations with a frequency of 23/60 and 44/80, respectively. Indicating that the resistance to LPS-dependent phages in clinical strains was more difficult to develop compared to those impaired by giant viruses, with 11/60 and 1/20 still sensitive to phiKZ and PA5oct phages, respectively. Almost all PAO1 (80/95) and clinical (127/140) clones treated with phages developed resistance to phage PA5oct, whereas the resistance to the entire phage panel emerged regardless of the single or cocktails application.To conclude, the selection of cross-resistance to other phages recognising the same bacterial receptor was mostly valid in the PAO1 model, whereas the other clinical strains primarily developed the cross-resistance to T4P-dependent phages.Do phages from different taxonomy groups recognising the same receptor cause the emergence of the same type of resistant mutants? Are the defence response and genome changes correlated with the receptor specificity of infecting phage?To assess the genetic basis of the resistance selected by phages, we performed single nucleotide polymorphisms (SNPs) and mapping analyses of 102 reference PAO1 clones and 156 clones derived from clinical strains (Figs. 2, 3, Table S2–S4). The wild-type P. aeruginosa strains were also re-sequenced with Illumina and PacBio technologies to ascertain their complete genomic background. Missense, nonsense, and frameshift mutation variants were taken into account in the analyses. Mutations that also occurred in control isolates were excluded from further consideration. The remaining mutations were divided into six groups: LPS-related genes, mucoidity-associated genes (EPS production, biofilm formation), T4P-related genes, global regulatory genes, and others (hypothetical or undefined function genes). The comparative analysis showed the presence of point mutations in 64 out of 95 examined PAO1 mutants. The frequency of mutations in PAO1 clones isolated after treatment with single or multiple phages was similar (73% vs 61%, respectively). In most of those isolates, only one gene mutation event was recorded (43%). However, in 23 cases SNPs occurred in two or three genes belonging to different metabolic gene groups. Five PAO1 isolates showed the presence of mutations in two genes from one gene group. The 33 cases of SNPs related to LPS synthesis were found in 29 mutants selected with single LPS-dependent phage preparations or in polyvalent combinations. Among these, the most frequent mutation (21/33 cases) was observed within the wzy gene, encoding the B-band O-antigen polymerase [30]. These frequent mutations in the LPS-biosynthesis cluster confirmed the phage resistance results emerging after LUZ7, KT28, and KTN6 phages propagation. In some cases, the LPS gene modification was accompanied by changes in EPS-related genes, leading to a mucoid phenotype. The T4P-dependent phage treatment also led to the selection of specific mutations in genes responsible for T4P expression, but also alterations in flagella-related genes (flgH, fliN, fliP, flhA). The mutations in global regulatory genes (most frequent yqjG and vfr) and “others” gene groups did not show any correlation to the type of phages used.Fig. 2: Graphical presentation of genetic changes occurring in the population of P. aeruginosa as a result of the infection by selected phages.The colour dots refer to particular gene groups where the point mutations (accumulated results) were recorded within the genomes of examined mutant clones. The lower line contains information on the maximum and minimum size of large deletions (grey bands) and the presence of intact prophages (light blue bands). * means mutation in promoter region of the gene.Full size imageFig. 3: The frequency of genetic changes per clone detected in P. aeruginosa strains.Panel (A) represents the PAO1 clones, and panel (B) represents the clinical strains populations. Populations were selected by specific phages targeting LPS (red dots) or T4P (blue dots) as a single treatment or in cocktails. The colour bars refer to particular gene groups where the point mutations were recorded within the genomes of examined mutant clones. N means the number of analysed clones for each strain.Full size imageApart from point mutations, 23% of phage-resistant PAO1 isolates contained large genomic deletions (23,983 bp–544,729 bp) appeared regardless of the phage-type and cocktail composition used as selective pressure agents. All deletions were located in the same region and despite different starting/ending points, they hold a core element of 19,038 bp. This core element carries the galU gene (responsible for LPS biosynthesis), as well as the hmgA gene, which causes the accumulation of brown pigment in bacterial cells when absent. Besides, the cumulated deletion range contained a total of 706,374 bp, including many key genes involved in the bacterial metabolism.Mutations detected in other clinical phage-resistant clones were classified according to the same criteria as in PAO1 (Figs. 2, 3, Table S3, S4). The genome-driven response to phage infection of A5803 was primarily located in global (71%, cpdA) and other genes (34%, PA2911); of AA43 in other genes (31%, PA2911); of CHA in T4P (34%) and global genes (34%, morA); and of PAK in T4P (25%) and other genes (23%, PA2911). Most of the mutations selected by LPS-dependent phage exposition were found in the global regulatory genes (9–11–25–54%) or “other” genes (17–23–31%), rather than in the LPS biosynthesis locus (0–3–6–17%) depending on the impacted strain (Table S3). That confirmed the phage-typing results where LUZ7 and KTN6 phages remained lytic towards surviving clones. In contrast, the application of phiKZ selected for the cross-resistance to T4P-dependent phages as well as for the genetic modifications in pili-related genes. Mutations in global regulatory and “others” genes show no correlation to the receptor specificity of phages used. Interestingly, a portion of phenotypically phage-resistant clones in each clinical P. aeruginosa population (5-9/35 clones) did not reveal any distinguishable genetic modifications. Consistent with PAO1, large genomic deletions were observed in A5803, AA43, and PAK strains ranging between 92,207 bp and 383,693 bp in size, encompassing the galU region. The MEME analysis of the regions flanking the deletions did not reveal specific motifs that would indicate recombination events. Interestingly, the unique large deletion found in CHA strain (15,126 bp) turned out to be the induced ssDNA filamentous Pf1-like phage.Summarising the analyses one might say that phages from different taxonomy groups recognising the same receptor generally cause the emergence of a similar type of resistance within a particular strain. However, the defence response and genome changes correlated with the receptor specificity of infecting phage differ in a strain-dependent manner.Do different strains of P. aeruginosa react similarly to a specific phage infection?The next step aimed to assess the impact of gaining phage resistance in terms of population growth efficiency as an indicator for bacterial fitness. Three of the examined wild-type strains (PAO1, A5803, and CHA) have a naturally rapid growth rate, while the other two (AA43 and PAK) display moderate growth rates. For this reason, the final results are expressed as the cumulated OD600 (Fig. 4, Table S2, S3). Overall, the majority of PAO1 mutants (61/95; 64%, p  0.001) for the clones resistant to 6–7 phages but only in the PAO1 group. Moreover, only the selection done by phage cocktails gave a statistically significant reduction of bacterial growth (p  > 0.001), while no differences were observed regarding groups treated with single LPS- or T4P- dependent phages. In contrast to the PAO1 reference strain, the statistical analyses conducted in the A5803, AA43, CHA, and PAK strains did not show any differences in terms of phage-typing profile nor phage-type selection pressure versus the population fitness reduction (growth rate).Fig. 4: The population growth efficiency as an indicator for bacterial fitness expressed as the cumulative OD600 values of 18 h culture at 37 °C measured at 20-minute intervals.Dots represent the growth of particular clones: the wild-type and control clones (green dots); mutants selected by LPS-dependent phage (red dots); mutants selected by T4P-dependent phage (blue dots); mutants selected by LPS/T4P-dependent PA5oct phage (orange dots); mutants selected by phage cocktail (black dots). * statistically different cumulative OD value compared to phage-untreated pool (p  More

1.Reynolds, M. P. et al. Improving global integration of crop research. Science 357, 359–360 (2017).CAS 
Article 

Google Scholar 
2.Braun, H., Atlin, G. & Payne, T. Multi-location testing as a tool to identify plant response to global climate change. in Climate Change and Crop Production (ed. Reynolds, M. P.) 115–138 (CABI, 2010).3.Ray, D. K., Ramankutty, N., Mueller, N. D., West, P. C. & Foley, J. A. Recent patterns of crop yield growth and stagnation. Nat. Commun. 3, 1293 (2012).Article 

Google Scholar 
4.Liu, B. et al. Similar estimates of temperature impacts on global wheat yield by three independent methods. Nat. Clim. Change 6, 1130–1136 (2016).Article 

Google Scholar 
5.Lobell, D. B. & Field, C. B. Global scale climate–crop yield relationships and the impacts of recent warming. Environ. Res. Lett. 2, 14–21 (2007).Article 

Google Scholar 
6.Asseng, S. et al. Rising temperatures reduce global wheat production. Nat. Clim. Change 5, 143–147 (2015).Article 

Google Scholar 
7.Crespo-Herrera, L. A. et al. Genetic yield gains in CIMMYT’s international elite spring wheat yield trials by modeling the genotype × environment interaction. Crop Sci. 57, 789–801 (2017).Article 

Google Scholar 
8.Tester, M. & Langridge, P. Breeding technologies to increase crop production in a changing world. Science 327, 818–822 (2010).CAS 
Article 

Google Scholar 
9.Rosegrant, M. W. & Cline, S. A. Global food security: challenges and polices. Science 302, 1917–1919 (2003).CAS 
Article 

Google Scholar 
10.Li, Y., Suontama, M., Burdon, R. D. & Dungey, H. S. Genotype by environment interactions in forest tree breeding: review of methodology and perspective on research and application. Tree Genet. Genomes 13, 60 (2017).Article 

Google Scholar 
11.Mishra, R. M. et al. Crossover interactions for grain yield in multienvironmental trials of winter wheat. Crop Sci. 46, 1291–1298 (2006).Article 

Google Scholar 
12.Allard, R. W. & Bradshaw, A. D. Implications of genotype–environmental interactions in applied plant breeding. Crop Sci. 4, 503–508 (1964).Article 

Google Scholar 
13.Reynolds, M. P., Hays, D. & Chapman, S. Breeding for adaptation to heat and drought stress. in Climate Change and Crop Production (ed. Reynolds, M. P.) 71–91 (CABI, 2010).14.Leon, N., Jannink, J., Edwards, J. W. & Kaeppler, S. M. Introduction to a special issue on genotype by environment interaction. Crop Sci. 56, 2081–2089 (2016).Article 

Google Scholar 
15.Reynolds, M. & Langridge, P. Physiological breeding. Curr. Opin. Plant Biol. 31, 162–171 (2016).Article 

Google Scholar 
16.Gourdji, S. M., Mathews, K. L., Reynolds, M., Crossa, J. & Lobell, D. B. An assessment of wheat yield sensitivity and breeding gains in hot environments. Proc. R. Soc. B. 2018, 20122190 (2012).
Google Scholar 
17.Pingali, P. L. Green revolution: impacts, limits, and the path ahead. Proc. Natl Acad. Sci. USA 109, 12302–12308 (2012).CAS 
Article 

Google Scholar 
18.Sharma, R. C. et al. Genetic gains for grain yield in CIMMYT spring bred wheat across international environment. Crop Sci. 52, 1522–1533 (2012).Article 

Google Scholar 
19.Boehm Jr, J. D., Ibba, M., Kiszonas, A. & Morris, C. F. End-use quality of CIMMYT-derived soft kernel durum wheat germplasm. II. Dough strength and pan bread quality. Crop Sci. 57, 1485–1498 (2017).Article 

Google Scholar 
20.Lillemo, M., van Ginkel, M., Trethowan, R. M., Hernandez, E. & Crossa, J. Differential adaptation of CIMMYT bread wheat to global high temperature environments. Crop Sci. 45, 2443–2453 (2005).Article 

Google Scholar 
21.Manes, Y. et al. Genetic yield gains of the CIMMYT international semi-arid wheat yield trials from 1994 to 2010. Crop Sci. 52, 1543–1552 (2012).Article 

Google Scholar 
22.You, L. et al. Spatial Production Allocation Model (SPAM) 2005 V3.2 International Food Policy Research Institute (IFPRI), International Institute fo Applied Systems Analysis (IIASA) (2017).23.Finlay, K. W. & Wilkinson, G. N. The analysis of adaptation in a plant-breeding programme. Aus. J. Agric. Res. 14, 742–754 (1963).Article 

Google Scholar 
24.De los Campos et al. A data-driven simulation platform to predict cultivars’ performance under uncertain weather conditions. Nat. Commun. 11, 4876 (2020).Article 

Google Scholar 
25.Lantican, M. A. et al. Impacts of International Wheat Improvement Research 1994–2014 (CIMMYT, 2016).26.Dreccer, M. F., Bonnett, D. & Lafarge, T. Plant breeding under a changing climate. in Encyclopedia of Sustainability Science and Technology (ed. Meyers, R. A.) 8013–8024 (Springer, 2012).27.Laiding, F., Drobek, T. & Meyer, U. Genotypic and environmental variability of yield for cultivars from 30 different crops in German official variety trials. Plant Breed. 127, 541–547 (2008).Article 

Google Scholar 
28.Allard, R. W. Principles of Plant Breeding 2nd edn (John Wiley & Sons, 1999).29.Kusmec, A., Srinivasan, S., Nettleton, D. & Schnable, P. S. Distinct genetic architectures for phenotype means and plasticities in Zea mays. Nat. Plants 3, 715–723 (2017).CAS 
Article 

Google Scholar 
30.Gauch, H. G. Statistical Analysis of Regional Yield Trials: AMMI Analysis of Factorial Designs (Elsevier, 1992). More

1.Lapola, D. M. et al. Proc. Natl Acad. Sci. USA 115, 11671–11679 (2018).CAS 
Article 

Google Scholar 
2.Phillips, O. L., Brienen, R. J. W. & the RAINFOR collaboration. Carbon Balance Manag. 12, 1 (2017)..3.Hubau, W. et al. Nature 579, 80–87 (2020).CAS 
Article 

Google Scholar 
4.Fleischer, K. et al. Nat. Geosci. 12, 736–741 (2019).CAS 
Article 

Google Scholar 
5.Huntingford, C. et al. Nat. Geosci. 6, 268–273 (2013).CAS 
Article 

Google Scholar 
6.Koch, A., Hubau, W. & Lewis, S. L. Earth’s Future 9, e2020EF001874 (2021).CAS 
Article 

Google Scholar 
7.Eyring, V. et al. Geosci. Model Dev. 9, 1937–1958 (2016).Article 

Google Scholar 
8.Climate Change and Land: an IPCC Special Report on Climate Change, Desertification, Land Degradation, Sustainable Land Management, Food Security, and Greenhouse Gas Fluxes in Terrestrial Ecosystems. (IPCC, 2019).9.Friend, A. D. et al. Proc. Natl Acad. Sci. USA 111, 3280–3285 (2014).CAS 
Article 

Google Scholar 
10.Pugh, T. A. M. et al. Biogeosciences https://doi.org/10.5194/bg-17-3961-2020 (2020).Article 

Google Scholar 
11.Hartmann, H., Adams, H. D., Anderegg, W. R. L., Jansen, S. & Zeppel, M. J. B. New Phytol. 205, 965–969 (2015).Article 

Google Scholar 
12.Scheiter, S., Langan, L. & Higgins, S. I. New Phytol. 198, 957–969 (2013).Article 

Google Scholar  More

Species-level dataSpecies range maps were derived from BirdLife International and NatureServe50 and the IUCN51. The threat data were from the IUCN Threats Classification Scheme (Version 3.2), which contains 11 primary threat classes and almost 50 subclasses52. In Red List assessments, assessors assign those threats that impact the species. For birds, the scope of the impact is also recorded categorically as the percentage of the species population that the threat impacts (unknown, negligible, 90%) and the severity, describing the scale of the impact on population declines: unknown, no decline, negligible declines, fluctuations, slow but significant declines (30%).Model development approachWe designed our analytical framework with three considerations in mind. First, the threat location information is limited: for each species, the data only describe whether a species is threatened by a given activity anywhere within its range (data on the timing, scope and severity of threats are available only for birds and are not spatially explicit). Second, we wanted to compare the spatial patterns of threat against independent data on spatial distributions of human activities. Third, for many activities, the relationship between human activity (for example, hunting or invasive species and diseases) and biodiversity response is poorly understood. We therefore chose not to incorporate known patterns of human activity as explanatory variables in our models.In the absence of global datasets on the spatial patterns of the impact probability of each threat, we used a simulation approach to develop our models and assess the ability of different model parameterizations to reproduce our simulated threat. This process had four steps (Extended Data Fig. 1).Simulated threat intensity mapsFirst, we simulated a continuous synthetic threat across sub-Saharan Africa. The concept behind this is that a credible model should be able to reproduce a ‘true’, synthetic threat pattern on the basis of information comparable to that available in the Red List. To test this, we generated a set of synthetic, continuous surfaces of threat intensity with different levels of spatial autocorrelation and random variation (Supplementary Fig. 1). This was achieved by taking a grid of 50 km × 50 km (2,500 km2) pixels across the Afrotropic biogeographic realm (i.e., sub-Saharan Africa). Threat intensity was modelled as a vector of random variables, Z, one for each pixel i, generated with a correlation structure given by the distance matrix between points weighted by a scalar value, r, indicating the degree of correlation (equations (1–3)). Four values of r were used: 1 × 10−6, which yields very strong autocorrelation; 1 × 10−4, which yields strong autocorrelation; 0.05, which yields moderate autocorrelation; and 0.3, which produces a low-correlation, localized pattern (Supplementary Fig. 1). The model included the following equations:$${mathbf{Z}}(r) = U^{mathrm{T}}{mathrm{Norm}}left( {n,0,1} right)$$
(1)
$$W = UU^{ast}$$
(2)
$$W = {mathrm{e}}^{left( { – rD} right)}$$
(3)
where r is a scalar determining the degree of spatial autocorrelation (as r decreases, the autocorrelation increases), D is the Euclidean distance matrix between each pair of pixels, W is the matrix of weights for the threat intensity, U and U* are the upper triangular factors of the Choleski decomposition of W and its conjugate transpose, UT is the transpose of U and n is the number of pixels.We chose the Afrotropic biogeographic realm (sub-Saharan Africa) as our geography within which to develop the modelling approach because it permitted more rapid iterations than a global-scale simulation while also retaining characteristics of importance for the model evaluation such as strong environmental gradients and heterogeneity in species richness. However, for the simulation, no information from the geography or overlapping species ranges was used, except the spatial configuration of the polygons. Thus, the use of the Afrotropic realm was purely to avoid generating thousands of complex geometries for the purpose of the simulation. Using a real geography and actual species ranges ensures that our simulation contains conditions that are observed in reality (for example, areas of high and low species richness also observed in the real world). We took the simulated threat maps generated through this process to be our ‘true’ likelihood of a randomly drawn species that occurs in that location being impacted by the synthetic threat (Supplementary Fig. 1).Simulating the red-listing processSecond, we wanted to simulate the red-listing process whereby experts evaluate whether a threat is impacting a species on the basis of the overall threat intensity within its range. For this, we used the range maps for all mammal species in Africa and assigned a binary threat classification (that is, affected or not affected) to each species on the basis of the values of the synthetic threat within each species’ range. We assumed that the binary assessment of threat for a species is based on whether the level of impact across a proportion of its range is judged as significant. This step was intended to replicate the real red-listing process, where assessors define threats that impact the species on the basis of an assessment of the information available on threatening mechanisms and species responses. In practice, this was done by overlaying the real range maps for mammals over the four simulated threat surfaces and assessing the intensity of synthetic threat within each species range map. We wanted to assign species impacts considering that species will be more likely to be impacted if a greater part of their range has a high threat intensity. Understanding how to set a threshold for what intensity would constitute sufficient threat to be assessed as affected is a complicated exercise. We thus tested three thresholds to capture different assumptions. These thresholds were chosen after discussion with leading experts on the red-listing process. More specifically, we calculated the 25th, 50th and 75th percentiles of threat intensity across pixels within the species range. We then used a stochastic test to convert these quantiles to binary threat class, C. For each species, we produced a set of ten draws from a uniform distribution bounded by 0 and 1. If over half of the draws were lower than the threat intensity quantile, the species was classified as threatened for that percentile.The above simulation assumes perfect knowledge of the threat intensities across the species range, which might not always be the case in the actual red-listing process. In real life, certain areas within species ranges are less well known for a suite of different reasons. To incorporate some uncertainty about the knowledge of the red-listing experts about the ‘true’ threat intensity, we constructed a layer to describe the spatial data uncertainty associated with the Red List. This aspect was intended to simulate the imperfect knowledge of the simulated ‘Red List assessors’. This layer was calculated as the proportion of species present in a given location that are categorized as Data Deficient—in other words, there is insufficient information known about the species to assess its extinction risk using the IUCN Red List Criteria (Extended Data Fig. 7). Then, when calculating the 25th, 50th and 75th percentiles of threat intensity across each range, we weighted this calculation by one minus the proportion of Data Deficient species, so that more uncertain places (those with a greater proportion of Data Deficient species) contributed less to the calculation than locations where knowledge was more certain. These were then converted to a binary threat class accounting for uncertainty in expert knowledge among the simulated ‘assessors’, CUncertain, using the same stochastic process described above for the calculation of C.This step produced, for each species, a threat classification analogous to the threat classification assigned by experts as part of the IUCN Red List process. Six sets of threat classifications were produced for each synthetic threat surface, on the basis of the 25th, 50th and 75th percentiles with perfect (C0.25, C0.5 and C0.75) or uncertain (CUncertain-0.25, CUncertain-0.5 and CUncertain-0.75) spatial knowledge.Model formulation and selectionThird, using all species polygons with assigned threat assessments from step 2 (that is, affected or not affected), we fitted nine candidate models and predicted the estimated probability of impact for each grid cell. Then, in a fourth step, we compared the predicted probabilities of impact produced in step 3 with the original synthetic threat maps created in step 1 to test the predictive ability of our models.The Red List threat assessment does not contain information on where in the range the impact occurs. Therefore, a species with a very small range provides higher spatial precision about the location of the impact, whereas a species with a large range may be impacted anywhere within a wide region. To address this lack of precision in the impact location, we took the area of each species range to serve as a proxy for the spatial certainty of the impact information. The certainty that a species was impacted or not impacted in a given cell depended on its range size, R. The models we evaluated therefore incorporated R in different ways (Supplementary Table 1).The models were fitted as a binomial regression with a logit link function. For each pixel, the model predicts the probability of impact, PTh—in other words, the probability that if you sampled a species at random from those that occur in that pixel, the species would be impacted by the activity being considered. To account for uncertainties in the simulation of the threat assessment process (thresholds for impact and perfect or imperfect knowledge), models were fitted to the six sets of threat codes (C0.25, C0.5, C0.75, CUncertain-0.25, CUncertain-0.5 and CUncertain-0.75), and the root mean squared error (RMSE) was calculated between PTh and the simulated threat intensity, Z(r), for each value of r. For each simulation, we ranked the different models according to their model fit as measured by the RMSE. We assessed these ranks across all simulations and sets of threat codes. We evaluated the models on the basis of the ranks of RMSE, across the threat code sets and threat intensity maps. Rank distributions for each model are shown in Extended Data Fig. 2, and the results from these models are shown in Supplementary Tables 1 and 2.All models were correlated (Pearson’s r2  > 0.5), albeit with some variation between model types and across the simulation parameters (Supplementary Fig. 2). However, some models had greater predictive accuracy when evaluated using the RMSE. The top four ranking models were, in order of decreasing summed rank, (1) inverse of cube root of range size as a weight, (2) inverse 2.5 root of range size as a weight, (3) inverse square root of range size as a weight and (4) inverse natural logarithm of range size as a weight. The fact that these four models showed good model fit suggests that the best model structure had a measure of range size as a weight but that the model was not particularly sensitive to the transformation of range size.The best-fitting model across the range of simulation parameters was an intercept-only logistic regression where the response variable was the binary threat code (1 = threatened, 0 = not threatened) for each species in the pixel and where the inverse cube root of the range size of each species was used as a weight. The model was concordant across the set of simulated datasets with a relationship that was predominantly linear with r2 between 0.47 and 0.7, depending on simulation parameters for Z(r) in 0.05, 10−4 and 10−6, centred around unity and with the RMSE ranging between 0.129 and 0.337 depending on simulation parameters (Supplementary Figs. 2 and 3). The choice of the inverse cube root range size weight was based on the performance of this against eight other model types (Supplementary Fig. 4 and Supplementary Table 1).We conducted a decomposition of variance in model performance using a binomial regression model, with RMSE as the dependent variable and model type, knowledge level and autocorrelation structure as the independent factorial variables. This showed that knowledge about the threats underlying each species range and how that threat information is used in the assessment explained the vast majority (94.7%) of the variation in RMSE outcomes (Supplementary Fig. 4).For birds, further information on the scope of the threat was available as an ordinal variable describing the fraction of range that the threat covers. We explored the use of scope in our models but concluded that to avoid arbitrary decisions about the scope of non-threatened species (where they are either not threatened anywhere or threatened in only a small part of their range), and for consistency with other taxonomic groups, we would model birds using the same model structure as used for mammals and amphibians (see the Supplementary Methods for further details).Mapping probability of impactOnce the best-performing model was identified using the simulated data, we then used this model on the actual Red List threat and range data to develop threat maps. This model produced threat maps for each taxonomic group (amphibians, birds and mammals) of the probability of impact, PTh, for each individual threat. For a given pixel, threat and taxonomic group, this estimates the probability that a randomly sampled species with a range overlapping with that pixel is being impacted by the threat, while taking into account spatial imprecision in the Red List data.Threat maps were generated using range map data and threat assessments from the IUCN Red List18. We intersected range maps for 22,898 extant terrestrial amphibians (n = 6,458), birds (n = 10,928; excluding the spatial areas within the range that are associated with ‘Passage’—where the species is known or thought very likely to occur regularly during relatively short periods of the year on migration) and mammals (n = 5,512; including those with uncertain ranges) with a global 50 km × 50 km (2,500 km2) resolution, equal-area grid for the terrestrial world. This provided, for each 50 km × 50 km pixel, a list of the species whose range overlapped it, along with the associated range size of each species. For each pixel and taxonomic group (amphibians, birds and mammals) independently, we then modelled the probability of impact, PTh,Activity (for example, PTh,Logging for logging, PTh,Agriculture for agriculture or PTh,Pollution for pollution), for each of the six threats: agriculture, hunting and trapping, logging, pollution, invasive species and diseases, and climate change. We focused on these as the six main threats as defined by the Intergovernmental Science-Policy Platform on Biodiversity and Ecosystem Services4, but our methodological framework is flexible and could be expanded to other threats in the IUCN classification19. We used only taxonomic groups with a sufficiently high total number of species and where they have been comprehensively assessed so that potential biases associated with the groups of species prioritized by experts are avoided.Calculating uncertainties for the threat probabilityWe estimated a measure of uncertainty associated with our impact probability predictions using maps of the proportions of Data Deficient species in each cell within each taxonomic class (amphibians, birds or mammals) as a measure of knowledge certainty in that cell. The rationale for this approach is that places with more Data Deficient species with unknown threatened status should have greater uncertainty in the probability of impact. We therefore created greater variation in the data where there were more Data Deficient species. We used the knowledge-certainty map to probabilistically draw a sample of 100 threat codes for each species, on the basis of the median Data Deficiency across the species range. The random sample changed the species threat code with a probability related to the proportion of Data Deficient species within its range. If the median proportion of Data Deficient species was zero, then we assumed that there was a small probability (0.005) that the species could have been incorrectly coded. Where the median proportion was greater than zero, the probability increased linearly. So, for a species with 5% Data Deficient species within its range, the sample changed the species threat code with a probability close to 5%; if the median proportion was equal to 0.5, then the probability of the species being incorrectly assigned was equal to 0.5. We then fitted the impact probability model with each of the 100 species threat codes and generated a distribution of predicted threat probabilities in each grid cell, from which we took the 95% confidence intervals as the uncertainty estimates (Extended Data Figs. 8–10).Evaluating modelled threat patternsWe evaluated the spatial patterns of threat on the basis of the real Red List threat assessment data against empirical data in two independent ways. First, we compared the probability of impact from logging and agriculture combined within forested biomes (that is, corresponding to remotely detected forest loss, which we refer to as the probability of impact from forest loss, PTh,Forest-loss) with data on forest cover change10. Forest cover change was aggregated from their native 30 m × 30 m (900 m2) resolution pixels to our 50 km × 50 km resolution pixels using Google Earth Engine. For each 50 km × 50 km pixel, we calculated the total area lost between 2000 and 2013 and the area lost as a proportion of the area in 2000. We restricted our analysis to forested biomes: (1) tropical and subtropical moist broadleaf forests, (2) tropical and subtropical dry broadleaf forests, (3) tropical and subtropical coniferous forests, (4) temperate broadleaf and mixed forests, (5) temperate coniferous forests and (6) boreal forests/taiga, following the World Wildlife Fund’s ecoregions classification53. The relationship between forest loss and the probability of impact from forest loss as captured by agriculture and logging overall showed a significant positive correlation: PTh,Forest-loss increased with increasing forest cover loss (P  More

1.Collins, J. P., & Crump, M. L. Extinction in Our Times: Global Amphibian Decline. (2009).2.Catenazzi, A. State of the world’s amphibians. Annu. Rev. Environ. Resour. 40, 91–119 (2015).Article 

Google Scholar 
3.González-del-Pliego, P. et al. Phylogenetic and trait-based prediction of extinction risk for data deficient amphibians. Curr. Biol. 29, 1557–1563 (2019).PubMed 
Article 
CAS 
PubMed Central 

Google Scholar 
4.Lips, K. R. et al. Emerging infectious disease and the loss of biodiversity in a Neotropical amphibian community. PNAS 102, 3165–3170 (2006).ADS 
Article 
CAS 

Google Scholar 
5.Lips, K. R., Diffendorfer, J., Mendelson, J. R. & Sears, M. W. Riding the wave: Reconciling the roles of disease and climate change in amphibian declines. PLoS Biol. 6, e72 (2008).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
6.Berger, L. et al. Chytridiomycosis causes amphibian mortality associated with population declines in the rain forests of Australia and Central America. PNAS 95, 9031–9036 (1998).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
7.Longcore, J. E., Pessier, A. P. & Nichols, D. K. Batrachochytrium dendrobatidis gen. et sp. nov., a chytrid pathogenic to amphibians. Mycologia 91, 219–227 (1999).Article 

Google Scholar 
8.Berger, L. et al. History and recent progress on chytridiomycosis in amphibians. Fungal Ecol. 19, 89–99 (2016).Article 

Google Scholar 
9.Martel, A. et al. Batrachochytrium salamandrivorans sp. nov. causes lethal chytridiomycosis in amphibians. PNAS 110, 15325–15329 (2013).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
10.Scheele, B. C. et al. Amphibian fungal panzootic causes catastrophic and ongoing loss of biodiversity. Science 363, 1459–1463 (2019).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
11.Lambert, M. R. et al. Comment on “Amphibian fungal panzootic causes catastrophic and ongoing loss of biodiversity”. Science https://doi.org/10.1126/science.aay1838 (2020).Article 
PubMed 
PubMed Central 

Google Scholar 
12.Scheele, B. C. et al. Response to Comment on “Amphibian fungal panzootic causes catastrophic and ongoing loss of biodiversity”. Science https://doi.org/10.1126/science.aay2905 (2020).Article 
PubMed 
PubMed Central 

Google Scholar 
13.Puschendorf, R. et al. Distribution models for the amphibian chytrid Batrachochytrium dendrobatidis in Costa Rica: Proposing climatic refuges as a conservation tool. Divers. Distrib. 15, 401–408 (2009).Article 

Google Scholar 
14.Zumbado-Ulate, H. et al. Endemic infection of Batrachochytrium dendrobatidis in Costa Rica: Implications for amphibian conservation at regional and species level. Diversity 11, 129 (2019).Article 

Google Scholar 
15.Crawford, A. J., Lips, K. R. & Bermingham, E. Epidemic disease decimates amphibian abundance, species diversity, and evolutionary history in the highlands of central Panama. PNAS 107, 13777–13782 (2010).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
16.Woodhams, D. C. et al. Chytridiomycosis and amphibian population declines continue to spread eastward in Panama. EcoHealth 5, 268–274 (2008).PubMed 
Article 
PubMed Central 

Google Scholar 
17.Catenazzi, A., Lehr, E., Rodriguez, L. & Vredenburg, V. Batrachochytrium dendrobatidis and the collapse of anuran species richness and abundance in the upper Manu National Park, southeastern Peru. Conserv. Biol. 25, 382–391 (2011).PubMed 
Article 
PubMed Central 

Google Scholar 
18.Voyles, J. et al. Pathogenesis of chytridiomycosis, a cause of catastrophic amphibian declines. Science 326, 582–585 (2009).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
19.James, T. et al. Disentangling host, pathogen, and environmental determinants of a recently emerged wildlife disease: Lessons from the first 15 years of amphibian chytridiomycosis research. Ecol. Evol. 5, 4079–4097 (2015).PubMed 
PubMed Central 
Article 

Google Scholar 
20.Soto-Azat, C. et al. Xenopus laevis and emerging amphibian pathogens in Chile. EcoHealth 13, 775–783 (2016).PubMed 
Article 
PubMed Central 

Google Scholar 
21.Ron, S. Predicting the distribution of the amphibian pathogen Batrachochytrium dendrobatidis in the new world. Biotropica 37, 209–221 (2005).Article 

Google Scholar 
22.Rödder, D., Kielgast, J. & Lötters, S. Future potential distribution of the emerging amphibian chytrid fungus under anthropogenic climate change. Dis. Aquat. Org. 92, 201–207 (2010).Article 

Google Scholar 
23.Murray, K. A. et al. Assessing spatial patterns of disease risk to biodiversity: Implications for the management of the amphibian pathogen, Batrachochytrium dendrobatidis. J. Appl. Ecol. 48, 163–173 (2011).Article 

Google Scholar 
24.Liu, X., Rohr, J. & Li, Y. Climate, vegetation, introduced hosts and trade shape a global wildlife pandemic. Proc. Biol. Sci. 280, 20122506 (2013).PubMed 
PubMed Central 

Google Scholar 
25.Olson, D. H. et al. Mapping the global emergence of Batrachochytrium dendrobatidis, the amphibian chytrid fungus. PLoS ONE 8, e56802 (2013).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
26.Penner, J. et al. West Africa—A safe haven for frogs? A sub-continental assessment of the chytrid fungus (Batrachochytrium dendrobatidis). PLoS ONE 8, e56236 (2013).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
27.Xie, G. Y., Olson, D. H. & Blaustein, A. R. Projecting the global distribution of the emerging amphibian fungal pathogen, Batrachochytrium dendrobatidis, based on IPCC climate futures. PLoS ONE 11, e0160746 (2016).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
28.Harris, R. N. et al. Skin microbes on frogs prevent morbidity and mortality caused by a lethal skin fungus. ISME J. 3, 818–824 (2009).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
29.Searle, C. L. et al. Differential host susceptibility to Batrachochytrium dendrobatidis, an emerging amphibian pathogen. Conserv. Biol. 25, 965–974 (2011).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
30.Farrer, R. et al. Multiple emergences of genetically diverse amphibian-infecting chytrids include a globalized hypervirulent recombinant lineage. PNAS 108, 18732–18736 (2011).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
31.Lips, K. R. Overview of chytrid emergence and impacts on amphibians. Philos. Trans. R. Soc. B. 371, 20150465 (2016).Article 

Google Scholar 
32.Bolom-Huet, R., Pineda, E., Díaz-Fleischer, F., Muñoz-Alonso, A. L. & Galindo-González, J. Known and estimated distribution in Mexico of Batrachochytrium dendrobatidis, a pathogenic fungus of amphibians. Biotropica 51, 731–746 (2019).Article 

Google Scholar 
33.Zumbado-Ulate, H., García-Rodríguez, A. & Searle, C. L. Species distribution models predict the geographic expansion of an enzootic amphibian pathogen. Biotropica 53, 221–231 (2021).Article 

Google Scholar 
34.Berger, L. et al. Effect of season and temperature on mortality in amphibians due to chytridiomycosis. Aust. Vet. J. 82, 434–439 (2004).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
35.Clare, F. C. et al. Climate forcing of an emerging pathogenic fungus across a montane multi-host community. Philos. Trans. R. Soc. B. 371, 20150454 (2016).Article 
CAS 

Google Scholar 
36.Bacigalupe, L. D., Soto-Azat, C., García-Vera, C., Barría-Oyarzo, I. & Rezende, E. L. Effects of amphibian phylogeny, climate and human impact on the occurrence of the amphibian-killing chytrid fungus. Glob. Change Biol. 23, 3543–3553 (2017).ADS 
Article 

Google Scholar 
37.Raffel, T., Michel, P., Sites, W. & Rohr, J. What drives chytrid infections in newt populations? Associations with substrate, temperature, and shade. EcoHealth 7, 526–536 (2010).PubMed 
Article 
PubMed Central 

Google Scholar 
38.Pounds, J. A. et al. Widespread amphibian extinctions from epidemic disease driven by global warming. Nature 439, 161–167 (2006).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
39.Hudson, M. et al. Reservoir frogs: Seasonality of Batrachochytrium dendrobatidis infection in robber frogs. PeerJ 7, e7021 (2019).PubMed 
PubMed Central 
Article 

Google Scholar 
40.Kriger, K. M. & Hero, J. M. Large-scale seasonal variation in the prevalence and severity of chytridiomycosis. J. Zool. 271, 352–359 (2007).
Google Scholar 
41.Longo, A. V., Burrowes, P. A. & Joglar, R. L. Seasonality of Batrachochytrium dendrobatidis infection in direct-developing frogs suggests a mechanism for persistence. Dis. Aquat. Org. 92, 253–260 (2010).Article 

Google Scholar 
42.Zumbado-Ulate, H., Bolaños, F., Gutiérrez-Espeleta, G. & Puschendorf, R. Extremely low prevalence of Batrachochytrium dendrobatidis in frog populations from Neotropical dry forest of Costa Rica supports the existence of a climatic refuge from disease. EcoHealth 11, 593–602 (2014).PubMed 
Article 
PubMed Central 

Google Scholar 
43.Bacigalupe, L. D. et al. The amphibian-killing fungus in a biodiversity hotspot: Identifying and validating high-risk areas and refugia. Ecosphere. 10, e02724 (2019).Article 

Google Scholar 
44.Flechas, S. V. et al. Current and predicted distribution of the pathogenic fungus Batrachochytrium dendrobatidis in Colombia, a hotspot of amphibian biodiversity. Biotropica 49, 685–694 (2017).Article 

Google Scholar 
45.Lampo, M. et al. Batrachochytrium dendrobatidis in Venezuela. Herpetol. Rev. 39, 449 (2008).
Google Scholar 
46.Valenzuela-Sánchez, A. et al. Genomic epidemiology of the emerging pathogen Batrachochytrium dendrobatidis from native and invasive amphibian species in Chile. Transbound. Emerg. Dis. 65, 309–314 (2018).PubMed 
Article 
PubMed Central 

Google Scholar 
47.O’Hanlon, S. et al. Recent Asian origin of chytrid fungi causing global amphibian declines. Science 360, 621–627 (2018).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
48.Soto-Azat, C. et al. The population decline and extinction of Darwin’s frogs. PLoS ONE 8, e66957 (2013).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
49.Soto-Azat, C. et al. ASG Chile leads update of the extinction risk of Chilean amphibians for the IUCN red list of threatened speciesTM. FrogLog 23, 6–7 (2015).
Google Scholar 
50.Mora, M. et al. High abundance of invasive African clawed frog Xenopus laevis in Chile: Challenges for their control and updated invasive distribution. Manag. Biol. Invasions. 10, 377–388 (2019).Article 

Google Scholar 
51.Solís, R., Penna, M., De la Riva, I., Fisher, M. & Bosch, J. Presence of Batrachochytrium dendrobatdis in anurans from the Andes highlands of northern Chile. Herpetol. J. 24, 55–59 (2015).
Google Scholar 
52.Soto-Azat, C. et al. Is Chytridiomycosis driving Darwin’s frogs to extinction?. PLoS ONE 8, e79862 (2013).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
53.Valenzuela-Sánchez, A. et al. Cryptic disease-induced mortality may cause host extinction in an apparently stable host–parasite system. Proc. Biol. Sci. 284, 20171176 (2017).PubMed 
PubMed Central 

Google Scholar 
54.Lips, K. R., Reeve, J. D. & Witters, L. R. Ecological traits predicting amphibian population declines in Central America. Conserv. Biol. 17, 1078–1088 (2003).Article 

Google Scholar 
55.Hero, J. M., Williams, S. E. & Magnusson, W. E. Ecological traits of declining amphibians in upland areas of eastern Australia. J. Zool. 267(3), 221–232 (2005).Article 

Google Scholar 
56.Kriger, K. M. & Hero, J. M. Altitudinal distribution of chytrid (Batrachochytrium dendrobatidis) infection in subtropical Australian frogs. Austral. Ecol. 33(8), 1022–1032 (2008).Article 

Google Scholar 
57.Skerratt, L. F. et al. Spread of chytridiomycosis has caused the rapid global decline and extinction of frogs. EcoHealth 4, 125 (2007).Article 

Google Scholar 
58.Langwig, K. et al. Context-dependent conservation responses to emerging wildlife diseases. Front. Ecol. Environ. 13, 195–202 (2015).Article 

Google Scholar 
59.Shaw, S. D. et al. The distribution and host range of Batrachochytrium dendrobatidis in New Zealand, 1930–2010. Ecology 94, 2108–2111 (2013).Article 

Google Scholar 
60.Ghirardi, R. et al. Endangered amphibians infected with the chytrid fungus Batrachochytrium dendrobatidis in austral temperate wetlands from Argentina. Herpetol. J. 24, 129–133 (2014).
Google Scholar 
61.Bielby, J., Cooper, N., Cunningham, A., Garner, T. & Purvis, A. Predicting susceptibility to future declines in the world’s frogs. Conserv. Lett. 1, 82–90 (2008).Article 

Google Scholar 
62.Barrionuevo, S. & Mangione, S. Chytridiomycosis in two species of Telmatobius (Anura: Leptodactylidae) from Argentina. Dis. Aquat. Org. 73, 171–174 (2006).Article 

Google Scholar 
63.Seimon, T. A. et al. Upward range extension of Andean anurans and chytridiomycosis to extreme elevations in response to tropical deglaciation. Glob. Change Biol. 13, 288–299 (2007).ADS 
Article 

Google Scholar 
64.Burrowes, P. A. & De la Riva, I. Unraveling the historical prevalence of the invasive chytrid fungus in the Bolivian Andes: Implications in recent amphibian declines. Biol. Invasions. 19, 1781–1794 (2017).Article 

Google Scholar 
65.Vredenburg, V. T., Knapp, R., Tunstall, T. & Briggs, C. Dynamics of an emerging disease drive large-scale amphibian population extinctions. PNAS 107, 9689–9694 (2010).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
66.Azat, C. et al. A flagship for Austral temperate forest conservation: an action plan for Darwin’s frogs bringing together key stakeholders. Oryx 55, 356–363 (2021). Article 

Google Scholar 
67.Pilliod, D. S. et al. Effects of amphibian chytrid fungus on individual survival probability in wild boreal toads. Conserv. Biol. 24, 1259–1267 (2010).PubMed 
Article 
PubMed Central 

Google Scholar 
68.Walker, S. F. et al. Factors driving pathogenicity vs. prevalence of amphibian panzootic chytridiomycosis in Iberia. Ecol. Lett. 13, 372–382 (2010).PubMed 
Article 
PubMed Central 

Google Scholar 
69.Kriger, K. M., Pereoglou, F. & Hero, J. M. Latitudinal variation in the prevalence and intensity of chytrid (Batrachochytrium dendrobatidis) infection in eastern Australia. Conserv. Biol. 21, 1280–1290 (2007).PubMed 
Article 
PubMed Central 

Google Scholar 
70.Petersen, C. E., Lovich, R. E., Phillips, C. A., Dreslik, M. J. & Lannoo, M. J. Prevalence and seasonality of the amphibian chytrid fungus Batrachochytrium dendrobatidis along widely separated longitudes across the United States. EcoHealth 13, 368–382 (2016).PubMed 
Article 
PubMed Central 

Google Scholar 
71.Thorpe, C. J. et al. Climate structuring of Batrachochytrium dendrobatidis infection in the threatened amphibians of the northern Western Ghats, India. R. Soc. Open Sci. 5, 180211 (2018).ADS 
PubMed 
PubMed Central 
Article 

Google Scholar 
72.Sonn, J. M., Utz, R. M. & Richards-Zawacki, C. L. Effects of latitudinal, seasonal, and daily temperature variations on chytrid fungal infections in a North American frog. Ecosphere 10, e02892 (2019).Article 

Google Scholar 
73.Raffel, T. R., Halstead, N. T., McMahon, T. A., Davis, A. K. & Rohr, J. R. Temperature variability and moisture synergistically interact to exacerbate an epizootic disease. Proc. Biol. Sci. 282, 20142039 (2015).PubMed 
PubMed Central 

Google Scholar 
74.Woodhams, D. C. & Alford, R. A. Ecology of chytridiomycosis in rainforest stream frog assemblages of tropical Queensland. Conserv. Biol. 19, 1449–1459 (2005).Article 

Google Scholar 
75.Adams, M. J. et al. Using occupancy models to understand the distribution of an amphibian pathogen, Batrachochytrium dendrobatidis. Ecol. Appl. 20, 289–302 (2010).PubMed 
Article 
PubMed Central 

Google Scholar 
76.Fisher, M., Garner, T. & Walker, S. F. Global emergence of Batrachochytrium dendrobatidis and amphibian chytridiomycosis in space, time, and host. Annu. Rev. Microbiol. 63, 291–310 (2009).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
77.Schloegel, L. M. et al. Novel, panzootic and hybrid genotypes of amphibia chytridiomycosis associated with the bullfrog trade. Mol. Ecol. 21, 5162–5177 (2012).PubMed 
Article 
PubMed Central 

Google Scholar 
78.Wilson, E. A., Briggs, C. J. & Dudley, T. L. Invasive African clawed frogs in California: A reservoir for or predator against the chytrid fungus?. PLoS ONE 13, e0191537 (2018).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
79.Becker, C. G., Longo, A. V., Haddad, C. F. & Zamudio, K. R. Land cover and forest connectivity alter the interactions among host, pathogen and skin microbiome. Proc. Biol. Sci. 284, 20170582 (2017).PubMed 
PubMed Central 

Google Scholar 
80.McCoy, K. A. & Peralta, A. L. Pesticides could alter amphibian skin microbiomes and the effects of Batrachochytrium dendrobatidis. Front. Microbiol. 9, 748 (2018).PubMed 
PubMed Central 
Article 

Google Scholar 
81.Ellis, E. & Ramankutty, N. Putting people in the map: Anthropogenic biomes of the world. Front. Ecol. Environ. 6, 439–447 (2008).Article 

Google Scholar 
82.Rohr, J., Halstead, N. & Raffel, T. Modelling the future distribution of the amphibian chytrid fungus: The influence of climate and human-associated factors. J. Appl. Ecol. 48, 174–176 (2011).Article 

Google Scholar 
83.Myers, N., Mittermeier, R. A., Mittermeier, C. G., Da Fonseca, G. A. & Kent, J. Biodiversity hotspots for conservation priorities. Nature 403, 853–858 (2000).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
84.Echeverria, C., Coomes, D., Hall, M. & Newton, A. Spatially explicit models to analyze forest loss and fragmentation between 1976 and 2020 in southern Chile. Ecol. Model. 212, 439–449 (2008).Article 

Google Scholar 
85.Rodriguez, D., Becker, C., Pupin, C., Haddad, F. & Zamudio, K. Long-term endemism of two highly divergent lineages of the amphibian-killing fungus in the Atlantic Forest of Brazil. Mol. Ecol. 23, 774–787 (2014).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
86.Puschendorf, R., Hodgson, L., Alfors, R. A., Skerrat, L. F. & VanDerWal, J. Underestimated ranges and overlooked refuges from amphibian chytridiomycosis. Divers. Distrib. 19, 1313–1321 (2013).Article 

Google Scholar 
87.Scheele, B. C. et al. After the epidemic: Ongoing declines, stabilizations and recoveries in amphibians afflicted by chytridiomycosis. Biol. Conserv. 206, 37–46 (2017).Article 

Google Scholar 
88.Mendelson, J. R. III., Whitfield, S. M. & Sredl, M. J. A recovery engine strategy for amphibian conservation in the context of disease. Biol. Conserv. 236, 188–191 (2019).Article 

Google Scholar 
89.Van Rooij, P., Martel, A., Haesebrouck, F. & Pasmans, F. Amphibian chytridiomycosis: A review with focus on fungus-host interactions. Vet. Res. 46, 1–22 (2015).Article 
CAS 

Google Scholar 
90.Christie, M. R. & Searle, C. L. Evolutionary rescue in a host–pathogen system results in coexistence not clearance. Evol. Appl. 11, 681–693 (2018).PubMed 
Article 
PubMed Central 

Google Scholar 
91.Bletz, M. C. et al. Mitigating amphibian chytridiomycosis with bioaugmentation: Characteristics of effective probiotics and strategies for their selection and use. Ecol. Lett. 16, 807–820 (2013).PubMed 
Article 
PubMed Central 

Google Scholar 
92.Bosch, J. et al. Successful elimination of a lethal wildlife infectious disease in nature. Biol. Lett. 11, 20150874 (2015).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
93.Olson, D. M. et al. Terrestrial ecoregions of the world: A new map of life on earth. Bioscience 51, 933–938 (2001).Article 

Google Scholar 
94.Pellet, J. & Schmidt, B. R. Monitoring distributions using call surveys: Estimating site occupancy, detection probabilities and inferring absence. Biol. Conserv. 123, 27–35 (2005).Article 

Google Scholar 
95.Drechsler, A. & Bock, D. Ortmann’s funnel trap—A highly efficient tool for monitoring amphibian species. Herpetol. Notes. 3, 13–21 (2010).
Google Scholar 
96.Hudson, M. et al. Dynamics and genetics of a disease-driven species decline to near extinction: Lessons for conservation. Sci. Rep. 6, 1–13 (2016).Article 
CAS 

Google Scholar 
97.R Development Core Team. R: A Language and Environment for Statistical Computing, https://www.R-project.org/ (2019).98.Sanderson, E. W. et al. The human footprint and the last of the wild. Bioscience 52, 891–904 (2002).Article 

Google Scholar 
99.Center for International Earth Science Information Network (CIESIN). Gridded Population of the World, Version 3 (GPWv3). https://doi.org/10.7927/H4639MPP (2005).100.Booth, T. H., Nix, H. A., Busby, J. R. & Hutchinson, M. F. BIOCLIM: The first species distribution modelling package, its early applications and relevance to most current MAXENT studies. Divers. Distrib. 20, 1–9 (2014).Article 

Google Scholar 
101.Center for International Earth Science Information Network (CIESIN). Gridded Species Distribution: Global Amphibian Richness Grids. https://doi.org/10.7927/H4RR1W66 (2015).102.Fick, S. & Hijmans, R. WorldClim 2: New 1-km spatial resolution climate surfaces for global land areas. Int. J. Climatol. 37, 4302–4315 (2017).Article 

Google Scholar 
103.ASTER. ASTER global digital elevation model V003. https://doi.org/10.5067/ASTER/ASTGTM (2018).104.QGIS Development Team. QGIS Geographic Information System. Open Source Geospatial Foundation Project. http://qgis.osgeo.org (2018).105.Zuur, A. F., Ieno, E. N. & Elphick, C. S. A protocol for data exploration to avoid common statistical problems. Methods Ecol. Evol. 1, 3–14 (2010).Article 

Google Scholar 
106.Fox, J. Effect displays in R for generalised linear models. J. Stat. Softw. 8, 1–27 (2003).Article 

Google Scholar 
107.Carpenter, T. E. Methods to investigate spatial and temporal clustering in veterinary epidemiology. Prev. Vet. Med. 48, 303–320 (2001).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
108.Kulldorff, M. A spatial scan statistic. Commun. Stat-Theor. M. 26, 1481–1496 (1997).MathSciNet 
MATH 
Article 

Google Scholar 
109.Kulldorff, M. Information Management Services, Inc. SaTScanTM v.9.4.4: software for the spatial and space-time scan statistics. http://www.satscan.org (2009). More

Photosynthetic characteristicsWater and nitrogen coupling treatment had a significant effect on the photosynthetic characteristics (Fig. 1). Generally, the net photosynthetic rates of the treatments were in the following order: CK, W1N1, W1N3, W3N1, W3N3, W2N1, W1N2, W3N2, W2N3, and W2N2. The treatments with low water and low nitrogen had significantly lower net photosynthetic rates than W2N2. The stomatal conductance and transpiration rate changed in similar patterns. The net photosynthetic rate showed a unimodal trend with the increase of nitrogen application at the same irrigation level. Under the same nitrogen application level, the net photosynthetic rate increased first and then decreased slowly with the increase of irrigation amount, with the highest photosynthetic rates in the order of W2  > W3  > W1. The net photosynthetic rate was the highest, with a mean value was 13.87 μmol m−2 s−1, in treatment W2N2. The results showed that severe water stress and excessive nitrogen were not conducive to the absorption and utilization of water and nutrients by crop roots, which led to the decrease of the photosynthetic rate. The effect of water and nitrogen treatment on the intercellular CO2 concentration was significant (Fig. 1). Under the condition of excessive water or nitrogen, the photosynthesis of Isatis indigotica decreased, and the intercellular CO2 concentration showed a trend opposite to that of the net photosynthetic rate.Compared with N, P, and K deficiency treatments, water–N coupling could increase the Pn of crops, which was the same as that of other fruit trees and vegetables13. Accumulated photoassimilates in the third internode of the upper part of the main stems, as well as in the flag leaf sheath, are mobilized in a higher proportion and can contribute to grain filling in rice plants subjected to water stress in the tillering phase14. The Pn, Gs, and Tr of maize leaves at the seedling stage decreased significantly, while the Ci increased significantly when the nitrogen application rate was low15.The experiments with Isatis indigotica demonstrate that the Pn, Gs, and Tr under the same irrigation level first increased and then decreased with the increase of the nitrogen application rate. The net photosynthetic rate, transpiration rate, and stomatal conductance of Isatis indigotica were improved by rational nitrogen application. Studies have reported similar findings in Isatis indigotica; with the decrease of N level, the net photosynthetic rate, transpiration rate, and stomatal conductance of leaves gradually decreased, while the intercellular CO2 concentration increased16,17. Under reasonable water and nitrogen coordination conditions, the synergistic effect of water and nitrogen increased, which effectively promoted the photosynthesis of Isatis indigotica. Under the condition of too much nitrogen or too little water, the antagonism of water and nitrogen was obvious, and the photosynthesis of Isatis indigotica was inhibited to a certain extent.Yield and water use efficiencyThe Isatis indigotica yield values presented are the average of two consecutive years of water–nitrogen trials (Fig. 2). The I. indigotica yields differed significantly between the water–nitrogen treatments; the W2N2 and W2N3 treatments had the highest yields at 7277.5 and 6820.5 kg hm−2, respectively. The lowest yield of 3264.5 kg hm−2 was recorded in the control treatment. The yields of all treatments were significantly higher than that of the control treatment. The yields of the W2N2 and W2N3 treatments were significantly higher than those of the W1N1 and the W3N1 treatments. With the increase of the nitrogen application rate, the yield first increased and then decreased under the same irrigation conditions.The water use efficiency values of Isatis indigotica presented are the average of 2 consecutive years of water–nitrogen trials (Fig. 2). The water use efficiency of Isatis indigotica differed significantly between the water–nitrogen treatments; the W1N2 and W2N2 treatments had the highest yields at 20.78 and 19.63 kg mm−1 hm−2, respectively. The lowest yield of 13.65 kg mm−1 hm−2 was recorded in the W3N1 treatment. The water use efficiency values of the W1N2 and W2N2 treatments were significantly higher than that of the W3N3 treatment, which was the treatment with excess water and nitrogen fertilizer. The water use efficiency decreased with the increase of irrigation under the same nitrogen application conditions. The water use efficiency first increased and then decreased with the increase in nitrogen application rate under the same irrigation conditions. The W2N2 treatment had the highest yield and water use efficiency. Therefore, the water–nitrogen coupling mode of medium water and medium nitrogen application achieved the highest yield and effectively saved water. This was mainly due to the moderate water and nitrogen to promote the photosynthesis of Isatis indigotica and lead to more dry matter accumulation, so as to increase the yield.Generally, appropriate water deficits can improve crop yield and water use efficiency18,19, and rational fertilization can increase crop yield, such as in fruit trees and vegetables20,21,22. The yield increase in the current experiment was probably related to reasonable water stress and reasonable nitrogen application; the W2N2 treatment had the highest yield and water use efficiency. However, excessive water and nitrogen reduced the yield and water use efficiency of Isatis indigotica. This was consistent with recent research reports23,24. Compared with the local flooding irrigation and excessive nitrogen fertilizer mode, the W2N2 treatment with moderate water and nitrogen application not only obtained a high yield but also significantly improved the water use efficiency. This method could reduce the effect of excessive water and fertilizer application on soil productivity and would be a better water and nitrogen management model for local Isatis indigotica production.QualityThe Isatis indigotica quality values presented are the average of two consecutive years of water–nitrogen trials (Fig. 3). These quality indicators mainly include the following content indicators: indigo, indirubin, (R, S)-goitrin, and polysaccharides. The Isatis indigotica quality indicators differed significantly between the water–nitrogen treatments. The CK treatment had the highest values of all quality indicators. Each quality indicator decreased gradually with the increase of water content under the same nitrogen application conditions. Each quality indicator decreased gradually with the increase of nitrogen application under the same water conditions. The (R, S)-goitrin content of the W2N2 treatment decreased by 6.5% compared with CK and by 3.9% compared with the W1N1 treatment.Water is the medium for improving crop quality. Generally, the crop quality was improved by a suitable water deficit25,26,27 and reasonable fertilization28,29,30. The quality of Isatis indigotica in the current experiment increased gradually with the decrease of water. The water deficit treatment increased the content of effective components and improved the quality of Isatis indigotica. The content of the effective components in all treatments reached the pharmacopoeia standard12. The quality indicator values of each treatment in the current experiment were significantly lower than those of the CK treatment, but there was little difference in the quality indicator values between each treatment. Moreover, the yield of the control treatment was much lower than that of other treatments. Therefore, the effective quality content of the control treatment was lower than other treatments. Excessive water and nitrogen inputs were not conducive to quality, which was not consistent with recent research reports31. Generally, the water-nitrogen coupling type of W2N60 was antagonism basing on the average yield of winter wheat in the 10 years32. Some scholars have studied the irrigation of jujube that WUE and ANUE of jujube cannot reach the maximum at the same time. Different ratio of water and nitrogen will produce coupling and antagonism33. The results showed that total N applications over 200 kg ha−1 did not increase yield or quality. Water deficit treatment could be increased the content of effective components and improve the quality of Isatis indigotica. Due to the high evaporation, moderate water stress and effective use of nitrogen fertilizer, the active components of Isatis indigotica were easier to accumulate in its roots. The synergistic effect of water and nitrogen will lead to the accumulation of active components in Isatis indigotica. More

1.Mimura, M. et al. Understanding and monitoring the consequences of human impacts on intraspecific variation. Evol. Appl. 10(2), 121–139 (2017).PubMed 
Article 
PubMed Central 

Google Scholar 
2.Leigh, D. M., Hendry, A. P., Vázquez-Domínguez, E. & Friesen, V. L. Estimated six per cent loss of genetic variation in wild populations since the industrial revolution. Evol. Appl. 12(8), 1505–1512 (2019).PubMed 
PubMed Central 
Article 

Google Scholar 
3.Habel, J. C., Husemann, M., Finger, A., Danley, P. D. & Zachos, F. E. The relevance of time series in molecular ecology and conservation biology. Biol. Rev. 89(2), 484–492 (2014).PubMed 
Article 
PubMed Central 

Google Scholar 
4.Klütsch, C. F. C. et al. Genetic changes caused by restocking and hydroelectric dams in demographically bottlenecked brown trout in a transnational subarctic riverine system. Ecol. Evol. 9(10), 6068–6081 (2019).PubMed 
PubMed Central 
Article 

Google Scholar 
5.Hansen, M. M., Fraser, D. J., Meier, K. & Mensberg, K.-L.D. Sixty years of anthropogenic pressure: A spatio-temporal genetic analysis of brown trout populations subject to stocking and population declines. Mol. Ecol. 18(12), 2549–2562 (2009).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
6.Savary, R. et al. Stocking activities for the Arctic char in Lake Geneva: Genetic effects in space and time. Ecol. Evol. 7(14), 5201–5211 (2017).PubMed 
PubMed Central 
Article 

Google Scholar 
7.Hughes, J. B., Daily, G. C. & Ehrlich, P. R. Population diversity: its extent and extinction. Science 278, 689–692 (1997).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
8.Perrier, C., Guyomard, R., Bagliniere, J.-L., Nikolic, N. & Evanno, G. Changes in the genetic structure of Atlantic salmon populations over four decades reveal substantial impacts of stocking and potential resiliency. Ecol. Evol. 3(7), 2334–2349 (2013).PubMed 
PubMed Central 
Article 

Google Scholar 
9.Vøllestad, L. A. & Hesthagen, T. Stocking of freshwater fish in Norway: management goals and effects. Nordic J. Freshwater Res. 75, 143–152 (2001).
Google Scholar 
10.Christie, M. R., Marine, M. L., French, R. A., Waples, R. S. & Blouin, M. S. Effective size of a wild salmonid population is greatly reduced by hatchery supplementation. Heredity 109, 254–260 (2012).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
11.Araki, H., Cooper, B. & Blouin, M. S. Carry-over effect of captive breeding reduces reproductive fitness of wild-born descendants in the wild. Biol. Lett. 5, 621–624 (2009).PubMed 
PubMed Central 
Article 

Google Scholar 
12.O’Sullivan, R. J. et al. Captive-bred Atlantic salmon released into the wild have fewer offspring than wild-bred fish and decrease population productivity. Proc. R. Soc. B 287, 20201671 (2020).PubMed 
PubMed Central 
Article 

Google Scholar 
13.Amundsen, P.-A. et al. Invasion of vendace Coregonus albula in a subarctic watercourse. Biol. Conserv. 88(3), 405–413 (1999).Article 

Google Scholar 
14.Jensen, H., Bøhn, T., Amundsen, P.-A. & Aspholm, P. E. Feeding ecology of piscivorous brown trout (Salmo trutta L.) in a subarctic watercourse. Ann. Zool. Fenn. 41(1), 319–328 (2004).
Google Scholar 
15.Jensen, H. et al. Predation by brown trout (Salmo trutta) along a diversifying prey community gradient. Can. J. Fish. Aquat. Sci. 65, 1831–1841 (2008).Article 

Google Scholar 
16.Jensen, H. et al. Food consumption rates of piscivorous brown trout (Salmo trutta) foraging on contrasting coregonid prey. Fish. Manag. Ecol. 22, 295–306 (2015).Article 

Google Scholar 
17.Haugland, Ø. Langtidsstudie av næringsøkologi og vekst hos storørret i Pasvikvassdraget. Mastergradsoppgave i biologi (Universitetet i Tromsø, Fakultet for Biovitenskap, fiskeri og økonomi, Institutt for arktisk og marin biologi, 2014).18.Gossieaux, P., Bernatchez, L., Sirois, P. & Garant, D. Impacts of stocking and its intensity on effective population size in Brook Charr (Salvelinus fontinalis) populations. Conserv. Genet. 20(4), 729–742 (2019).Article 

Google Scholar 
19.Pinter, K., Epifanio, J. & Unfer, G. Release of hatchery-reared brown trout (Salmo trutta) as a threat to wild populations? A case study from Austria. Fish. Res. 219, 105296 (2019).Article 

Google Scholar 
20.Wringe, B. F., Purchase, C. F. & Fleming, I. A. In search of a “cultured fish phenotype”: A systematic review, meta-analysis and vote-counting analysis. Rev. Fish Biol. Fish. 26(3), 351–373 (2016).Article 

Google Scholar 
21.Gossieaux, P. et al. Effects of genetic origin on phenotypic divergence in Brook Trout populations stocked with domestic fish. Ecosphere 11(5), e03119 (2020).Article 

Google Scholar 
22.Fleming, I. A., Jonsson, B. & Gross, M. R. Phenotypic divergence of sea-ranched, farmed, and wild salmon. Can. J. Fish. Aquat. Sci. 51, 2808–2824 (1994).Article 

Google Scholar 
23.Heath, D. D., Heath, J. W., Bryden, C. A., Johnson, R. M. & Fox, C. W. Rapid evolution of egg size in captive salmon. Science 299, 1738–1740 (2003).ADS 
CAS 
PubMed 
Article 

Google Scholar 
24.Naish, K. A., Seamons, T. R., Dauer, M. B., Hauser, L. & Quinn, T. P. Relationship between effective population size, inbreeding and adult fitness-related traits in a steelhead (Oncorhynchus mykiss) population released in the wild. Mol. Ecol. 22, 1295–1309 (2013).CAS 
PubMed 
Article 

Google Scholar 
25.Van Oosterhout, C., Weetman, D. & Hutchinson, W. F. Estimation and adjustment of microsatellite null alleles in nonequilibrium populations. Mol. Ecol. Notes 6(1), 255–256 (2006).Article 

Google Scholar 
26.Rousset, F. Genepop’007: a complete reimplementation of the Genepop software for Windows and Linux. Mol. Ecol. Resour. 8(6), 103–106 (2008).PubMed 
Article 

Google Scholar 
27.Peakall, R. & Smouse, P. E. GENALEX 6: Genetic analysis in Excel. Population genetic software for teaching and research. Mol. Ecol. Notes 6, 288–295 (2006).Article 

Google Scholar 
28.Peakall, R. & Smouse, P. E. GenAlEx 6.5: Genetic analysis in Excel. Population genetic software for teaching and research-an update. Bioinformatics 28(19), 2537–2539 (2012).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
29.Szpiech, Z. A., Jacobsson, M. & Rosenberg, N. A. ADZE: A rarefaction approach for counting alleles private to combinations of populations. Bioinformatics 24(21), 2498–2504 (2008).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
30.Waples, R. S. & Anderson, E. C. Purging putative siblings from population genetic data sets: A cautionary view. Mol. Ecol. 26(5), 1211–1224 (2017).PubMed 
Article 
PubMed Central 

Google Scholar 
31.Pritchard, J. K., Stephens, M. & Donnelly, P. Inference of population structure using multilocus genotype data. Genetics 155(2), 945–959 (2000).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
32.Jombart, T., Devillard, S. & Balloux, F. Discriminant analysis of principal components: A new method for the analysis of genetically structured populations. BMC Genet. 11, 94 (2010).PubMed 
PubMed Central 
Article 

Google Scholar 
33.Pew, J., Muir, P. H., Wang, J. & Frasier, T. R. Related: An R package for analysing pairwise relatedness from codominant molecular markers. Mol. Ecol. Resour. 15(3), 557–561 (2015).PubMed 
Article 
PubMed Central 

Google Scholar 
34.Piry, S., Luikart, G. & Cornuet, J.-M. Bottleneck: A computer program for detecting recent reductions in the effective population size using allele frequency data. J. Heredity 90(4), 502–503 (1999).Article 

Google Scholar 
35.Cornuet, J. M. & Luikart, G. Description and power analysis of two tests for detecting recent population bottlenecks from allele frequency data. Genetics 144(4), 2001–2014 (1996).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
36.Peery, M. Z. et al. Reliability of genetic bottleneck tests for detecting recent population declines. Mol. Ecol. 21(14), 3403–3418 (2012).PubMed 
Article 
PubMed Central 

Google Scholar 
37.Luikart, G. Usefulness of molecular markers for detecting population bottlenecks and monitoring genetic change. Ph. D. Thesis. (University of Montana, 1997).38.Do, C. et al. NEESTIMATOR v2: Re-implementation of software for the estimation of contemporary effective population size (Ne) from genetic data. Mol. Ecol. Resour. 14, 209–214 (2014).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
39.Waples, R. S. & Do, C. LDNE: A program for estimating effective population size from data on linkage disequilibrium. Mol. Ecol. Resour. 8, 753–756 (2008).PubMed 
Article 
PubMed Central 

Google Scholar 
40.Zhdanova, O. L. & Pudovkin, A. I. Nb_HetEx: A program to estimate the effective number of breeders. J. Hered. 99(6), 694–695 (2008).PubMed 
Article 

Google Scholar 
41.Nomura, T. Estimation of effective number of breeders from molecular coancestry of single cohort sample. Evol. Appl. 1, 462–474 (2008).PubMed 
PubMed Central 
Article 

Google Scholar 
42.Jones, O. R. & Wang, J. COLONY: A program for parentage and sibship inference from multilocus genotype data. Mol. Ecol. Resour. 10, 551–555 (2010).PubMed 
Article 
PubMed Central 

Google Scholar 
43.Wang, J. A. comparison of single-sample estimators of effective population sizes from genetic data. Mol. Ecol. 25, 4692–4711 (2016).PubMed 
Article 
PubMed Central 

Google Scholar 
44.Nei, M. & Chesser, R. K. Estimation of fixation indexes and gene diversities. Ann. Hum. Genet. 47(3), 253–259 (1983).CAS 
PubMed 
MATH 
Article 
PubMed Central 

Google Scholar 
45.Jost, L. Gst and its relatives do not measure differentiation. Mol. Ecol. 17(18), 4015–4026 (2008).PubMed 
Article 
PubMed Central 

Google Scholar 
46.Benjamini, Y. & Hochberg, Y. Controlling the false discovery rate: A practical and powerful approach to multiple testing. J. R. Stat. Soc. Ser. B (Methodological) 57(1), 289–300 (1995).MathSciNet 
MATH 

Google Scholar 
47.R Core Team. R: A Language and Environment for Statistical Computing. https://www.R-project.org/ (R Foundation for Statistical Computing, 2019).48.White, T., van der Ende, J. & Nichols, T. E. Beyond Bonferroni revisited: Concerns over inflated false positive research findings in the fields of conservation genetics, biology, and medicine. Conserv. Genet. 20, 927–937 (2019).Article 

Google Scholar 
49.Falush, D., Stephens, M. & Pritchard, J. K. Inference of population structure using multilocus genotype data: Linked loci and correlated allele frequencies. Genetics 164(4), 1567–1587 (2003).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
50.Hubisz, M. J., Falush, D., Stephens, M. & Pritchard, J. K. Inferring weak population structure with the assistance of sample group information. Mol. Ecol. Resour. 9(5), 1322–1332 (2009).PubMed 
PubMed Central 
Article 

Google Scholar 
51.Miller, M. A., Pfeiffer, W. & Schwartz, T. Creating the CIPRES science gateway for inference of large phylogenetic trees. in 2010 Gateway Computing Environments Workshop (GCE) 1–8 (2010).52.Besnier, F. & Glover, K. A. ParallelStructure: A R package to distribute parallel runs of the population genetics program STRUCTURE on multi-core computers. PLoS ONE 8(7), e70651 (2013).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
53.Li, Y.-L. & Liu, J.-X. StructureSelector: A web-based software to select and visualize the optimal number of clusters using multiple methods. Mol. Ecol. Resour. 18(1), 176–177 (2018).PubMed 
Article 

Google Scholar 
54.Puechmaille, S. J. The program structure does not reliably recover the correct population structure when sampling is uneven: Subsampling and new estimators alleviate the problem. Mol. Ecol. Resour. 16(3), 608–627 (2016).PubMed 
Article 

Google Scholar 
55.Kopelman, N. M., Mayzel, J., Jakobsson, M., Rosenberg, N. A. & Mayrose, I. CLUMPAK: A program for identifying clustering modes and packaging population structure inferences across K. Mol. Ecol. Resour. 15(5), 1179–1191 (2015).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
56.Anderson, E. C. & Dunham, K. K. The influence of family groups on inferences made with the program structure. Mol. Ecol. Resour. 8, 1219–1229 (2008).CAS 
PubMed 
Article 

Google Scholar 
57.Dray, S. & Dufour, A. The ade4 package: Implementing the duality diagram for ecologists. J. Stat. Softw. 22(4), 1–20 (2007).Article 

Google Scholar 
58.Levene, H. Robust tests for equality of variances. in Contributions to Probability and Statistics: Essays in Honor of Harold Hotelling (Olkin, I., Hotelling, H. et al. eds.). 278–292 (Stanford University Press, 1960).59.Kassambara, A. rstatix: Pipe-Friendly Framework for Basic Statistical Tests. R package version 0.4.0. https://CRAN.R-project.org/package=rstatix (2020).60.Wang, J. An estimator for pairwise relatedness using molecular markers. Genetics 160, 1203–1215 (2002).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
61.White, S. L., Miller, W. L., Dowell, S. A., Bartron, M. L. & Wagner, T. Limited hatchery introgression into wild brook trout (Salvelinus fontinalis) populations despite reoccurring stocking. Evol. Appl. 11(9), 1567–1581 (2018).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
62.Lehnert, S. J. et al. Multiple decades of stocking has resulted in limited hatchery introgression in wild brook trout (Salvelinus fontinalis) populations of Nova Scotia. Evol. Appl. 13(5), 1069–1089 (2020).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
63.Knudsen, C. M. et al. Comparison of life history traits between first-generation hatchery and wild upper Yakima River spring Chinook salmon. Trans. Am. Fish. Soc. 135, 1130–1144 (2006).Article 

Google Scholar 
64.Hansen, M. M. & Mensberg, K.-L.D. Admixture analysis of stocked brown trout populations using mapped microsatellite DNA markers: Indigenous trout persist in introgressed populations. Biol. Lett. 5, 656–659 (2009).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
65.Christie, M. R., Ford, M. J. & Blouin, M. S. On the reproductive success of early-generation hatchery fish in the wild. Evol. Appl. 7, 883–896 (2014).PubMed 
PubMed Central 
Article 

Google Scholar 
66.Fraser, D. J. et al. Population correlates of rapid captive-induced maladaptation in a wild fish. Evol. Appl. 12, 1305–1317 (2019).PubMed 
Article 

Google Scholar 
67.Fischer, J. R. et al. Growth, condition, and trophic relations of stocked trout in southern Appalachian mountain streams. Trans. Am. Fish. Soc. 148, 771–784 (2019).CAS 
Article 

Google Scholar 
68.Hendry, A. P. & Day, T. Population structure attributable to reproductive time: Isolation by time and adaptation by time. Mol. Ecol. 14, 901–916 (2005).CAS 
PubMed 
Article 

Google Scholar 
69.Gauthey, Z. et al. Brown trout spawning habitat selection and its effects on egg survival. Ecol. Freshwater Fish 26, 133–140 (2017).Article 

Google Scholar 
70.Dupont, P.-P., Bourret, V. & Bernatchez, L. Interplay between ecological, behavioural and historical factors in shaping the genetic structure of sympatric walleye populations (Sander vitreus). Mol. Ecol. 16, 937–951 (2007).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
71.Sandoval-Castillo, J. et al. SWINGER: A user-friendly computer program to establish captive breeding groups that minimize relatedness without pedigree information. Mol. Ecol. Resour. 17, 278–287 (2017).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar  More