Philippa Kaur

1.Dell, J. et al. Interaction diversity maintains resiliency in a frequently disturbed ecosystem. Front. Ecol. Evol. 7, 145 (2019).Article 

Google Scholar 
2.White, P. S. & Pickett, S. T. A. In The Ecology of Natural Disturbance and Patch Dynamics (eds S. T. A. Pickett & P. S. White) 3–13 (Academic Press, 1985).3.Newman, E. A. Disturbance ecology in the anthropocene. Front. Ecol. Evolut. https://doi.org/10.3389/fevo.2019.00147 (2019).Article 

Google Scholar 
4.Barnosky, A. D. et al. Approaching a state shift in Earth’s biosphere. Nature 486, 52–58 (2012).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
5.Yuan, Z., Jiao, F., Li, Y. & Kallenbach, R. L. Anthropogenic disturbances are key to maintaining the biodiversity of grasslands. Sci. Rep. 6, 22132 (2016).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
6.Hughes, A. R., Byrnes, J. E., Kimbro, D. L. & Stachowicz, J. J. Reciprocal relationships and potential feedbacks between biodiversity and disturbance. Ecol. Lett. 10, 849–864. https://doi.org/10.1111/j.1461-0248.2007.01075.x (2007).Article 
PubMed 
PubMed Central 

Google Scholar 
7.Connell, J. H. & Slatyer, R. O. Mechanisms of succession in natural communities and their role in community stability and organization. Am. Nat. 111, 1119–1144 (1977).Article 

Google Scholar 
8.Connell, J. H. Diversity in tropical rain forests and coral reefs. Science 199, 1302–1310 (1978).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
9.Fox, J. W. The intermediate disturbance hypothesis should be abandoned. Trends Ecol. Evol. 28, 86–92. https://doi.org/10.1016/j.tree.2012.08.014 (2013).Article 
PubMed 
PubMed Central 

Google Scholar 
10.Sheil, D. & Burslem, D. F. Disturbing hypotheses in tropical forests. Trends Ecol. Evol. 18, 18–26 (2003).Article 

Google Scholar 
11.Teixidó, N., Garrabou, J., Gutt, J. & Arntz, W. E. Recovery in Antarctic benthos after iceberg disturbance: Trends in benthic composition, abundance and growth forms. Mar. Ecol. Prog. Ser. 278, 1–16. https://doi.org/10.3354/meps278001 (2004).ADS 
Article 

Google Scholar 
12.Teixidó, N., Garrabou, J., Gutt, J. & Arntz, W. Iceberg disturbance and successional spatial patterns: the case of the shelf Antarctic benthic communities. Ecosystems 10, 143–158 (2007).Article 

Google Scholar 
13.Johst, K., Gutt, J., Wissel, C. & Grimm, V. Diversity and disturbances in the Antarctic megabenthos: Feasible versus theoretical disturbance ranges. Ecosystems 9, 1145–1155 (2006).Article 

Google Scholar 
14.Mackey, R. L. & Currie, D. J. The diversity-disturbance relationship: Is it generally strong and peaked?. Ecology 82, 3479–3492. https://doi.org/10.1890/0012-9658(2001) (2001).Article 

Google Scholar 
15.Huston, M. A. Disturbance, productivity, and species diversity: Empiricism vs. logic in ecological theory. Ecology 95, 2382–2396. https://doi.org/10.1890/13-1397.1 (2014).Article 

Google Scholar 
16.Smale, D. A., Brown, K. M., Barnes, D. K., Fraser, K. P. & Clarke, A. Ice scour disturbance in Antarctic waters. Science 321, 371. https://doi.org/10.1126/science.1158647 (2008).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
17.Griffiths, H. J., Danis, B. & Clarke, A. Quantifying Antarctic marine biodiversity: The SCAR-MarBIN data portal. Deep Sea Res. Part II 58, 18–29. https://doi.org/10.1016/j.dsr2.2010.10.008 (2011).ADS 
Article 

Google Scholar 
18.Grange, L. J. & Smith, C. R. Megafaunal communities in rapidly warming fjords along the West Antarctic Peninsula: Hotspots of abundance and beta diversity. PLoS ONE 8, e77917 (2013).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
19.Gutt, J., Griffiths, H. J. & Jones, C. D. Circumpolar overview and spatial heterogeneity of Antarctic macrobenthic communities. Mar. Biodivers. 43, 481–487. https://doi.org/10.1007/s12526-013-0152-9 (2013).Article 

Google Scholar 
20.Potthoff, M., Johst, K. & Gutt, J. How to survive as a pioneer species in the Antarctic benthos: Minimum dispersal distance as a function of lifetime and disturbance. Polar Biol. 29, 543–551 (2006).Article 

Google Scholar 
21.Convey, P. et al. The spatial structure of Antarctic biodiversity. Ecol. Monogr. 84, 203–244 (2014).Article 

Google Scholar 
22.Peck, L. S., Brockington, S., Vanhove, S. & Beghyn, M. Community recovery following catastrophic iceberg impacts in a soft-sediment shallow-water site at Signy Island, Antarctica. Mar. Ecol Progr. Ser. 186, 1–8 (1999).ADS 
Article 

Google Scholar 
23.Lee, H., Vanhove, S., Peck, L. & Vincx, M. Recolonisation of meiofauna after catastrophic iceberg scouring in shallow Antarctic sediments. Polar Biol. 24, 918–925. https://doi.org/10.1007/s003000100300 (2001).Article 

Google Scholar 
24.Armstrong, T. World Meteorological Organization. WMO sea-ice nomenclature. Terminology, codes and illustrated glossary. Edition 1970. Geneva, Secretariat of the World Meteorological Organization, 1970. [ix], 147 p. [including 175 photos]+ corrigenda slip. (WMO/OMM/BMO, No. 259, TP. 145.). J. Glaciol. 11, 148–149 (1972).25.Robinson, B. J., Barnes, D. K. & Morley, S. A. Disturbance, dispersal and marine assemblage structure: A case study from the nearshore Southern Ocean. Mar. Environ. Res. 160, 105025 (2020).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
26.Gutt, J., Starmans, A. & Dieckmann, G. Impact of iceberg scouring on polar benthic habitats. Mar. Ecol. Prog. Ser. 137, 311–316 (1996).ADS 
Article 

Google Scholar 
27.Barnes, D. K. A. & Conlan, K. E. Disturbance, colonization and development of Antarctic benthic communities. Philos. Trans. R. Soc. Lond. B Biol. Sci. 362, 11–38. https://doi.org/10.1098/rstb.2006.1951 (2007).Article 
PubMed 
PubMed Central 

Google Scholar 
28.Smale, D. A. Ecological traits of benthic assemblages in shallow Antarctic waters: Does ice scour disturbance select for small, mobile, secondary consumers with high dispersal potential?. Polar Biol. 31, 1225–1231. https://doi.org/10.1007/s00300-008-0461-9 (2008).Article 

Google Scholar 
29.Barnes, D. K. A. The influence of ice on polar nearshore benthos. J. Mar. Biol. Assoc. U.K. 79, 401–407 (1999).Article 

Google Scholar 
30.Gutt, J. On the direct impact of ice on marine benthic communities, a review. Polar Biol. 24, 553–564 (2001).Article 

Google Scholar 
31.Barnes, D. K. A. & Tarling, G. A. Polar oceans in a changing climate. Curr. Biol. 27, R454–R460. https://doi.org/10.1016/j.cub.2017.01.045 (2017).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
32.Barnes, D. K. A., Fleming, A., Sands, C. J., Quartino, M. L. & Deregibus, D. Icebergs, sea ice, blue carbon and Antarctic climate feedbacks. Philos. Trans. R. Soc. A Math. Phys. Eng. Sci. 376, 20170176. https://doi.org/10.1098/rsta.2017.0176 (2018).ADS 
Article 

Google Scholar 
33.Cook, A. J., Fox, A. J., Vaughan, D. G. & Ferrigno, J. G. Retreating glacier fronts on the Antarctic Peninsula over the past half-century. Science 308, 541–544. https://doi.org/10.1126/science.1104235 (2015).ADS 
CAS 
Article 

Google Scholar 
34.Cook, A. et al. Ocean forcing of glacier retreat in the western Antarctic Peninsula. Science 353, 283–286 (2016).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
35.Clarke, A. et al. Climate change and the marine ecosystem of the western Antarctic Peninsula. Philos. Trans. R. Soc. Lond. B Biol. Sci. 362, 149–166. https://doi.org/10.1098/rstb.2006.1958 (2007).Article 
PubMed 

Google Scholar 
36.Turner, J. & Comiso, J. Solve Antarctica’s sea-ice puzzle. Nat. News 547, 275 (2017).CAS 
Article 

Google Scholar 
37.Meredith, M. P. & King, J. C. Rapid climate change in the ocean west of the Antarctic Peninsula during the second half of the 20th century. Geophys. Res. Lett. https://doi.org/10.1029/2005GL024042 (2005).Article 

Google Scholar 
38.Barnes, D. K. A. & Souster, T. Reduced survival of Antarctic benthos linked to climate-induced iceberg scouring. Nat. Clim. Chang. 1, 365–368. https://doi.org/10.1038/nclimate1232 (2011).ADS 
Article 

Google Scholar 
39.Parkinson, C. L. Global sea ice coverage from satellite data: Annual cycle and 35-yr trends. J. Clim. 27, 9377–9382. https://doi.org/10.1175/jcli-d-14-00605.1 (2014).ADS 
Article 

Google Scholar 
40.Rogers, A. et al. Antarctic futures: An assessment of climate-driven changes in ecosystem structure, function, and service provisioning in the Southern Ocean. Ann. Rev. Mar. Sci. 12, 87–120 (2020).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
41.Morley, S. A. et al. Global drivers on Southern Ocean ecosystems: Changing physical environments and anthropogenic pressures in an Earth system. Front. Mar. Sci. 7, 1097 (2020).Article 

Google Scholar 
42.Barnes, D. K. et al. Blue carbon gains from glacial retreat along Antarctic fjords: What should we expect?. Glob. Change Biol. 26, 2750–2755 (2020).ADS 
Article 

Google Scholar 
43.Barnes D. K. A. Blue carbon on polar and subpolar seabeds. In Carbon capture, utilization and sequestration (InTech, 2018). https://doi.org/10.5772/intechopen.78237.44.Bowler, D. et al. The geography of the Anthropocene differs between the land and the sea. bioRxiv https://doi.org/10.1101/432880 (2019).Article 

Google Scholar 
45.Arntz, W., Brey, T. & Gallardo, V. Antarctic zoobenthos. Oceanogr. Mar. Biol. 32, 241–304 (1994).
Google Scholar 
46.Clarke, A. Marine benthic populations in Antarctica: Patterns and processes. Antarct. Res. Ser. 70, 373–388 (1996).Article 

Google Scholar 
47.Fillinger, L., Janussen, D., Lundälv, T. & Richter, C. Rapid glass sponge expansion after climate-induced Antarctic ice shelf collapse. Curr. Biol. 23, 1330–1334 (2013).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
48.Clarke, A., Meredith, M. P., Wallace, M. I., Brandon, M. A. & Thomas, D. N. Seasonal and interannual variability in temperature, chlorophyll and macronutrients in northern Marguerite Bay, Antarctica. Deep Sea Res. Part II 55, 1988–2006. https://doi.org/10.1016/j.dsr2.2008.04.035 (2008).ADS 
Article 

Google Scholar 
49.Barnes, D. K. A. Iceberg killing fields limit huge potential for benthic blue carbon in Antarctic shallows. Glob. Chang. Biol. 23, 2649–2659. https://doi.org/10.1111/gcb.13523 (2017).ADS 
Article 
PubMed 
PubMed Central 

Google Scholar 
50.Pinkerton, M., Bradford-Grieve, J., Bowden, D. & Cummings, V. Benthos: Trophic modelling of the Ross Sea. Support. Docum. CCAMLR Sci. 17, 1–31 (2010).
Google Scholar 
51.Pielou, E. Shannon’s formula as a measurement of species diversity: It’s use and disuse. Am. Nat. 100, 463–465 (1966).Article 

Google Scholar 
52.Fisher, R. A., Corbet, A. S. & Williams, C. B. The relation between the number of species and the number of individuals in a random sample of an animal population. J. Anim. Ecol. 1, 42–58 (1943).Article 

Google Scholar 
53.Everitt, B. & Skrondal, A. The Cambridge Dictionary of Statistics Vol. 106 (Cambridge University Press, Cambridge, 2002).MATH 

Google Scholar 
54.Smale, D. A., Barnes, D. K. A. & Fraser, K. P. P. The influence of ice scour on benthic communities at three contrasting sites at Adelaide Island, Antarctica. Aust. Ecol. 32, 878–888. https://doi.org/10.1111/j.1442-9993.2007.01776.x (2007).Article 

Google Scholar 
55.Peck, L. S., Convey, P. & Barnes, D. K. A. Environmental constraints on life histories in Antarctic ecosystems: Tempos, timings and predictability. Biol. Rev. 81, 75–109. https://doi.org/10.1017/s1464793105006871 (2006).Article 
PubMed 
PubMed Central 

Google Scholar 
56.Waller, C., Worland, M., Convey, P. & Barnes, D. Ecophysiological strategies of Antarctic intertidal invertebrates faced with freezing stress. Polar Biol. 29, 1077–1083 (2006).Article 

Google Scholar 
57.Barnes, D. K. A. Polar zoobenthos blue carbon storage increases with sea ice losses, because across-shelf growth gains from longer algal blooms outweigh ice scour mortality in the shallows. Glob. Chang Biol. 23, 5083–5091. https://doi.org/10.1111/gcb.13772 (2017).ADS 
Article 
PubMed 
PubMed Central 

Google Scholar 
58.Smith, C. R., Mincks, S. & DeMaster, D. J. A synthesis of bentho-pelagic coupling on the Antarctic shelf: Food banks, ecosystem inertia and global climate change. Deep Sea Res. Part II 53, 875–894 (2006).ADS 
Article 

Google Scholar 
59.Jansen, J. et al. Abundance and richness of key Antarctic seafloor fauna correlates with modelled food availability. Nat. Ecol. Evolut. 2, 71–80 (2018).Article 

Google Scholar 
60.Henley, S. F. et al. Changing biogeochemistry of the Southern Ocean and its ecosystem implications. Front. Mar. Sci. 7, 581 (2020).Article 

Google Scholar 
61.Marshall, G. J. et al. Causes of exceptional atmospheric circulation changes in the Southern Hemisphere. Geophys. Res. Lett. 31, 14 (2004).Article 

Google Scholar 
62.Ashton, G. V., Morley, S. A., Barnes, D. K., Clark, M. S. & Peck, L. S. Warming by 1 C drives species and assemblage level responses in Antarctica’s marine shallows. Curr. Biol. 27, 2698-2705e2693 (2017).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
63.Riesgo, A. et al. Some like it fat: Comparative ultrastructure of the embryo in two demosponges of the genus Mycale (order poecilosclerida) from Antarctica and the Caribbean. PLoS ONE 10, e0118805 (2015).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
64.Toszogyova, A. & Storch, D. Global diversity patterns are modulated by temporal fluctuations in primary productivity. Glob. Ecol. Biogeogr. 28, 1827–1838 (2019).Article 

Google Scholar 
65.Clark, G. F. et al. Light-driven tipping points in polar ecosystems. Glob. Change Biol. 19, 3749–3761 (2013).ADS 
Article 

Google Scholar 
66.Brockington, S., Clarke, A. & Chapman, A. Seasonality of feeding and nutritional status during the austral winter in the Antarctic sea urchin Sterechinus neumayeri. Mar. Biol. 139, 127–138 (2001).Article 

Google Scholar 
67.Fratt, D. B. & Dearborn, J. Feeding biology of the Antarctic brittle star Ophionotus victoriae (Echinodermata: Ophiuroidea). Polar Biol. 3, 127–139 (1984).Article 

Google Scholar 
68.Sahade, R., Tatián, M. & Esnal, G. B. Reproductive ecology of the ascidian Cnemidocarpa verrucosa at Potter Cove, South Shetland Islands, Antarctica. Mar. Ecol. Progr. Ser. 272, 131–140 (2004).ADS 
Article 

Google Scholar 
69.Dayton, P. K. et al. Recruitment, growth and mortality of an Antarctic hexactinellid sponge, Anoxycalyx joubini. PLoS ONE 8, e56939 (2013).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
70.Vacchi, M., Cattaneo-Vietti, R., Chiantore, M. & Dalù, M. Predator-prey relationship between the nototheniid fish Trematomus bernacchii and the Antarctic scallop Adamussium colbecki at Terra Nova Bay (Ross Sea). Antarct. Sci. 12, 64–68 (2000).ADS 
Article 

Google Scholar 
71.Sheil, D. & Burslem, D. F. Defining and defending Connell’s intermediate disturbance hypothesis: a response to Fox. Trends Ecol. Evol. 28, 571–572. https://doi.org/10.1016/j.tree.2013.07.006 (2013).Article 
PubMed 
PubMed Central 

Google Scholar  More

A. Strain sampling and isolationBradyrhizobium is a commonly occurring genus in soil [21]. Closely related Bradyrhizobium diazoefficiens (previously Bradyrhizobium japonicum) strains were isolated from soil, as previously described [20, 22]. In brief, Bradyrhizobium isolates that formed symbiotic associations with a foundational legume species in the censused region, Acacia acuminata, were isolated from soil sampled along a large region spanning ~300,000 km2 in South West Australia, a globally significant biodiversity hotspot [23]. In total 60 soil samples were collected from twenty sites (3 soil samples per site; Supplementary Fig. S1) and 380 isolates were sequenced (19 isolates per site, 5 or 6 isolates per soil sample, each isolate re-plated from a single colony at least 2 times). Host A. acuminata legume plants were inoculated with field soil in controlled chamber conditions and isolates were cultured on Mannitol Yeast agar plates from root nodules (see [20, 22] for details). A total of 374 strains were included in this study after removing 5 contaminated samples and one sample that was a different Bradyrhizobium species; non- Bradyrhizobium diazoefficiens sample removal was determined from 16S rRNA sequences extracted from draft genome assemblies (Method C) using RNAmmer [24].B. Environmental variation among sampled sitesIn this study, I focus on environmental factors (temperature, rainfall, soil pH and salinity) previously identified to impact either rhizobia growth performance, functional fitness or persistence in soil [25,26,27,28] and where a directionality of rhizobial stress response could be attributed with respect to environmental variation present in the sampled region (i.e. stress occurs at high temperatures, low rainfall, high acidity and high salinity). Each environmental factor was standardised to a mean of 0 and a standard deviation of 1, and pH and rainfall scales were reversed to standardise stress responses directions so that low stress is at low values and high stress is at high values for all factors (Supplementary Fig. S2). Additionally, salinity was transformed using a log transformation (log(x + 0.01) to account for some zeroes) prior to standardisation.C. Isolate sequencing and pangenome annotationIllumina short reads (150 bp paired-end) were obtained and draft genome assemblies were generated using Unicycler from a previous study [29]. Resulting assemblies were of good assembly quality (99.2% of all strains had >95.0% genome completeness score according to BUSCO [30]; Table S1; assembled using reads that contained nominal 0.016 ± 0.00524% non-prokaryotic DNA content across all 374 isolates, according to Kraken classification [31]). Protein coding regions (CDS regions) were identified using Prokka [32] and assembled into a draft pangenome using ROARY [33], which produced a matrix of counts of orthologous gene clusters (i.e. here cluster refers to a set of protein-coding sequences containing all orthologous variants from all the different strains, grouped together and designated as a single putative gene). Gene clusters that occurred in 99% of strains were designated as ‘core genes’ and used to calculate the ‘efficiency of selection’ [34, 35] (measured as dN/dS, Method G.2) and population divergence measured as Fixation Index ‘Fst’, Method H) across each environmental stress factor. The identified gene clusters were then annotated using eggNOG-mapper V2 [36] and the strain by gene cluster matrix was reaggregated using the Seed ortholog ID returned by eggNOG-mapper as the protein identity. Out of the total 2,744,533 CDS regions identified in the full sample of 374 strains, eggNOG-mapper was able to assign 2,612,345 of them to 91,230 unique Seed orthologs. These 91,230 protein coding genes constituted the final protein dataset for subsequent analyses.D. Calculation and statistical analysis of gene richness and pangenome diversity along the stress gradientGene richness was calculated as the total number of unique seed orthologues for each strain (i.e. genome). Any singleton genes that occurred in only a single strain, as well as ‘core’ genes that occurred in every strain (for symmetry, and because these are equally uninformative with respect to variation between strains) were removed, leaving 74,089 genes in this analysis. Gene richness (being count data) was modelled on a negative binomial distribution (glmer.nb function) as a function of each of the four environmental stressors as predictors using the lme4 package in R [37], also accounting for hierarchical structure in the data by including site and soil sample as random effects.To rule out potentially spurious effects of assembly quality (i.e. missed gene annotations due to incomplete draft genomes) on key findings, I confirmed no significant association between gene richness and genome completeness (r = 0.042, p = 0.4224, Fig. S3).Finally, pangenome diversity was calculated as the total number of unique genes that occurred in each soil sample (since multiple strains were isolated from a single soil sample). Pangenome diversity was modelled the same as gene richness, except here soil sample was not included as a random effect.E. Calculation of network and duplication traits for each geneI used the seed orthologue identifier from eggNOG-mapper annotations to query matching genes within StringDB ([38]; https://string-db.org/), which collects information on protein-protein interactions. Out of 91,230 query seed orthologues, 73,126 (~80%) returned a match in STRING. For matching seed orthologue hits, a network was created by connecting any proteins that were annotated as having pairwise interactions in the STRING database using the igraph package in R [39]. Two vertex-based network metrics were calculated for each gene: betweenness centrality, which measures a genes tendency to connect other genes in the gene network, and mean cosine similarity, which is a measure of how much a gene’s links to other genes are similar to other genes.Betweenness centrality was calculated using igraph (functional betweenness). For mean cosine similarity, a pairwise cosine similarity was first calculated between all genes. To do this, the igraph network object was converted into a (naturally sparse yet large) adjacency matrix and the cosSparse function in qlcMatrix in R [40] was used to calculate cosine similarity between all pairs of genes. To obtain an overall cosine similarity trait value for each gene, the average pairwise cosine similarity to all other genes in the network was calculated.Finally, gene duplication level was calculated for each gene as one additional measure of ‘redundancy’, by calculating the average number of gene duplicates found within the same strain. Duplicates were identified as CDS regions with the same Seed orthologue ID.F. Gene distribution modelsTo determine how gene traits predict accessory genome distributions patterns along the stress gradients, I first calculated a model-based metric (hereafter and more specifically a standardised coefficient, ‘z-score’) of the relative tendency of each gene to be found in different soil samples across the four stress gradients (heat, salinity, acidity, and aridity). This was achieved by modelling each gene’s presence or absence in a strain as a function of the four stress gradients, with site and soil sample as a random effect, using a binomial model in lme4 (the structure of the model being the same as the gene richness model, only the response is different). To reduce computational overhead, these models were only run for the set of genes that were used in the gene richness analysis (e.g. after removing singletons and core genes), and which had matching network traits calculated (e.g. they occurred in the STRING database; n = 64,867 genes). Distribution models were run in tandem for each gene using the manyany function in the R package mvabund [41]. Standardised coefficients, or z-scores (coefficient/standard error) indicating the change in the probability of occurrence for each gene across each of the stress gradients were extracted. More negative coefficients correspond to genes that are more likely to be absent in high stress (and vice versa for positive coefficients).To determine how network and duplication traits influence the distribution of genes across the stress gradient, I performed a subsequent linear regression model where the gene’s z-scores was the response and gene traits as predictors. The environmental stress type (i.e. acidity, aridity, heat and salinity) was included as a categorical predictor, and the interaction between stress category and gene function traits were used to infer the influence of gene function traits on gene distributions in a given stress type (see Supplementary Methods S1 for z-score transformation).G. Quantifying molecular signals of natural selection on accessory and core genesTo examine molecular signatures of selection in accessory and core genes, I calculated dN/dS for a subsample of the total pool (n=74,089 genes), which estimates the efficiency of selection [34, 35]. Two major questions relevant to dN/dS that are addressed here require a different gene subsampling approach:(1) Do variable environmental stress responses lead to different dN/dS patterns among accessory genes?Here, I subsampled accessory genes (total accessory gene pool across 374 strains, 74,089) to generate and compare dN/dS among 3 categorical groups, each representing a different level of stress response based on their z-scores (accessory genes that either strongly increase, decrease or have no change in occurrence as stress increases; n = 1000 genes/category; see Supplementary Methods S2 for subsample stratification details).For each gene, sequences were aligned using codon-aware alignment tool, MACSE v2 [42]. dN/dS was estimated by codon within each gene using Genomegamap’s Bayesian model-based approach [43], which is a phylogeny-free method optimised for within bacterial species dN/dS calculation (see Supplementary Methods S3 for dN/dS calculation/summarisation; S9 for xml configuration). The proportion of codons with dN/dS that were credibly less than 1 (purifying selection) and those credibly greater than 1 (positive selection) were analysed, with respect to the genes’ occurrence response to stress. Specifically, I modelled the proportion of codons with dN/dS  1 was overall too low to analyse in this way, so the binary outcome (a gene with any codons with dN/dS  > 1 or not) was modelled using a binomial response model with the response categories as predictors (see Supplementary Methods S4 for details of both models).(2) Does dN/dS among microbial populations vary across environmental stress?Here, I compared the average change in dN/dS in core genes present across all environments at the population level (i.e. all isolates from one soil sample), which is used here to measure the change in the efficiency of selection across each stress gradient. Core genes were used since they occur in all soil samples, allowing a consistent set and sample size of genes to be used in the population-level dN/dS calculation. This contrasts with the previous section, which quantifies gene-level dN/dS on extant accessory genes that intrinsically have variable presence or absence across soil samples. For computational feasibility, 500 core genes were subsampled (total core 1015 genes) and, for each gene, individual strain variants were collated into a single fasta file based on soil sample membership, such that dN/dS could be calculated separately for each gene within each soil sample (i.e. 60 soil samples × 500 genes = 30,000 fasta files). Each fasta file was then aligned in MACSE and dN/dS calculated using the same methodology for accessory genes (Supplementary Method S3). This enabled the average dN/dS in a sample to be associated with soil-sample level environmental stress variables. Specifically, I modelled the mean proportion of codons with dN/dS  1 due to overall rarity of positive selection (average proportion of genes where at least 1 codon with dN/dS  > 1 was ~0.006). This low level of positive selection is expected for core genes which tend to be under strong selective constraint.H. Calculation and analysis of Fixation index (Fst) along stress gradientsUsing the core genome alignment (all SNPs among 1015 core genes) generated previously with ROARY, I computed pairwise environmentally-stratified Fst. Each soil sample (n = 60) was first placed into one of 5 bins based on their distances in total environmental stress space (using all four stress gradients), with the overall goal of generating roughly evenly sized bins such that the environmental similarity of stress was greater within bins than between (see Supplementary Methods S6 and Fig. S4 for clustering algorithm details). Next, fasta alignments were converted to binary SNPs using the adegenet package. Pairwise Fst between all strains (originating from a particular soil sample) within a single bin was calculated using StAMPP in R [44]. For each strain pair, the average of the two stress gradient values was assigned.The relationship between pairwise Fst and the average stress value was evaluated using a linear regression model with each of the four stress values as predictors. Since the analysis uses pairwise data (thus violating standard independence assumptions), the significance of the relationship was determined using a permutation test (see Supplementary Methods S7 for details).I. Chromosomal structure analysis of gene loss patternsTo gain insight into structural variation and test for regional hotspots in gene loss along the chromosome, I mapped each gene’s stress response (i.e. probability of loss or gain indicated by each genes z-score) onto a completed Bradyrhizobium genome (strain ‘36_1’ from the same set of 374 strains (Genbank CP067102.1; [45]). Putative CDS positions on the complete genome were determined using Prokka and annotated with SEED orthologue ID’s using eggNOG-mapper. Matching accessory genes derived from the full set of 374 incomplete draft genomes (n = 74,089) were mapped to positions on the complete genome (6274 matches). The magnitude of gene loss or gain (as measured by their standardised z-scores for each environment from the gene distribution models; see Method F) was then modelled across the genome using a one-dimensional spatial smoothing model. This model was implemented in R INLA [46] (www.r-inla.org), and models a response in a one-dimensional space using a Matern covariance-based random effect. The method uses an integrated nested Laplace approximation to a Bayesian posterior distribution, with a cyclical coordinate system to accommodate circular genomes. The model accounts for spatial non-independence among sites and produces a continuous posterior distribution of average z-score predictions along the entire genome, which was then used to visualise potential ‘hotspots’ of gene loss or gain. The modelling procedure was repeated, instead with gene network traits, such that model predictions of similarity and betweenness could be visualised on the reference chromosome. More

State Key Laboratory of Tree Genetics and Breeding, Chinese Academy of Forestry, Beijing, ChinaZhilin Yuan, Huanshen Wei & Long PengResearch Institute of Subtropical Forestry, Chinese Academy of Forestry, Hangzhou, ChinaZhilin Yuan, Xinyu Wang, Huanshen Wei & Long PengFungal Genomics Laboratory (FungiG), College of Resources and Environmental Sciences, Nanjing Agricultural University, Nanjing, ChinaIrina S. Druzhinina & Feng CaiDepartment of Food Science, University of Massachusetts, Amherst, MA, USAJohn G. GibbonsState Key Laboratory of Ecological Pest Control for Fujian and Taiwan Crops, College of Plant Protection, Fujian Agriculture and Forestry University, Fuzhou, ChinaZhenhui ZhongDepartment of Molecular, Cell and Developmental Biology, University of California, Los Angeles, CA, USAZhenhui ZhongDepartment of Plant Biotechnology and Bioinformatics, Ghent University, Ghent, BelgiumYves Van de PeerVIB Center for Plant Systems Biology, Ghent, BelgiumYves Van de PeerCentre for Microbial Ecology and Genomics, Department of Biochemistry, Genetics and Microbiology, University of Pretoria, Hatfield, South AfricaYves Van de PeerAdaptive Symbiotic Technologies, University of Washington, Seattle, WA, USARussell J. RodriguezKey Laboratory of National Forestry and Grassland Administration for Orchid Conservation and Utilization at College of Landscape Architecture, Fujian Agriculture and Forestry University, Fuzhou, ChinaZhongjian LiuState Key Laboratory of Mycology, Institute of Microbiology, Chinese Academy of Sciences, Beijing, ChinaQi Wu & Guohui ShiKey Laboratory of Plant Resources Conservation and Sustainable Utilization, South China Botanical Garden, Chinese Academy of Sciences, Guangzhou, ChinaJieyu WangBeijing Advanced Innovation Center for Tree Breeding by Molecular Design, Beijing Forestry University, Beijing, ChinaFrancis M. MartinUniversité de Lorraine, INRAE, UMR Interactions Arbres/Micro-Organismes, Centre INRAE Grand Est Nancy, Champenoux, FranceFrancis M. Martin More

Detecting alternative attractors in ecosystem dynamicsWe use empirical dynamical modeling, a set of equation-free tools for analyzing non-linear time series (for a review and assumptions see25,26, respectively), to test if the temporal dynamics of alternative dynamical regimes are qualitatively different. Empirical dynamic modeling builds fundamentally on Takens embedding theorem, which shows that attractors of multi-dimensional dynamical systems can be reconstructed using higher order lags of its embedded time series27. However, if a dynamical system has gone through a bifurcation, or switched to an alternative basin of attraction, attractors are qualitative dissimilar in the two regimes. Theoretically, this infers that it should be possible to reconstruct the attractor of one regime using information from the same regime, but not from the other regime. In practice, this implies that if a model (attractor reconstruction) based on one dynamical regime is used to predict the dynamics of variables from the same dynamical regime predictions should be accurate (i.e. low prediction errors), whereas if an attractor reconstruction based on one dynamical regime is used to predict the dynamics of variables of another attractor predictions should be less accurate (i.e. high prediction errors). We make use of this idea by specifically testing if prediction errors of across and within regime predictions are different. As explained below this idea can be used for both univariate and multivariate time series data.Univariate approachUnivariate attractor reconstructions can be found using the simplex algorithm28,29. First, for a given dynamical regime, a time series can be split into a library of vectors, and each vector is described by$${underline{y}}_{A}(t)= < {Y}_{A}(t),{Y}_{A}(t-1),{Y}_{A}(t-2),ldots ,{Y}_{A}(t-(E-1)) > ,$$
(1)
where ({Y}_{A}(t)) is an observation of variable Y at time t in dynamical regime A and E is the reconstructed attractors embedding dimension. Using the simplex projection algorithm, a one-step ahead forecast is produced as follows:$${hat{Y}}_{A}(t+1)|{{{{{{boldsymbol{M}}}}}}}_{B}=mathop{sum}limits_{m=1ldots E+1}{w}_{m}{Y}_{B}({t}_{m}+1),$$
(2)
where tm is a time index of an observation in dynamical regime B, E is the embedding dimension of regime B, and wm is an exponential weighting described by:$${w}_{m}={u}_{m}/mathop{sum}limits_{n=1,ldots ,E+1}{u}_{n},$$
(3)
where n and m belongs to the set of the E+1 nearest neighbors of vector ({underline{y}}_{A}(t)) in the set of vectors ({{underline{y}}_{B}({t}_{m})}), ({u}_{m}=exp {-d[{underline{y}}_{A}(t),{underline{y}}_{B}({t}_{m})]/d[{underline{y}}_{A}(t),{underline{y}}_{B}({t}_{1})]}), and (d[{underline{y}}_{A}(t),{underline{y}}_{B}({t}_{1})],)is the Euclidean distance between the prediction vector ({underline{y}}_{A}(t)) and its nearest neighbor ({underline{y}}_{B}({t}_{1})) in the set ({{underline{y}}_{B}({t}_{m})}).The only parameter that is estimated using the simplex algorithm is the embedding dimension E. This parameter is estimated by optimizing the correlation between observations (({Y}_{A}(t+1))) and predictions (({hat{Y}}_{A}(t+1)|{{{{{{boldsymbol{M}}}}}}}_{A})) using a leave-one-out cross validation approach (See Supplementary Discussion). The embedding dimension E and its corresponding set of E-dimensional vectors (Eq. 1) constitutes the reconstructed attractor, MA, of a given dynamical regime A. This reconstructed attractor (MA) is then used to predict data for both the same dynamical regime (({hat{Y}}_{A}(t+1)|{{{{{{boldsymbol{M}}}}}}}_{A})), and the contrasting dynamical regime ({hat{Y}}_{B}(t+1)|{{{{{{boldsymbol{M}}}}}}}_{A}). Likewise, the reconstructed attractor MB can be used to predict time series dynamics from both dynamical regimes; that is, ({hat{Y}}_{A}(t+1)|{{{{{{boldsymbol{M}}}}}}}_{B}) and ({hat{Y}}_{B}(t+1)|{{{{{{boldsymbol{M}}}}}}}_{B}), respectively.Multivariate approachA multivariate time series describes a number of simultaneously evolving variables. For example, a bivariate time series can be described by variables X and Y. For such time series, Sugihara et al.30 developed an approach for testing if two variables (time series) are dynamically coupled. Their methodology builds on the fact that a reconstructed attractor should map 1:1 to the original attractor on which the reconstruction is based. This infers that two attractor reconstructions (based on two different variables) should also map 1:1 to each other30. Practically, this means that if two variables are dynamically coupled one-time series should be predictable based on an attractor reconstruction of another variable. However, if a dynamical system has gone through a bifurcation, or potentially switched to an alternative basin of attraction, a new set of rules will govern the dynamics of the system. Hence, a new attractor should have emerged. Now, since this new attractor is most likely governed by a new set of rules it should be difficult to predict the dynamics of this new alternative attractor based on information from the former attractor. Thus, if one variable in one dynamical regime is used to predict another variable in another dynamical regime, predictions should be biased. Yet, if one variable from one dynamical regime is used to predict another variable from the same regime predictions should be more accurate.The simplex algorithm can be used to make predictions of a variable Y using a time series of another variable X30. Predictions are produced as follows:$${hat{Y}}_{{{{{{boldsymbol{A}}}}}}}(t)|{{{{{{boldsymbol{M}}}}}}}_{B}=mathop{sum}limits_{m=1ldots E+1}{w}_{m}{Y}_{B}({t}_{m}),$$
(4)
where tm is the time series index of a vector of variable X of dynamical regime B, wm is an exponential weighting based on variable X:$${w}_{m}={u}_{m}/mathop{sum}limits_{n=1,ldots ,E+1}{u}_{n},$$
(5)
where n and m belongs to the set of the E+1 nearest neighbors of ({underline{x}}_{A}(t)) in ({{underline{x}}_{B}({t}_{m})}), ({u}_{m}=exp {-d[{underline{x}}_{A}(t),{underline{x}}_{B}({t}_{m})]/d[{underline{x}}_{A}(t),{underline{x}}_{B}({t}_{1})]}), and (d[{underline{x}}_{A}(t),{underline{x}}_{B}({t}_{1})],)is the Euclidean distance between the prediction vector(,{underline{x}}_{A}(t)) and its nearest neighbor ({underline{x}}_{B}({t}_{1})) in dynamical regime (B).The reconstructed attractors, MA and MB, for each variable and regime are found using the univariate simplex algorithm described above28,29,30. Similar to the univariate case, the reconstructed attractor (MA) is used to predict data from the same dynamical regime (({hat{Y}}_{{{{{{boldsymbol{A}}}}}}}(t)|{{{{{{boldsymbol{M}}}}}}}_{A})), and to predict time series of a contrasting dynamical regime (({hat{Y}}_{{{{{{boldsymbol{A}}}}}}}(t)|{{{{{{boldsymbol{M}}}}}}}_{B})). Yet, it is important to stress that MA here reflects an attractor reconstruction based on a variable that is not being predicted (that is, variable X is used to predict variable Y). This prediction approach thus infers that predictions are made on data that was not used to fit the model (X predicts Y and vice versa). Thus, neither across nor within regime predictions are made on data used to fit a model.Test statisticWe used mean absolute prediction errors to test for difference between across and within regime predictions. Alternative metrics, such as mean sum of square errors, can also be used. However, since our approach gives skewed prediction errors we used mean absolute prediction errors to reduce the impact of extreme values. Further, since the absolute prediction errors are non-normally distributed we used a permutation test. The null hypothesis that is tested reads:$$H0:{{{{{rm{MAP{E}}}}}}}_{A} < {{{{{rm{MAP{E}}}}}}}_{w},$$ (6) where MAPEA is the mean absolute prediction error for across regime predictions (that is, ({{{{{rm{MAP{E}}}}}}}_{A}=frac{1}{n}mathop{sum}limits_{t=1:n}{{{{{rm{abs}}}}}}({hat{Y}}_{{{{{{{boldsymbol{M}}}}}}}_{{{{{{boldsymbol{A}}}}}}}}(t)|{{{{{{boldsymbol{M}}}}}}}_{B}-{Y}_{{{{{{{boldsymbol{M}}}}}}}_{{{{{{boldsymbol{A}}}}}}}}(t))), and ({{{{{rm{MAP{E}}}}}}}_{w}) is the mean absolute prediction error for within regime predictions (that is, ({{{{{rm{MAP{E}}}}}}}_{w}=frac{1}{n}mathop{sum}limits_{t=1:n}{{{{{rm{abs}}}}}}({hat{Y}}_{{{{{{{boldsymbol{M}}}}}}}_{{{{{{boldsymbol{A}}}}}}}}(t)|{{{{{{boldsymbol{M}}}}}}}_{A}-{Y}_{{{{{{{boldsymbol{M}}}}}}}_{{{{{{boldsymbol{A}}}}}}}}(t))). A test is consider significant if observed difference in across and within regime mean prediction errors is larger than the 95th percentile of 1000 permuted data sets.Food-chain modelWe used a food-chain model parameterized as in McCann and Yodzis31 to simulate food-chain dynamics:$$frac{{{{{{rm{d}}}}}}R}{{{{{{rm{d}}}}}}t}=Rleft(1-frac{R}{K}right)-frac{{x}_{c}{y}_{c}CR}{R+{R}_{0}}$$ (7) $$frac{{{{{{rm{d}}}}}}C}{{{{{{rm{d}}}}}}t}={x}_{c}Cleft(-1+frac{{y}_{C}R}{R+{R}_{0}}right)-frac{{x}_{P}{y}_{P}PC}{C+{C}_{0}}$$$$,frac{{{{{{rm{d}}}}}}P}{{{{{{rm{d}}}}}}t}={x}_{P}Pleft(-1+frac{{y}_{P}C}{C+{C}_{0}}right),$$where R is the resource density, C consumer density, and P predator density. All parameters, except half-saturation constants R0 (here set to 0.16129) and C0 (here set to 0.5), and resource carrying capacity K, are derived from bioenergetics and body size allometry30 (xc = 0.4, yc = 2.009, yp = 2.876, R0, r = 1, xp = 0.08).This model can display a rich set of dynamics depending on parameter values31. Here we alter resource carrying capacity K in order to simulate the dynamics (using the deSolve package32 in R) of qualitatively different attractors (See Supplementary Fig. 1; K = 0.78, equilibrium; K = 0.85; two-point limit cycle; K = 0.92, four-point limit cycle; K = 0.997, chaotic dynamics). Every fifth time step of the simulated dynamics, corresponding to a sampling frequency of ≈10 samples per cycle for the 2-point limit cycle, was sampled. Observation noise was thereafter added to the deterministic dynamics produced by the model:$${N}_{l}(t)={N}_{l}^{prime}(t)+rho * e(t);e(t) sim N(0,{sigma }_{N^{prime_{l}}}),$$ (8) where (N_{l}^{prime}(t)) is the abundance of species l (P, C or R) simulated by the food-chain model at time point t, (rho) is the level of observation noise and ({sigma }_{N_{l}^{prime}}) is the standard deviation of the deterministic dynamics of species l produced by the food chain model.In order to investigate how time series length and observation noise affects the probability of detecting alternative attractors we derived probability landscapes. These were derived by testing the null-hypothesis (H0:(|{hat{{{{{{boldsymbol{Y}}}}}}}}_{{{{{{boldsymbol{i}}}}}}}|{{{{{{boldsymbol{M}}}}}}}_{{{{{{boldsymbol{i}}}}}}}-{{{{{{boldsymbol{Y}}}}}}}_{{{{{{boldsymbol{i}}}}}}}| > |{hat{{{{{{boldsymbol{Y}}}}}}}}_{{{{{{boldsymbol{i}}}}}}}|{{{{{{boldsymbol{M}}}}}}}_{{{{{{boldsymbol{j}}}}}}}-{{{{{{boldsymbol{Y}}}}}}}_{{{{{{boldsymbol{i}}}}}}}|); See Test statistic above) across 100 replicates for each combination of time series length and level of observation noise, (rho). Time-series length was varied from 10 to 100 in steps of 10, and observation noise, (rho), was varied from 0.01 to 0.3 in steps of 0.01, in total yielding 300 combinations of observation noise and time series length, for each combination of dynamical regimes i and j. Predator dynamics was used to predict consumer and resource dynamics using the multivariate approach described above (results for the cases where consumer or resource dynamics are used to predict the other species´ dynamics are presented in Supplementary Figs. 2, 3). All time series were standardized ((mu =0;sd=1)) prior testing for dynamical difference.Experimental data setThe experimental data set was given by Fussman et al.7. This data set contains 14 time series of a predator Brachionus calyciflorus and its prey Chlorella vulgaris derived from chemostat experiments. Time series for different dilution rates were produced by keeping the dilution rate fixed in different chemostats (Supplementary Figs. 3–11). Brachionus calyciflorus and Chlorella vulgaris time series were used to predict Chlorella vulgaris and Brachionus calyciflorus time series, respectively, using the multivariate approach described above. We tested for qualitative difference in the temporal dynamics across all time series, which were standardized ((mu =0;sd=1)) prior testing.Alternative stable state modelWe used a stochastic version of a well-known alternative stable state model4,33 to produce alternative stochastic dynamical regimes. The model is described by:$${{{{{rm{d}}}}}}x=left(xleft(1-frac{x}{{{{{{rm{K}}}}}}}right)+frac{c{x}^{2}}{1-{x}^{2}}right){{{{{rm{d}}}}}}t+sigma {{{{{rm{d}}}}}}w,$$
(9)
where K is the carrying capacity (here set to 11), c is a harvest rate, and σ (here set to 0.01) is the magnitude of noise which is described by a Wiener process (dw).The model was simulated for fixed harvest rates (c) assuming that the system state resides in either of its two basins of attraction. The initial value for the simulation was set to the equilibrium of the noise-free model skeleton for fixed harvest rates c, and σ is set low in order to avoid stochastic flips, so-called flickering, between alternative basins of attraction. Dynamics was integrated (Δt = 0.01) using the matlab-package SDE-Tools34.In order to investigate how time-series length and harvest rate, c, affects the probability of detecting alternative attractors in stochastic regimes we derived probability landscapes.These were derived by testing the null-hypothesis H0:(|{hat{{{{{{boldsymbol{Y}}}}}}}}_{{{{{{boldsymbol{i}}}}}}}|{{{{{{boldsymbol{M}}}}}}}_{{{{{{boldsymbol{i}}}}}}}-{{{{{{boldsymbol{Y}}}}}}}_{{{{{{boldsymbol{i}}}}}}}| > |{hat{{{{{{boldsymbol{Y}}}}}}}}_{{{{{{boldsymbol{i}}}}}}}|{{{{{{boldsymbol{M}}}}}}}_{{{{{{boldsymbol{j}}}}}}}-{{{{{{boldsymbol{Y}}}}}}}_{{{{{{boldsymbol{i}}}}}}}|) (permutation test p = 0.05) across 100 simulated data sets for each combination of time series length and harvest rate, c. Time-series length was varied between 50 and 150 in steps of 10, and c was varied between 1.83 and 2.73 in steps of 0.05, in total yielding 209 combinations of time series length and harvest rate. Each time series was standardized ((mu =0;sd=1)) prior testing for difference in temporal dynamics of contrasting regimes.Natural time-series dataIn a previous study on early warning signals of impending regime shifts, Gsell et al.18 used breakpoint analysis to identify two potential alternative dynamical regimes. We here test if these two-time series segments constitute alternative dynamical attractors. Prior analysis, we imputed a few missing observations (n = 24) using a kalman smoother35. The two time series segments, i.e. pre- and post-breakpoint time series, were standardized ((mu =0;sd=1)) prior testing for dynamical difference.Reporting summaryFurther information on research design is available in the Nature Research Reporting Summary linked to this article. More

1.Watson, J. E. M. et al. The exceptional value of intact forest ecosystems. Nat. Eco. Evo. 2, 599–610 (2018).Article 

Google Scholar 
2.European Commission. in COM(2020) 380 final (Brussels, 2020).3.Vandekerkhove, K. et al. Reappearance of Old-Growth Elements in Lowland Woodlands in Northern Belgium: Do the Associated Species Follow? Silva Fenn. 45, 909–935 (2011).Article 

Google Scholar 
4.Di Marco, M., Ferrier, S., Harwood, T. D., Hoskins, A. J. & Watson, J. E. Wilderness areas halve the extinction risk of terrestrial biodiversity. Nature 573, 582–585 (2019).ADS 
PubMed 
Article 
CAS 
PubMed Central 

Google Scholar 
5.Frey, S. J. K. et al. Spatial models reveal the microclimatic buffering capacity of old-growth forests. Sci. Adv. 2, e1501392 (2016).ADS 
PubMed 
PubMed Central 
Article 

Google Scholar 
6.Zhou, G. Y. et al. Old-growth forests can accumulate carbon in soils. Science 314, 1417–1417 (2006).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
7.Burrascano, S., Keeton, W. S., Sabatini, F. M. & Blasi, C. Commonality and variability in the structural attributes of moist temperate old-growth forests: A global review. For. Ecol. Manag. 291, 458–479 (2013).Article 

Google Scholar 
8.Bauhus, J., Puettmann, K. & Messier, C. Silviculture for old-growth attributes. For. Ecol. Manag. 258, 525–537 (2009).Article 

Google Scholar 
9.Moore, K. D. In the shadow of the cedars: the spiritual values of old-growth forests. Conserv. Biol. 21, 1120–1123 (2007).PubMed 
Article 
PubMed Central 

Google Scholar 
10.FOREST EUROPE. State of Europe’s Forests 2015. (Ministerial Conference on the Protection of Forests in Europe, Madrid, 2015).11.Ceccherini, G. et al. Abrupt increase in harvested forest area over Europe after 2015. Nature 583, 72–77 (2020).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
12.Levers, C. et al. Drivers of forest harvesting intensity patterns in Europe. For. Ecol. Manag. 315, 160–172 (2014).Article 

Google Scholar 
13.Potapov, P. et al. The last frontiers of wilderness: Tracking loss of intact forest landscapes from 2000 to 2013. Sci. Adv. 3 (2017).14.Schickhofer, M. & Schwarz, U. Inventory of Potential Primary and Old-Growth Forest Areas in Romania (PRIMOFARO). Identifying the largest intact forests in the temperate zone of the European Union. (Euronatur Foundation, 2019).15.Knorn, J. et al. Continued loss of temperate old-growth forests in the Romanian Carpathians despite an increasing protected area network. Environ. Conserv. 40, 182–193 (2013).Article 

Google Scholar 
16.Court of Justice of the European Union. C-441/17 – Commission v Poland (Forêt de Białowieża) Judgment of the Court (Grand Chamber) of 17 April 2018 (2018).17.Chylarecki, P. & Selva, N. Ancient forest: spare it from clearance. Nature 530, 419–419 (2016).ADS 
PubMed 
Article 
CAS 
PubMed Central 

Google Scholar 
18.Earthsight. Complicit in corruption. How billion-dollar firms and EU governments are failing Ukraine’s forests. (Earthsight, 2018).19.Mikoláš, M. et al. Primary forest distribution and representation in a Central European landscape: Results of a large-scale field-based census. For. Ecol. Manag. 449, 117466 (2019).Article 

Google Scholar 
20.Hance, J. IKEA Logging Old-growth Forest for Low-price Furniture in Russia. https://news.mongabay.com/2012/05/ikea-logging-old-growth-forest-for-low-price-furniture-in-russia/ (2012).21.Sabatini, F. M. et al. Protection gaps and restoration opportunities for primary forests in Europe. Divers. Distrib. 26, 1646–1662 (2020).Article 

Google Scholar 
22.Barredo Cano, J. I. et al. Mapping and assessment of primary and old-growth forests in Europe. (EUR 30661 EN, Publications Office of the European Union, 2021).23.Adam, D. & Vrška, T. Important localities of old-growth forests in Landscape Atlas of the Czech Republic (eds T Hrnčiarová, P Mackovčin, & I Zvara) (Ministry of Environment and Silva Tarouca Research Institute, Prague–Silva Tarouca Research Institute for Landscape and Ornamental Gardening, 2009).24.Diaci, J. Virgin forests and forest reserves in Central and East European countries-History, present status and future development. Proceedings of the invited lecturers’ reports presented at the COST E4 management committee and working groups meeting in Ljubljana, Slovenia (1999).25.Kirchmeir, H. & Kovarovics, A. Nomination Dossier “Primeval Beech Forests of the Carpathians and Other Regions of Europe“ as extension to the existing Natural World Heritage Site “Primeval Beech Forests of the Carpathians and the Ancient Beech Forests of Germany” (1133bis). (2016).26.García Feced, C., Berglund, H. & Strnad, M. Scoping document: information related to European old growth forests. (ETC/BD report to the EEA, 2015).27.Veen, P. et al. Virgin forests in Romania and Bulgaria: results of two national inventory projects and their implications for protection. Biodivers. Conserv. 19, 1805–1819 (2010).Article 

Google Scholar 
28.Ibisch, P. L. & Ursu, A. Potential primary forests of Romania. (Greenpeace CEE Romania; Centre for Econics and Ecosystem Management, Eberswalde University for Sustainable Development; Geography Department, A. I. Cuza University of Iași, 2017).29.Spracklen, B. D. & Spracklen, D. V. Identifying European Old-Growth Forests using Remote Sensing: A Study in the Ukrainian Carpathians. Forests 10, 127 (2019).Article 

Google Scholar 
30.Frank, G. et al. COST Action E27. Protected Forest Areas in Europe-analysis and harmonisation (PROFOR): results, conclusions and recommendations. (Federal Research and Training Centre for Forests, Natural Hazards and Landscape (BFW), 2007).31.Sabatini, F. M. et al. Where are Europe’s last primary forests? Divers. Distrib. 24, 1426–1439 (2018).Article 

Google Scholar 
32.McRoberts, R. E., Susanne, W., Gherardo, C. & Elizabeth, L. Assessing Forest Naturalness. For. Sci. 58, 294–309 (2012).Article 

Google Scholar 
33.Sabatini, F. M. et al. European Primary Forest Database. figshare https://doi.org/10.6084/m9.figshare.13194095.v1 (2020).34.FAO. Global Forest Resources Assessment 2015. Terms and definitions. (FAO, 2015).35.Buchwald, E. A hierarchical terminology for more or less natural forests in relation to sustainable management and biodiversity conservation in Proceedings: Third expert meeting on harmonizing forest-related definitions for use by various stakeholders (Food and Agriculture Organization of the United Nations, 2005).36.Blasi, C., Burrascano, S., Maturani, A. & Sabatini, F. M. Old-growth forests in Italy. (Palombi Editori, 2010).37.Cateau, E. et al. Le patrimoine forestier des réserves naturelles. Focus sur les forêts à caractère naturel. Cahier n°7. (Réserves Naturelles de France, 2017).38.Svoboda, M. et al. Landscape-level variability in historical disturbance in primary Picea abies mountain forests of the Eastern Carpathians, Romania. J. Veg. Sci. 25, 386–401 (2014).Article 

Google Scholar 
39.Potapov, P. et al. Mapping the world’s intact forest landscapes by remote sensing. Ecol. Soc. 13 (2008).40.Britz, H. et al. Nomination of the “Ancient Beech Forests of Germany” as Extension to the World Natural heritage “Primeval Beech Forests of the Carpathians”. Nationale Naturlandschaften, Federal Republic of Germany. Nieden-stein: Specialised editing Cognitio Kommunikation & Planung (2009).41.UNEP-WCMC & IUCN. Protected Area Profile for Ancient and Primeval Beech Forests of the Carpathians and Other Regions of Europe from the World Database of Protected Areas https://www.protectedplanet.net/903141 (2019).42.EEA. Biogeographical regions of Europe https://www.eea.europa.eu/data-and-maps/data/biogeographical-regions-europe-3 (2016).43.EEA. European forest types. Categories and types for sustainable forest management reporting and policy. (EEA Technical Report No 9/2006. EEA, Copenhagen, 2006).44.Bohn, U. et al. Map of the natural vegetation of Europe. Explanatory text with CD-ROM, (German Federal Agency for Nature Conservation, Bonn, Germany, 2003).45.Hansen, M. C. et al. High-resolution global maps of 21st-century forest cover change. Science 342, 850–853 (2013).ADS 
CAS 
Article 

Google Scholar 
46.R Development Core Team. R: A language and environment for statistical computing v. 3.6.1. R Foundation for Statistical Computing http://www.R-project.org/ (2019).47.Miljødirektoratet. (2016).48.FAO. Global Forest Resources Assessment 2015. Desk reference. 245 (FAO, Rome, 2015).49.FOREST EUROPE. Quantitative Indicators Country reports 2015 https://foresteurope.org/state-europes-forests-2015-report/#1476295965372-d3bb1dd0-e9a0 (2015).50.Cohen, W. B., Yang, Z., Healey, S. P., Kennedy, R. E. & Gorelick, N. A LandTrendr multispectral ensemble for forest disturbance detection. Remote Sens. Environ. 205, 131–140 (2018).ADS 
Article 

Google Scholar 
51.Kennedy, E. R. et al. Implementation of the LandTrendr Algorithm on Google Earth Engine. Remote Sensing 10 (2018).52.Gorelick, N. et al. Google Earth Engine: Planetary-scale geospatial analysis for everyone. Remote Sens. Environ. 202, 18–27 (2017).ADS 
Article 

Google Scholar 
53.Kennedy, R. E., Yang, Z. & Cohen, W. B. Detecting trends in forest disturbance and recovery using yearly Landsat time series: 1. LandTrendr—Temporal segmentation algorithms. Remote Sens. Environ. 114, 2897–2910 (2010).ADS 
Article 

Google Scholar 
54.Griffiths, P., Van Der Linden, S., Kuemmerle, T. & Hostert, P. A pixel-based landsat compositing algorithm for large area land cover mapping IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing 6, 2088–2101 (2013).55.Cohen, W. B. & Spies, T. A. Estimating structural attributes of Douglas-fir/western hemlock forest stands from Landsat and SPOT imagery. Remote Sens. Environ. 41, 1–17 (1992).ADS 
Article 

Google Scholar 
56.Czerwinski, C. J., King, D. J. & Mitchell, S. W. Mapping forest growth and decline in a temperate mixed forest using temporal trend analysis of Landsat imagery, 1987–2010. Remote Sens. Environ. 141, 188–200 (2014).ADS 
Article 

Google Scholar 
57.Cohen, W. B., Yang, Z. & Kennedy, R. Detecting trends in forest disturbance and recovery using yearly Landsat time series: 2. TimeSync—Tools for calibration and validation. Remote Sens. Environ. 114, 2911–2924 (2010).ADS 
Article 

Google Scholar 
58.Grogan, K., Pflugmacher, D., Hostert, P., Kennedy, R. & Fensholt, R. Cross-border forest disturbance and the role of natural rubber in mainland Southeast Asia using annual Landsat time series. Remote Sens. Environ. 169, 438–453 (2015).ADS 
Article 

Google Scholar 
59.De Marzo, T. et al. Characterizing forest disturbances across the Argentine Dry Chaco based on Landsat time series. International Journal of Applied Earth Observation and Geoinformation 98, 102310 (2021).Article 

Google Scholar 
60.Frank, A. Inventering av nyckelbiotoper: resultat till och med 2003. (Skogsstyr., 2004).61.Länsstyrelsen Västerbotten. LstAC Skogar med höga naturvärden ovan gränsen för fjällnära skog 2003–2015 https://ext-geodatakatalog.lansstyrelsen.se/GeodataKatalogen/ (2019).62.Naturvårdsverket. Skyddsvärda statliga skogar http://mdp.vic-metria.nu/miljodataportalen/GetMetaDataById?UUID=3919E66E-2E09-440D-9171-B5074DF0C0ED (2017).63.Naturvårdsverket. Skogliga värdekärnor http://gpt.vic-metria.nu/data/land/skogliga_vardekarnor_2016.zip (2016).64.Naturvårdsverket. Preciserad kartering av kontinuitetsskog i Jämtlands län http://gpt.vic-metria.nu/data/land/Preciserad_kskog_jamtland.zip (2019).65.Ahlkrona, E., Giljam, C. & Wennberg, S. Kartering av kontinuitetsskogi boreal region. Metria AB på uppdrag av Naturvårdsverket (2017).66.Naturvårdsverket. Skyddad fjallbarrskog https://gpt.vic-metria.nu/data/land/NMD/Skyddad_Fjallbarrskog.zip (2019).67.Trotsiuk, V. et al. A mixed severity disturbance regime in the primary Picea abies (L.) Karst. forests of the Ukrainian Carpathians. For. Ecol. Manag. 334, 144–153 (2014).Article 

Google Scholar 
68.Kozák, D. et al. Profile of tree-related microhabitats in European primary beech-dominated forests. For. Ecol. Manag. 429, 363–374 (2018).Article 

Google Scholar 
69.Garbarino, M. et al. Gap disturbances and regeneration patterns in a Bosnian old-growth forest: a multispectral remote sensing and ground-based approach. Ann. For. Sci. 69, 617–625 (2012).Article 

Google Scholar 
70.Keren, S. et al. Comparative Structural Dynamics of the Janj Mixed Old-Growth Mountain Forest in Bosnia and Herzegovina: Are Conifers in a Long-Term Decline? Forests 5, 1243–1266 (2014).Article 

Google Scholar 
71.Motta, R. et al. Structure, spatio-temporal dynamics and disturbance regime of the mixed beech–silver fir–Norway spruce old-growth forest of Biogradska Gora (Montenegro). Plant Biosyst. 149, 966–975 (2015).Article 

Google Scholar 
72.Motta, R. et al. Development of old-growth characteristics in uneven-aged forests of the Italian Alps. Eur. J. For. Res. 134, 19–31 (2015).Article 

Google Scholar 
73.Panayotov, M. et al. Mountain coniferous forests in Bulgaria – structure and natural dynamics. (University of Forestry and Geosoft, 2016).74.Lõhmus, A. & Kraut, A. Stand structure of hemiboreal old-growth forests: Characteristic features, variation among site types, and a comparison with FSC-certified mature stands in Estonia. For. Ecol. Manag. 260, 155–165 (2010).Article 

Google Scholar 
75.EEA. Developing a forest naturalness indicator for Europe. Concept and methodology for a high nature value (HNV) forest indicator. (EEA Technical report No 13/2014, Luxembourg: Publications Office of the European Union, 2014).76.Rossi, M., Bardin, P., Cateau, E. & Vallauri, D. Forêts anciennes de Méditerrané e et des montagnes limitrophes: références pour la naturalité régionale. WWF France, Marseille, France, 144 (2013).77.Myhre, T. Skogkur 2020. redningsplan for Norges unike skoger. WWF Verdens villmarksfond, Norges naturvernforbund, SABIMA (2012).78.Ruete, A., Snäll, T. & Jönsson, M. Dynamic anthropogenic edge effects on the distribution and diversity of fungi in fragmented old-growth forests. Ecol. Appl. 26, 1475–1485 (2016).PubMed 
Article 
PubMed Central 

Google Scholar 
79.Heiri, C., Wolf, A., Rohrer, L., Brang, P. & Bugmann, H. Successional pathways in Swiss mountain forest reserves. Eur. J. For. Res. 131, 503–518 (2012).Article 

Google Scholar 
80.Brang, P., Heiri, C. & Bugmann, H. Waldreservate: 50 Jahre natürliche Waldentwicklung in der Schweiz. (Haupt, 2011).81.Pantić, D. et al. Structural, production and dynamic characteristics of the strict forest reserve’Račanska šljivovica’on Mt. Tara. Glasnik Šumarskog fakulteta, 93–114 (2011).82.Savoie, J. M. et al. Vieilles forêts pyrénéennes de Midi-Pyrénées. Deuxième phase. Evaluation et cartographie des sites. Recommandations. Rapport final. (Ecole d’Ingénieurs de PURPAN/DREAL Midi-Pyrénées, 2015).83.Savoie, J. M. et al. Forêts pyrénéennes anciennes de Midi-Pyrénées. Rapport d’Etude de projet FEDER 2008–2011. 320 (Ecole d’Ingénieurs de PURPAN/DREAL Midi-Pyrénées, 2011).84.WWF Finland. Kansallisomaisuus turvaan – valtion omistamia suojelun arvoisia metsä- ja suoalueita, (WWF Suomen raportteja, 2012).85.Kitenberga, M. et al. A mixture of human and climatic effects shapes the 250-year long fire history of a semi-natural pine dominated landscape of Northern Latvia. For. Ecol. Manag. 441, 192–201 (2019).Article 

Google Scholar 
86.Baders, E., Senhofa, S., Purina, L. & Jansons, A. Natural succession of Norway spruce stands in hemiboreal forests: case study in Slitere national park, Latvia. Baltic Forestry 23, 522–528 (2017).
Google Scholar 
87.Kokarēviča, I. et al. Vegetation changes in boreo–nemoral forest stands depending on soil factors and past land use during an 80 year period of no human impact. Can. J. For. Res. 46, 376–386 (2016).Article 

Google Scholar 
88.Fernandez López, A. B. Parque Nacional de Garajonay, Patrimonio Mundial. (Organismo Autonomo Parques Nacionales, 2009).89.TRAGSATEC. Segundo inventario ecológico del Parque Nacional de Garajonay. (Parque Nacional de Garajonay, 2006).90.Fernández, A. B. & Gómez, L. Qué son los bosques antiguos de laurisilva. Su valor y situación en Canarias. La Gomera, entre bosques y taparuchas, 177–236 (2016).91.Matović, B. et al. Comparison of stand structure in managed and virgin european beech forests in Serbia. Šumarski list 142, 47–57 (2018).Article 

Google Scholar 
92.Kiš, A., Stojšić, V., & Dinić, A. In 2nd International Symposium on Nature Conservation. Proceedings 373–382 (Institute for Nature Conservation of Vojvodina Province, Novi Sad, 2016).93.Kobyakov, K. & Jakolev, J. Atlas of high conservation value areas, and analysis of gaps and representativeness of the protected area network in northwest Russia. (Finnish Environment Institute, 2013).94.Diku, A. & Shuka, L. Pyjet e vjetër të ahut në shqipëri (Old Beech forests in Albania). (PSEDA – ILIRIA, 2017).95.Burrascano, S. et al. It’s a long way to the top: Plant species diversity in the transition from managed to old-growth forests. J. Veg. Sci. 29, 98–109 (2018).Article 

Google Scholar  More

1.Zarzycki, K. Paprotniki i rośliny kwiatowe (rośliny naczyniowe). In: Flora i Fauna Pienin. (ed. Razowski J). Monogr. Pienińskie 1, 75–79 (2000).
Google Scholar 
2.Środoń, W. Pieniny w historii szaty roślinnej Podhala [Pieniny in the history of plant cover in Podhale region]. In : K. Zarzycki (ed.). Przyroda Pienin w obliczu zmian [The nature of the Pieniny Mts (West Carpathians) in face of the coming changes]. Stud. Nat. 30B, 115–126 (1982).
Google Scholar 
3.Deptuła, C. Nad rekonstrukcją dziejów regionu czartoryskiego w XIII I XIV wieku [On the reconstruction of the history of the Czorsztyn region from the 13th to 16th centuries]. Pieniny—Człowiek Przyroda 5, 21–35 (1997) (in Polish with English summary).
Google Scholar 
4.Kierś, M. (ed.) Wołosi: Nomadzi Bałkanów (Uniwersytet Jagielloński, 2013).
Google Scholar 
5.Oravcová, M. & Krupa, E. Pedigree analysis of the former Valachian sheep. Slovak. J. Anim. Sci. 44, 6–12 (2011).
Google Scholar 
6.Wace, A.J.B. & Thompson, M.S. The Nomads of the Balkans. Vol. 6 (Methuen & Co., 1914). https://archive.org/stream/nomadsofbalkansa00wace#page/n9/mode/2up. Accessed 28 June 2021.7.Stachurska-Swakoń, A. Phytogeographical aspects of grasses occuring in tall-herb vegetation in the Carpathians. in Grasses in Poland and Elsewhere (ed. Frey, L.). 39–47. (W. Szafer Institute of Botany, Polish Academy of Sciences, 2009).8.Stachurska-Swakoń, A. Syntaxonomical revision of the communities with Rumex alpinus L. in the Carpathians. Phytocoenologia 39, 217–234. https://doi.org/10.1127/0340-269X/2009/0039-0217 (2009).Article 

Google Scholar 
9.Ralska-Jasiewiczowa, M., Nalepka, D. & Goslar, T. Some problems of forest transformation at the transition to the oligocratic/Homo sapiens phase of the Holocene interglacial in northern lowlands of central Europe. Veg. Hist. Archaeobot. 12, 233–247. https://doi.org/10.1007/s00334-003-0021-8 (2003).Article 

Google Scholar 
10.Pawłowski, B., Pawłowska, S. & Zarzycki, K. Zespoły roślinne kośnych łąk północnej części Tatr i Podtatrza. Fragm. Flor. Geobot. Pol. 6(2), 95–222 (1960).
Google Scholar 
11.Korzeniak, J. 6520* Mountain Yellow Trisetum and Bent-Grass Hay Meadows 55–67 (Methodological guide. GIOŚ, 2013).
Google Scholar 
12.Wróbel, I. Pasterstwo w regionie pienińskim [Sheep farming in the Pieniny region]. Pieniny Człowiek Przyroda 5, 43–52 (1997) (in Polish with English summary).
Google Scholar 
13.Kostrakiewicz-Gierałt, K., Palic, C. C., Stachurska-Swakoń, A., Nedeff, V. & Sandu, I. The causes of disappearance of sward lily Gladiolus imbricatus L from natural stands—Synthesis of current state of knowledge. Int. J. Conserv. Sci. 9, 821–834 (2018).
Google Scholar 
14.Wróbel, I. Szata roślinna Pienińskiego Parku Narodowego – podsumowanie Planu Ochrony na lata 2001–2020 [Plant cover of the Pieniny National Park – summing up the Protection Plan for the years 2001–2020]. Pieniny Człowiek Przyroda 8, 63–69 (2003).
Google Scholar 
15.Kubíková, P. & Zeidler, M. Habitat demands and population characteristics of the rare plant species Gladiolus imbricatus L. in the Frenštát region (NE Moravia, the Czech Republic). Čas. Slez. Muz. Opava 60(A), 154–164 (2011).
Google Scholar 
16.Mirek, Z., Piękoś-Mirkowa, H., Zając, A. & Zając, M. Flowering Plants and Pteridophytes of Poland, a Checklist (W. Szafer Institute of Botany, Polish Academy of Sciences, 2002).
Google Scholar 
17.Hamilton, A. P. The European Gladioli. Quart. Bull. Alp. Gard. Soc. 44, 140–146 (1976).
Google Scholar 
18.Kornaś, J. M. & Medwecka-Kornaś, A. Zespoły roślinne Gorców. I. Naturalne i na wpół naturalne zespoły nieleśne. Fragm. Flor. Geobot. Polon. 13(2), 167–316 (1967).
Google Scholar 
19.Ascherson, P. & Engler, A. Beiträge zur Flora Westgaliziens und der Central-Karpaten. Osterr. Bot. Z. 15, 273–285. https://doi.org/10.1007/BF01623075 (1865).Article 

Google Scholar 
20.Wołoszczak, E. Zapiski botaniczne z Karpat Sądeckich. Spraw. Komis. Fizjogr. AU 30, 174–206 (1895).
Google Scholar 
21.Zapałowicz, H. Conspectus Florae Galiciae Criticus Vol. 1 (Nakł. Akad. Umiej., 1906).
Google Scholar 
22.Piękoś-Mirkowa, H. & Mirek, Z. Flora Polski. Rośliny Chronione (Oficyna Wydawnicza Multico, 2006).
Google Scholar 
23.Dembicz, I. et al. New locality of Trollius europaeus L. and Gladiolus imbricatus L. near Sochocin by Płońsk (Central Poland). Opole Sci. Soc. Nat. J. 44, 36–46 (2011).
Google Scholar 
24.Kropač, Z. & Mochnacký, S. Contribution to the segetal communities of Slovakia, Thaiszia. J. Bot. 19, 145–211 (2009).
Google Scholar 
25.Mirek, Z., Nikel, A. & Wilk, Ł. Ozdoba łąk reglowych. Tatry 4(50), 50–51 (2014).
Google Scholar 
26.Kołos, A. A new locality of Gladiolus imbricatus (Iridaceae) in the North Podlasie Lowland. Fragm. Florist. Geobot. Polon. 22(2), 390–395 (2015).
Google Scholar 
27.Falkowski, M. Nowe stanowisko Gladiolus imbricatus (Iridaceae) w dolinie środkowej Wisły. Fragm. Florist. Geobot. Polon. 9, 369–370 (2002).
Google Scholar 
28.Nowak, A. & Antonin, A. Interesujące stanowiska Gladiolus imbricatus (Iridaceae) w Bramie Morawskiej. Fragm. Florist. Geobot. Polon. 13(1), 17–22 (2006).
Google Scholar 
29.Stepansky, A., Kovalski, I. & Perl-Treves, R. Interspecific classification of melons (Cucumis melo L.) in view of their phenotypic and molecular variation. Plant Syst. Evol. 271, 313–332. https://doi.org/10.1007/BF00984373 (1999).Article 

Google Scholar 
30.Gupta, M., Chyi, Y.-S., Romero-Sverson, J. & Owen. J.L. Amplification of DNA markers from evolutionarily diverse genomes using single primers of simple-sequence repeats. Theor. Appl. Genet. 89, 998–1006. https://doi.org/10.1007/BF00224530 (1994).31.Sutkowska, A., Pasierbiński, A., Warzecha, T., Mandal, A. & Mitka, J. Refugial pattern of Bromus erectus in Central Europe based on ISSR fingerprinting. Acta Biol. Cracov. Ser. Bot. 55(2), 107–119. https://doi.org/10.2478/abcsb-2013-0026 (2013).Article 

Google Scholar 
32.Bonin, A. et al. How to track and assess genotyping errors in population genetic studies. Mol. Ecol. 3, 3261–3273. https://doi.org/10.1111/j.1365-294X.2004.02346.x (2004).CAS 
Article 

Google Scholar 
33.Vekemans, X. AFLP-surv 1.0: A Program for Genetic Diversity Analysis with AFLP (and RAPD) Population Data. https://ebe.ulb.ac.be/ebe/AFLP-SURV.html (Laboratoire de Génétique et d’Ecologie Végétales, Université Libre de Bruxelles, 2002).34.Yeh, F., Yang, R. & Boyle, T. POPGENE Version 1.32. Microsoft-Based Freeware for Population Genetic Analysis. https://www.softpedia.com/get/Science-CAD/Popgene-Population-Genetic-Analysis.shtml (Molecular Biology and Biotechnology Center, University of Alberta, 1999).35.Schönswetter, P. & Tribsch, A. Vicariance and dispersal in the Alpine perennial Bupleurum stellatum L (Apiaceae). Taxon 54, 725–732. https://doi.org/10.2307/25065429 (2005).Article 

Google Scholar 
36.Ehrich, D. AFLPdat: A collection of r functions for convenient handling of AFLP data. Mol. Ecol. Notes 6, 603–604. https://doi.org/10.1111/j.1471-8286.2006.01380.x. https://mybiosoftware.com/tag/aflpdat (2006).37.Paun, O., Schönswetter, P. & Winkler, M., Intrabiodiv Consortium & Tribsch, A. Historical divergence versus contemporary gene flow: Evolutionary history of the calcicole Ranunculus alpestris group (Ranunculaceae) in the European Alps and the Carpathians. Mol. Ecol. 17, 4263–4275. https://doi.org/10.1111/j.1365-294x.2008.03908.x (2008).38.Excoffier, L., Smouse, P. E. & Quattro, J. M. Analysis of molecular variance inferred from metric distances among DNA haplotypes: Application to human mitochondrial DNA restriction data. Genetics 131, 479–491. https://doi.org/10.1093/genetics/131.2.479 (1992).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
39.Excoffier, L., Laval, G. & Schneider, S. Arlequin (version 3.0): An integrated software package for population genetics data analysis. Evol. Biol. 1, 47–50. http://cmpg.unibe.ch/software/arlequin3/. https://doi.org/10.1177/117693430500100003 (2005).40.Lynch, M. & Milligan, B. Analysis of population-genetic structure using RAPD markers. Mol. Ecol. 3, 91–99. http://cmpg.unibe.ch/software/arlequin3/. https://doi.org/10.1111/j.1365-294x.1994.tb00109.x (1994).41.Saitou, N. & Nei, M. The neighbour-joining method: A new method for reconstructing phylogenetic trees. Mol. Biol. Evol. 4, 406–425. https://doi.org/10.1093/oxfordjournals.molbev.a040454 (1987).CAS 
Article 
PubMed 

Google Scholar 
42.Makarenkov, V. T-Rex: Reconstructing and visualizing phylogenetic trees and reticulation networks. Bioinformatics 17, 664–668. http://www.fas.umontreal.ca/biol/casgrain/en/labo/t-rex. https://doi.org/10.1093/bioinformatics/17.7.664 (2001).43.Makarenkov, V. & Legendre, P. The fitting of a tree metric to a given dissimilarity with the weighted least squares criterion. J. Classif. 16, 3–26. https://doi.org/10.1007/s003579900040 (1999).Article 

Google Scholar 
44.Felsenstein, J. Phylip (Phylogeny Inference Package) Version 3.6. https://evolution.genetics.washington.edu/phylip.html (University of Washington, 2005).45.Nei, M. & Li, W. H. Mathematical model for studying genetic variation in terms of restriction endonucleases. Proc. Natl. Acad. Sci. USA 76, 5269–5273. https://doi.org/10.1073/pnas.76.10.5269 (1979).ADS 
CAS 
Article 
PubMed 
PubMed Central 
MATH 

Google Scholar 
46.Kruskal, J. B. Nonmetric multidimensional scaling: A new numerical method. Psychometrika 29, 115–129 (1964).MathSciNet 
Article 

Google Scholar 
47.Rohlf, F. J. NTSYS-pc. Numerical Taxonomy and Multivariate Analysis, Version 2.1. https://ntsyspc.software.informer.com/ (Exeter Software, 2002).48.Pritchard, J. K, Stephens, M. & Donelly P. Inference of population structure using multilocus genotype data. Genetics 155, 945–959. http://web.stanford.edu/group/pritchardlab/structure.html (2000).49.Falush, D., Stephens, M. & Pritchard, J. K. Inference of population structure using multilocus genotype data: Dominant markers and null alleles. Mol. Ecol. Notes 7, 574–578. https://doi.org/10.1111/j.1471-8286.2007.01758.x (2007).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
50.Evanno, G., Regnaut, S. & Goudet, J. Detecting the number of clusters of individuals using the software structure: A simulation study. Mol. Ecol. 14, 2611–2620. https://doi.org/10.1111/j.1365-294X.2005.02553.x (2005).CAS 
Article 
PubMed 

Google Scholar 
51.Nordborg, M., Hu, T. T., Ishino, Y., Jhaveri, J. & Toomajian, C. The pattern of polymorphism in Arabidopsis thaliana. PLOS Biol. 3(7), e196. https://doi.org/10.1371/journal.pbio.0030196 (2005).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
52.Dybova-Jachowicz, S. & Sadowska, A. (eds) Palinologia (Inst. Botaniki im. W. Szafera, Polska Akademia Nauk, 2003).
Google Scholar 
53.Cieślak, E., Szczepaniak, M., Kamiński, R. & Heine, W. Stan zachowania krytycznie zagrożonego gatunku Gladiolus paluster (Iridaceae) w Polsce – Analiza zmienności genetycznej osobników w uprawie Ogrodu Botanicznego Uniwersytetu Wrocławskiego w kontekście prowadzonych działań ochronnych. Fragm. Florist. Geobot. Polon. 21(1), 49–66 (2014).
Google Scholar 
54.Kutlunina, N., Permyakova, M. & Belyaev, A. Genetic diversity and reproductive traits in triploid and tetraploid populations of Gladiolus tenuis (Iridaceae). Plant Syst. Evol. 303, 1–10. https://doi.org/10.1007/s00606-016-1347-x (2017).Article 

Google Scholar 
55.Sutkowska, A., Pasierbiński, A., Warzecha, T. & Mitka, J. Multiple cryptic refugia of forest grass Bromus benekenii in Europe as revealed by ISSR fingerprinting and species distribution modelling. Plant Syst. Evol. 300, 1437–1452. https://doi.org/10.2478/abcsb-2013-0026 (2014).Article 

Google Scholar 
56.Gajewski, Z, Boroń, P, Lenart-Boroń, A, Nowak, B., Sitek, E. & Mitka, J. Conservation of Primula farinosa in Poland with respect to the genetic structure of populations. Acta Soc. Bot. Pol. 87(2), 3577 (2018). https://doi.org/10.5586/asbp.3577.Article 

Google Scholar 
57.Stojak, J., McDevitt, A. D., Herman, J. S., Searle, J. B. & Wójcik, J. M. Post-glacial colonization of eastern Europe from the Carpathian refugium: evidence from mitochondrial DNA of the common vole Microtus arvalis. Biol. J. Linn. Soc. 115, 927–939. https://doi.org/10.1111/bij.1253541 (2015).Article 

Google Scholar 
58.Szczepaniak, M. & Cieślak, E. Low level of genetic variation within Melica transsilvanica populations from the Kraków-Częstochowa Upland and the Pieniny Mts revealed by AFLPs analysis. Acta Soc. Bot. Pol. 76(4), 321–331. https://doi.org/10.5586/asbp.2007.036 (2007).Article 

Google Scholar 
59.Bennett, K. D. & Provan, J. What do we mean by ‘refugia’?. Quatern. Sci. Rev. 27, 27–28. https://doi.org/10.1016/j.quascirev.2008.08.019 (2008).Article 

Google Scholar 
60.Petit, R. J. et al. Glacial refugia: Hotspots but not melting pots of genetic diversity. Science 300(5625), 1563–1565. https://doi.org/10.1126/science.1083264 (2003).ADS 
CAS 
Article 

Google Scholar 
61.Brus, R. Growing evidence for the existence of glacial refugia of European beech (Fagus sylvatica L.) in the south-eastern Alps and north-western Dinaric Alps. Periodicum Biol. 112, 239–246 (2010).
Google Scholar 
62.Jŏgar, Ü. & Moora, M. Reintroduction of a rare plant (Gladiolus imbricatus) population to a river floodplain—How important is meadow management?. Restor. Ecol. 16, 382–385. https://doi.org/10.1111/j.1526-100X.2008.00435.x (2008).Article 

Google Scholar 
63.Mitka, J., Boroń, P., Wróblewska, A. & Bąba, W. AFLP analysis reveals intraspecific phylogenetic relationships and population genetic structure of two species of Aconitum in Central Europe. Acta Soc. Bot. Pol. 84(2), 267–276. https://doi.org/10.5586/asbp.2015.012 (2015).CAS 
Article 

Google Scholar 
64.Biernacka, M. Dawne oraz współczesne formy organizacji pasterstwa w Bieszczadach. Etnogr. Polska 6, 41–61. http://webcache.googleusercontent.com/search?q=cache:JDjzqMdApxIJ:cyfrowaetnografia.pl/Content/454+&cd=1&hl=pl&ct=clnk&gl=pl (1962).65.Stachurska-Swakoń, A., Cieślak, E. & Ronikier, M. Phylogeography of subalpine tall-herb species in Central Europe: the case of Cicerbita alpina. Preslia 84, 121–140. https://doi.org/10.1111/j.1095-8339.2012.01323.x (2012).Article 

Google Scholar  More

When introducing ARs as a fisheries management tool to Senegal, the Japanese management had the mindset of Japanese stakeholders, i.e., introducing fishing rights. However, after discussions with Senegalese stakeholders, it was decided that no-take areas would be delineated around ARs because the establishment of a strong fishing rights regime was not socially acceptable to the Senegalese fishing community. Japanese governance is based on the acceptance and respect of fishers towards individual, private AR concessions. In contrast, fishers in Senegal, and more widely in West Africa, are characterized by high mobility, particularly in the context of climate change and overexploitation18,19. Consequently, respect for local management regulations is lower, with open access being generally assumed. The basic concept of implementing a no-take area on the AR was not easily accepted by fishers. The immersion of AR concrete blocks was set as a top priority by managers at the expense of more complex socio-economic considerations, such as consciousness-raising activities and self-sustaining participative monitoring of the AR.The clear contradiction between the ecological knowledge of fishers and their behavior was explained by the well-known effects of open access resources on individual behavior. This phenomenon was also observed in our mathematical model. The processes in the mathematical model are in accordance with those perceived by the fishers, so that the results are also those expected by fisher’s local ecological knowledge. It is interesting to notice that the theoretical results presented here are the mathematical solutions of the model at equilibrium between fishing effort and fish population growth, i.e. after an oscillation period. It is obvious that short-term effect of fishing on the AR is always to increase the catch, but many fishers did perceive the longer-term effect of decreasing catches. The potential negative effect of the AR on catch when there is high fish attraction combined with high fishing pressure on the AR might explain the reluctance of a part of the fishers community to AR deployment (Fig. 2). In particular, the model illustrates that the AR attraction effect strongly determines the impact of the management. In general, fish attraction is the most immediate effect perceived after AR deployment11, as was true for our study16. Though the AR volume was relatively small (70 m3), the empty space between the higher blocks also contributes approximately 280 to 570 m3 of good habitat/refuge for schooling fish; therefore, it is actually difficult to accurately describe the volume that affects fish. Thus, it is difficult to say whether this AR is below or above the forecasted optimal volume in absence of fishing (120m3 with model parameters). The existence of an optimal volume for AR was also suggested by field studies as a trade off between food supply and refuge20, in line with our results. For management purposes, it is interesting to determine whether the AR is above or below this optimal level because if the volume is too small, the model predicted that any level of fishing on the AR would, in the long term, decrease the catch in the considered area. On the other hand, if the volume is above the optimal level, a small fishing effort on the AR could be authorized and would increase the total catch in the area.Field observation showed that the fish attraction effect was strong16 but precise estimation of this parameter cannot be inferred, as this would need, ideally, individual fish trajectories. Future field research on the attraction effect may permit estimating the AR attraction parameters. The model sensitivity test showed that the stronger the attraction parameter, the better the impact of the AR for the fisheries in case of no or small fishing effort on the AR (Fig. 3). But at the same time, the attraction is a strong incentive for fishers to fish on the AR, and the predicted benefit for fisheries in the fishing area rapidly vanishes when fishing effort on AR increases. This in turn provides further incentive for fishers to fish the AR, challenging the surveillance capacity. If fish attractiveness is strong and too many fishers fish on the AR, catch in the area will be concentrated on the AR, while the adjacent fishing area will be depleted, with catch levels lower than those prior to AR deployment.Specifically, in the context of generalized overfishing in Senegal21, deciding not to fish on the AR represents significant individual loss, despite being recognized as beneficial, globally22. It has been argued that this situation would rarely occur in small-scale fisheries, due to existing arrangements between individuals23. However, in the context of the highly mobile Senegalese artisanal fishing fleet and its overcapacity, as soon as the AR in Yenne was no longer subject to surveillance, it rapidly attracted fishers from other villages. Also, pre-existing arrangements between fishers might be overruled when new ARs are created, changing the structure of existing fishing grounds.At the time of the survey, the surveillance system set up by the co-management entities was not operational in our case study, because it was dependent on temporally limited external financing. These limitations are typical of short-term projects that focus on a single restricted area for a pre-determined duration, usually up to two years (e.g., NGOs, World Bank). Local fishers perceptions were globally in line with the model prediction that this AR fails to improve fisheries yield when surveillance is not in place to ensure AR regulations are observed, despite effective fish attraction and production existing in the AR.The model predicted that enhanced production on ARs could not keep pace with unrestricted access, which might be particularly true in Senegal where fishing effort rapidly reorganizes itself according to local yields24. Enhanced production due to the AR largely increases the catch if the fishing pressure on the AR remains null or very low, but it has no effect on the catch for higher fishing pressures on the AR (Fig. 3). These results were stable even if fish population growth, fish catchability, mobility and economic parameters could modulate the predicted amplitude of the catch and AR optimal volume. These results are consistent with existing theoretical studies of the impact of fisher movement to high production areas in and around MPAs25. Taking into account several species and their interactions (predation, competition) would lead to a very complex ecosystem model specific to the area (e.g. 26), with necessarily more assumptions. This model would necessarily be more difficult to share with fishers and other stakeholders. Both to simplify model structure and facilitate communication of results to stakeholders, we assumed in our model that the balance of entries exits and is in equilibrium, so that the migratory species did not affect the long-term equilibrium between fishing effort and fish abundance.The design of ARs could be adjusted to reduce the effect of illegal fishing by passively preventing both industrial and artisanal fishing activity. Complex structures are more effective for fish production and attraction27. We showed that, although production might have a limited effect on total catch, attraction can largely increase AR efficiency (total catch) if the rate of illegal fishing rate is very low or absent. Complex structures protect fish more effectively from small scale fishing gear28, including divers (Pers. Comm., Mamadou Sarr, Ouakam fishers committee). Thus, ARs should be appropriately designed to help mitigate potential issues28. Such designs might be more costly, and do not exclude the need for surveillance, but would enhance fisheries management, especially when surveillance cannot capture low levels of illegal fishing.Finally, if socio-economic and governance conditions are not met, well-intentioned AR projects will likely disturb the existing equilibrium among fishers that have different levels of access to the AR. Poor governance of marine resources has previously been described in West Africa, particularly in Senegal29, as has the failure of AR projects in a number of other developing countries9, which further deteriorate fishers trust and management plans efficiency30. In order to avoid that, NGO and governmental agencies driving ARs projects must consider that AR management induces collective costs before providing potentially collective gains. Thus, co-management that involves governmental institutions and fisher communities is required. Future management and adaptation plans for fishers, particularly in developing countries, should, therefore, focus efforts on raising long-term awareness of actors in both government institutions and fishing communities. At the level of institutional or development partners, long-term management costs should be included in the set-up of AR projects. For example, the local fishers committee of Yenne recently reported the establishment of a collective ship chandler whose profits are used to finance AR surveillance during the daytime. Subsequently, fishers noted an improvement in catches around the AR, even though illegal fishing likely continues on the AR at night (Pers. Comm. chair of local fishers committee). These observations support model predictions that low levels of illegal fishing might not disturb the positive impact of the AR. Alternatively surveillance effort could be supported by the community if benefits were managed according to ancestral traditions. Indeed, “no take area” regime on the AR would be in line with some past West African tribal laws, applied before the colonization era, which set marine area where fishing activities were restricted for occasional community celebrations. Collective processes where fishers and other stakeholders can design temporary no-take zones around the AR could increase fishers trust and compliance to the rules, fostering a positive socio-ecological feedback loop30.Hybridization of local and scientific knowledge, through the integration of natural sciences and social sciences, is key point for governance setting31,32,33. Indeed, the communication of the resulting hybrid knowledge in specific events gathering local stakeholders helps strengthen fisheries co-management for the establishment of surveillance and regulatory frameworks. This phenomenon was experienced during the public restitution of the present study with the community, fishers, children’s from local schools and governmental stakeholders. Science popularization of the study results was in French and local language (Wolof) retransmitted on national news (available at https://www.youtube.com/watch?v=yQqFU2P4XZU). Posters were exposed during the event, including pictures of local fishers interviewed and statements reflecting their own perception of how the artificial reef interacts with ecological processes and fisheries dynamics. Straightaway, stakeholders and local promoters of AR publicly expressed their concern and willingness to prioritize the setting up an efficient AR surveillance independent from external resources prior to increase AR deployments. Knowledge hybridization could produce more specific models that could be used for warning and advice, for example by considering potential impacts of ARs on species compositions3,34,35, environmental parameters36, and cascade effects on the trophic food web37. However this approach would need to be adapted to local social-ecological governance, which might require dedicated political-anthropological studies (see concept of adaptive co-management32).In summary, best practices should involve all stakeholders, consider local specificities, such as site configuration, governance, ecosystem, availability of ad hoc human and financial resources for AR surveillance, and define AR volume and design accordingly to these parameters. Thus, if plans exist to deploy ARs at large scales we recommend that legislation is strengthened, with detailed Environmental and social Impact Assessments38 to implement ARs, including considerations of long-term governance. More

Bagworm walking method using a ladder-like silk footholdWhen bagworms are reared in a plastic or glass cage, they walk not only on the floor but also on the walls or ceiling using only their three pairs of thoracic legs. The method by which they achieve this was clarified by placing a bagworm on black paper. Where the bagworm had walked, a ladder-like silk trace was observed on the black paper (Fig. 2a). Scanning electron microscopy (SEM) observation of one of the steps (or rungs) of the ladder-like trace revealed that each step was made up of a zigzag pattern of silk threads (Fig. 2b). Further magnified SEM observations revealed that the folded parts of the zigzag-spun thread were glued selectively to the substrate with adhesive whereas the remaining straight parts (hereafter, termed ‘bridges’ or ‘bridge threads’) were unglued (Fig. 2c–e).Figure 2Architecture of the ladder-like foothold. (a) A typical ladder-like foothold constructed by a bagworm on black paper, (b) an enlarged image showing one of the steps in the foothold and (c) a scanning electron microscopy image of the step shown in (b). The unglued bridge threads and a glued turn in the step shown in (c) are magnified in (d) and (e), respectively. (f) An enlarged image of four continuous steps in the foothold shown in (a). The neighbouring steps are connected via a single thread indicated by the arrows. (g) A schematic depiction of the basic architecture of the foothold; blue lines and green circles correspond to the silk thread and glued parts, respectively. (h) A photograph of a bagworm constructing a foothold on a transparent plastic board.Full size imageNotably, the steps of the foothold were not independent but rather always connected with neighbouring steps via a single thread (Fig. 2f). The overall basic construction of the foothold is schematically depicted in Fig. 2g. We found that the foothold was constructed in one continuous movement and always made of a single thread regardless of walking distance or time; therefore, a continuous thread exceeding a length of 100 m could be collected from one foothold14. We also observed bagworm climbing behaviour on a transparent plastic board, which clarified the important role of the silk trace as a foothold (Fig. 2h). During this behaviour, the bagworm used its sickle claws (Fig. 1e) to hook its second and third pairs of thoracic legs onto the first and second newest steps, respectively, and constructed the next step by spinning silk with a zigzag motion of the head and the skilful use of the first pair of thoracic legs. When the bagworm advanced one step, it always first shifted its third pair of thoracic legs to the next step before then shifting its second pair of thoracic legs to the newest step to avoid overloading this step, which may not yet be fully adhered to the surface (see Supplementary Movie S1). Because of this construction method, the interval distance between neighbouring steps is automatically determined by the interval between the thoracic legs. By repeating this process, the bagworm can advance forward slowly but steadily. This walking method was commonly observed on a horizontal floor surface, vertical wall, or horizontal ceiling. Although we have mainly described and shown observations from E. variegate here, with the exception of Supplementary Fig. S4 and Movie S1, we also observed instances of walking behaviour in other species, namely Eumeta minuscula, Mahasena aurea, Nipponopsyche fuscescens and Bambalina sp. (for a movie on E. minuscula walking behaviour, wherein it climbs a vertical wall, see Supplementary Movie S2). For at least 100 individuals of these bagworm species, we observed essentially identical walking behaviour to that described in the present study without exceptions for locomotion on substrates with slippery surfaces.Based on our observations, we asked the following question: how do bagworms selectively glue the folded parts of the foothold onto the substrate? Real-time observation of the tip of the spinneret (i.e. the spigot) through a transparent plastic board during the construction of the foothold revealed that adhesive was selectively discharged to attach the folded parts to the substrate; this process could be distinguished from the continuous spinning of the silk thread (for a movie showing construction behaviour, see Supplementary Movie S3). Figure 3a–g shows a time-sequence of foothold construction with enlarged images in the vicinity of the spinneret provided, whereas Fig. 3h depicts a schematic trace of the construction process. It was clearly noted that the bagworm discharged the adhesive only at the folded parts (shown in Fig. 3a–c,e,f; termed the ‘glued turn’) and not at the straight bridge parts (shown in Fig. 3d,g; termed the ‘unglued bridge thread’). From these time-sequence observations, we concluded that the bagworm controls the discharge of adhesive in an ‘on and off’ manner as necessary (essentially the same construction behaviours were confirmed for at least 20 individuals).Figure 3Foothold construction. (a–g) (left side) Time-sequence images taken during foothold construction and (right side) enlarged images of the vicinity of the spinneret (corresponding to the yellow rectangular area in each left-side image). The time-sequence images correspond to the parts of the schematic trace of foothold construction depicted by the red line in (h). In each right-side image and the schematic trace, the part of silk thread at which the adhesive was discharged is traced with a light-blue line. Green arrows in the right-side images show the direction of travel of the spinneret.Full size imagePassages of fibroin brins and adhesiveWe next investigated the spinning mechanism that enables continuous spinning of silk thread together with the selective discharge of adhesive via a single spigot. To this end, we observed the morphology of the bagworm from the silk gland to the spigot. Figure 4a shows the area in the vicinity of the spinneret, dissected and isolated from an E. variegata bagworm, which included a pair of silk glands and plural adhesive glands. As we previously reported21, the exterior shape of the silk gland in E. variegata (see Supplementary Fig. S1) is almost the same shape as that in the silkworm Bombyx mori and it is subdivided into three parts: the anterior (ASG), middle (MSG) and posterior (PSG) silk glands. We also previously confirmed that fibroin heavy chain (h-fib), fibroin light chain (l-fib) and fiboinhexamerin genes are expressed dominantly in the PSG, while sericin is expressed in the MSG, which strongly suggests that division-selective production of each protein exists in E. variegata (as has been shown in B. mori22). Figure 4b shows a magnified image of the spinneret including the end of the ASG. Beyond the pair of ASGs, which are merged into a common tube, a silk press and spinning tube appear before the spigot. This basic passage of silk fibroin from the ASG to the spigot is essentially the same as the passage observed in B. mori23. However, more detailed morphological observations of the inner structure of the passage revealed several obvious differences between E. variegata and B. mori.Figure 4Structural examination of the passages of fibroin brins and adhesive. (a) An optical microscope image of the area in the vicinity of a spinneret isolated from a female bagworm in the final instar stage. Indicated by arrows is a pair of silk glands (SG), one of the adhesive glands (ADG) and the spinneret (SP). (b) An optical microscope image of the passage including the (1) end of the anterior SGs (ASGs), (2) common tube, (3) silk press, (4) spinning tube and (5) spigot. (c–j) Optical microscope images showing cross-sections of the passage of fibroin brins obtained from the corresponding positions (c–j) in image (b). To focus on the fibroin brins and its passage, the surrounding outer part was removed so that a pair of fibroin brins was revealed in each image (except for image (c), which shows only one side of the ASG). Unmagnified images of (f–j), including the outer part, are shown in Supplementary Fig. S2. (k–n) 3D X-ray CT images of the spinneret: (k) overview, (l) cross-sectional top view, (m) cross-sectional side view and (n) passage of the fibroin brins and corresponding cross-sectional images at various positions. In the cross-sectional side view (m), the sheath and core parts are coloured blue and pink, respectively. (o) Image of the tip of a spigot from which adhesive is overflowing and a silk thread is emerging.Full size imageCross-sectional images along the spinneret are shown in Fig. 4c–j; these focus on the silk brins and their passage (unmagnified versions of the images in Fig. 4f–j are shown in Supplementary Fig. S2). The fibroin brins have an approximately round cross-sectional shape at the end of the ASG (Fig. 4c) and are merged at a common tube, which deforms their round shape slightly (Fig. 4d). The fibroin brins seem to be coated with a thin layer of sericin after the MSG, similar to B. mori; however, we omit the presence of the sericin layer here for convenience. The paired brins are gradually pressed between the ventral and dorsal hard cuticle plates at the silk press, and a gradual diameter decrease and shape deformation follows (Fig. 4e,f). At the exit of the silk press, each brin becomes elliptic and the diameter in the major axis decreases. Interestingly, the elliptical shape and 1.7-axial ratio for the major and minor axes of the fibroin brin cross-section in bagworm silk, which we previously reported14, are already determined at this stage in the silk press; afterwards, the diameter decreases without any change in the axial ratio of the elliptical cross-section. Notably, the two elliptical fibroin brins are aligned side-by-side so that their major axes are in line horizontally (to resemble a figure of ‘∞’) at the spinning press, and these are followed by the spinning tube (Fig. 4e–h). However, the alignment is twisted by 90° in one direction (to resemble a figure of ‘8’) before the brins are spun from the spigot (Fig. 4i,j).We found that the spinning tube was surrounded by a hard exoskeleton. Using 3D-X-ray CT observations, we produced clear images of the exterior and interior morphologies of the spinning tube enveloped by exoskeleton (Fig. 4k–m; the exterior shape observed from the dorsal-, ventral- and lateral-sides by optical microscopy is provided in Supplementary Fig. S3). The spigot was not cut perpendicularly to the spinning tube but rather with a slope of around 20°; consequently, it was elliptic. X-ray CT clearly showed the core-sheath structure of the spinneret and a wide expanse of sheath parts (Fig. 4m) between the exterior shell and interior spinning tube (Fig. 4l,m). Using optical microscope observations of the cross-sections, we found that at least three pairs of adhesive ducts were running in the sheath space (Supplementary Fig. S2E). Therefore, while the silk brins pass through the central narrow spinning tube, the plural adhesive ducts pass through the outer space independently of the silk thread. Finally, the adhesive enters a ladle-like reservoir located at the spigot and is released together with the silk thread (Fig. 4o). The presence of definitive routes connecting the adhesive passage and the spigot were not clearly observed in our X-ray CT images, probably due to the small structural scale relative to the space resolution used in our analysis (i.e. 0.31 μm). We speculate that the adhesive merges into the spigot via a fine, porous sponge-like structure, and we indicate assumed routes in Fig. 4l,m. X-ray CT observations also revealed a sophisticated structural design involving gradual twists in the silk brins by 90° from ‘∞’ to ‘8’ (Fig. 4n and Supplementary Movie S4). Essentially identical spinneret structures were observed by X-ray CT images for all of eight observed individuals from the third to final instars of E. variegata. More