Ecology

These authors contributed equally: Marc Palahí, Ruben ValbuenaEuropean Forest Institute, Joensuu, FinlandMarc Palahí, Lauri Hetemäki, Pieter Johannes Verkerk & Minna KorhonenSchool of Natural Sciences, Bangor University, Bangor, UKRubén ValbuenaEcosystem Dynamics and Forest Management Group, Technical University of Munich, Munich, GermanyCornelius Senf & Rupert SeidlSchool of Geography, Earth and Environmental Sciences, University of Birmingham, Birmingham, UKNezha Acil, Thomas A. M. Pugh & Jonathan SadlerBirmingham Institute of Forest Research, University of Birmingham, Birmingham, UKNezha Acil, Thomas A. M. Pugh & Jonathan SadlerDepartment of Physical Geography and Ecosystem Science, Lund University, Lund, SwedenThomas A. M. PughDepartment of Geographical Sciences, University of Maryland, College Park, MD, USAPeter PotapovInstitut Européen de la Forêt Cultivée, Cestas, FranceBarry GardinerDipartimento di Scienze e Tecnologie Agrarie, Alimentari, Ambientali e Forestali, Università degli Studi di Firenze, Florence, ItalyGherardo Chirici & Saverio FranciniDipartimento per l’Innovazione dei Sistemi Biologici, Agroalimentari e Forestali, Università degli Studi della Tuscia, Viterbo, ItalySaverio FranciniFaculty of Forestry and Wood Sciences, Czech University of Life Sciences Prague, Prague, Czech RepublicTomáš Hlásny & Róbert MarušákWageningen Environmental Research, Wageningen University and Research, Wageningen, The NetherlandsBas Jan Willem Lerink & Gert-Jan NabuursDepartment of Forest Resource Management, Swedish University of Agricultural Sciences (SLU), Umeå, SwedenHåkan Olsson & Jonas FridmanJoint Research Unit CTFC – AGROTECNIO, Solsona, SpainJosé Ramón González OlabarriaDepartment of Agricultural, Forest and Food Sciences, University of Turin, Grugliasco, ItalyDavide AscoliNatural Resources Institute Finland, Joensuu, FinlandAntti AsikainenAlbert-Ludwigs-University of Freiburg, Freiburg, GermanyJürgen Bauhus & Marc HanewinkelDepartment of Space, Earth and Environment, Chalmers University of Technology, Gothenburg, SwedenGöran BerndesLatvian State Forest Research Institute Silava, Salaspils, LatviaJanis DonisINRAE, University of Bordeaux, BIOGECO, Cestas, FranceHervé JactelEuropean Forest Institute, Bonn, GermanyMarcus LindnerUniversity of Molise, Campobasso, ItalyMarco MarchettiForest Ecology and Forest Management Group, Wageningen University and Research, Wageningen, The NetherlandsDouglas Sheil & Gert-Jan NabuursCentro de Estudos Florestais, Instituto Superior de Agronomia, Universidade de Lisboa, Lisbon, PortugalMargarida ToméForest Science and Technology Centre of Catalonia, CTFC, Solsona, SpainAntoni TrasobaresM.P. and G.-J.N. conceived and initiated the study. R.V., C.S., N.A., T.A.M.P., G.C., S.F., T.H., B.J.W.L. and D.A. ran different parts of the analyses and demonstrations. M.P., R.V., G.-J.N., C.S., T.A.M.P., J.S., R.S., B.G. and L.H. drafted the initial version of the manuscript. P.P. provided first-hand experience of the algorithms involved in the production of GFC data. All authors offered insights from their own national statistics and local knowledge, which focused the analyses and the argumentation, and contributed critically to the interpretation of the results, revising and approving the final version of the manuscript. More

Identification of ultra-small cells associated with blooms of anoxygenic phototrophic gammaproteobacteriaThe Salada de Chiprana (NE Spain) is the only permanent athalassic hypersaline lake in Western Europe. It harbors thick, conspicuous microbial mats covering its bottom (Fig. 1a, b) and exhibits periodic stratification during which the deepest part of the water column becomes anoxic and sulfide-rich, favoring the massive development of sulfide-dependent anoxygenic photosynthetic bacteria16. We collected microbial mat fragments that were maintained in culture in the laboratory. After several weeks, we observed a bloom of anoxygenic photosynthetic bacteria containing numerous intracellular sulfur granules (Fig. 1c). Many of these bacterial cells showed one or several much smaller and darker non-flagellated cells attached to their surface (Fig. 1d, e). The infected photosynthetic cells were highly mobile, swimming at high speed with frequent changes of direction (Supplementary Movie 1), in contrast with the non-infected cells that displayed slower (approximately half speed) and more straight swimming. Sometimes, two or more photosynthetic cells were connected through relatively short filaments formed by stacked epibiont cells (Fig. 1f). Although the photosynthetic cells carrying these epibionts were often actively swimming, in some cases the epibionts were associated to empty ghost cells where only the photosynthesis-derived sulfur granules persisted. Closer scrutiny of the epibionts revealed that they actually consisted of piles of up to 10 very flattened cells of 550 ± 50 nm diameter and 220 ± 20 nm height (n = 100). These characteristics (size, morphology, and specific attachment to sulfide-dependent anoxygenic photosynthetic bacteria) perfectly matched the morphological description of the genus Vampirococcus observed over forty years ago in sulfidic freshwater lakes15.Fig. 1: Sampling site and microscopy observation of Vampirococcus cells.a General view of the microbial mat covering the shore of the Salada de Chiprana lake. b Closer view of a microbial mat section. c Natural population of blooming sulfide-dependent anoxygenic photosynthetic bacteria in waters of microbial mat containers after several weeks of growth in the laboratory; note the conspicuous refringent intracellular sulfur inclusions. d–f Closer microscopy view of anoxygenic photosynthetic bacteria infected by epibiotic Vampirococcus cells and few-cell filaments (indicated by yellow arrows). g Scanning electron microscopy image of a host cell infected by two stacking Vampirococcus cells (yellow arrow). h Transmission electron microscopy (TEM) image of a thin section of a host cell infected by Vampirococcus (yellow arrow). i Closer TEM view of a thin section of Vampirococcus cells, notice the fibrous rugose cell surface and the large space separating contiguous cells. Scale bars: 5 cm (b), 5 µm (c), 1 µm (d–h), 0.5 µm (i).Full size imageSince the first Vampirococcus description included transmission electron microscopy (TEM) images, to further ascertain this identification we examined our Chiprana Lake samples under TEM and scanning electron microscopy (SEM). SEM images confirmed the peculiar structure of the epibionts, with multiple contiguous cells separated by deep grooves (Fig. 1g). Thin sections observed under TEM confirmed that the cells were actually separated by a space of ~20–50 nm filled by a fibrous material (Fig. 1h, i). The space between epibiont and host cells was larger (~100 nm) and also filled by dense fibrous material (Fig. 1h). The sections also showed that, in contrast with the typical Gram-negative double membrane structure of the host, the epibiont cells had a single membrane surrounded by a thick layer of fibrous material that conferred a rugose aspect to the cells (Fig. 1i). In sharp contrast with the often highly vacuolated cytoplasm of the host, the epibiont cells showed a dense, homogeneous content. These observations were also in agreement with those published for Vampirococcus, reinforcing our conclusion that the epibionts we observed belonged to this genus, although most likely to a different species, as the first described Vampirococcus occurred in a non-hypersaline lake15.Using a micromanipulator coupled to an inverted microscope, we collected cells of the anoxygenic photosynthetic bacterium carrying Vampirococcus attached to their surface (Supplementary Fig. 1) and proceeded to amplify, clone, and sequence their 16 S rRNA genes. We were able to obtain sequences for both the epibiont and the host for ten infected cells and, in all cases, we retrieved the same two sequences. The host was found to be a Halochromatium-like gammaproteobacterium (Supplementary Fig. 2). Phylogenetic analysis of the epibiont sequence showed that it branched within the CPR radiation close to the Absconditabacteria (Supplementary Fig. 3), previously known as candidate phylum SR12. Since all host and epibiont cells we analyzed had identical 16 S rRNA gene sequences, suggesting that they were the result of a clonal bloom, we collected three sets of ca. 10 infected cells and carried out whole genome amplification (WGA) before sequencing (Illumina HiSeq; see Methods). This strategy allowed us to assemble the nearly complete genome sequence of the Vampirococcus epibiont (see below). In contrast with the completeness of this genome, we only obtained a very partial assembly (~15%) of the host genome, probably because of the consumption of the host DNA by the epibiont. To make more robust phylogenetic analyses of Vampirococcus, we retrieved the protein sequence set used by Hug et al. to reconstruct a multi-marker large-scale phylogeny of bacteria2. The new multi-gene maximum likelihood (ML) phylogenetic tree confirmed the affiliation of our Vampirococcus species to the Absconditabacteria with maximum support, and further placed this clade within a larger well-supported group also containing the candidate phyla Gracilibacteria and Peregrinibacteria (Fig. 2a and Supplementary Fig. 4). Therefore, our epibiotic bacterium represents the first characterized member of this large CPR clade and provides a phylogenetic identity for the predatory bacterial genus Vampirococcus described several decades ago. We propose to call this new species Candidatus Vampirococcus lugosii (see Taxonomic appendix).Fig. 2: Phylogeny and global gene content of the Vampirococcus genome.a Maximum likelihood phylogenetic tree of bacteria based on a concatenated dataset of 16 ribosomal proteins showing the position of Vampirococcus lugosii close to the Absconditabacteria (for the complete tree, see Supplementary Fig. 4). Histograms on the right show the proportion of genes retained in each species from the ancestral pool inferred for the last common ancestor of Absconditabacteria, Gracilibacteria and Peregrinibacteria. b Percentage of Vampirococcus genes belonging to the different Clusters of Orthologous Groups (COG) categories. c Genes shared by Vampirococcus and the three Absconditabacteria genomes shown in the phylogenetic tree. COG categories are: Energy production and conversion [C]; Cell cycle control, cell division, chromosome partitioning [D]; Amino acid transport and metabolism [E]; Nucleotide transport and metabolism [F]; Carbohydrate transport and metabolism [G]; Coenzyme transport and metabolism [H]; Lipid transport and metabolism [I]; Translation, ribosomal structure and biogenesis [J]; Transcription [K]; Replication, recombination and repair [L]; Cell wall/membrane/envelope biogenesis [M]; Secretion, motility and chemotaxis [N]; Posttranslational modification, protein turnover, chaperones [O]; Inorganic ion transport and metabolism [P]; General function prediction only [R]; Function unknown [S]; Intracellular trafficking, secretion, and vesicular transport [U]; Defense mechanisms [V]; Mobilome: prophages, transposons [X]; Secondary metabolites biosynthesis, transport and catabolism [Q]. Source data are provided as a Source Data file.Full size imageGenomic evidence of adaptation to predatory lifestyleWe sequenced DNA from three WGA experiments corresponding each to ~10 Halochromatium-Vampirococcus consortia. Many of the resulting (57.2 Mb) raw sequences exhibited similarity to those of available Absconditabacteria/SR1 metagenome-assembled genomes (MAGs) and, as expected, some also to Gammaproteobacteria (host-derived sequences) as well as a small proportion of potential contaminants probably present in the original sample (Bacillus- and fungi-like sequences). To bin the Vampirococcus sequences out of this mini-metagenome, we applied tetranucleotide frequency analysis on the whole sequence dataset using emergent self-organizing maps (ESOM)6. One of the ESOM sequence bins was enriched in Absconditabacteria/SR1-like sequences and corresponded to the Vampirococcus sequences, which we extracted and assembled independently. This approach yielded an assembly of 1,310,663 bp. We evaluated its completeness by searching i) a dataset of 40 universally distributed single-copy genes17 and ii) a dataset of 43 single-copy genes widespread in CPR bacteria8. We found all them as single-copy genes in the Vampirococcus genome, with the exception of two signal recognition particle subunits from the first dataset which are absent in many other CPR bacteria18. These results supported that the Vampirococcus genome assembly was complete and did not contain multiple strains or other sources of contamination. Manually curated annotation predicted 1151 protein-coding genes, a single rRNA gene operon, and 38 tRNA coding genes. As already found in other Absconditabacteria/SR1 genomes19, the genetic code of Vampirococcus is modified, with the stop codon UGA reassigned as an additional glycine codon.A very large proportion of the predicted proteins (48.9%) had no similarity to any COG class and lacked any conserved domain allowing their functional annotation (Fig. 2b). Thus, as for other CPR bacteria, a significant part of their cellular functions remains inaccessible. A comparison with three other Absconditabacteria genomes revealed a very small set of only 390 genes conserved in all them (Fig. 2c), suggesting a highly dynamic evolution of gene content in these species. Comparison with more distantly related CPR groups (Gracilibacteria and Peregrinibacteria) showed that gene loss has been a dominant trend in all these organisms, which have lost 30–50% of the 1124 genes inferred to have existed in their last common ancestor (Fig. 2a). Nevertheless, this loss of ancestral genes was accompanied by the acquisition of new ones by different mechanisms, including horizontal gene transfer (HGT). In the case of Vampirococcus, we detected, by phylogenetic analysis of all individual genes that had homologs in other organisms, the acquisition of 126 genes by HGT from various donors (Supplementary Data 1).The set of genes that could be annotated provided interesting clues about the biology and lifestyle of Vampirococcus. The most striking feature was its oversimplified energy and carbon metabolism map (Fig. 3). ATP generation in this CPR species appeared to depend entirely on substrate-level phosphorylation carried out by the phosphoenolpyruvate kinase (EC 2.7.1.40). In fact, Vampirococcus only possesses incomplete glycolysis, which starts with 3-phosphoglycerate as first substrate. This molecule is the major product of the enzyme RuBisCO and, therefore, most likely highly available to Vampirococcus from its photosynthetic host. Comparison with nearly complete MAG sequences available for other Absconditabacteria/SR1 showed that Vampirococcus has the most specialized carbon metabolism, with 3-phosphoglycerate as the only exploitable substrate, whereas the other species have a few additional enzymes that allow them to use other substrates (such as ribulose-1,5 P and acetyl-CoA) as energy and reducing power (NADH) sources (Supplementary Fig. 5). This metabolic diversification probably reflects their adaptation to other types of hosts where these substrates are abundant. Vampirococcus also lacks all the enzymes involved in some Absconditabacteria/SR1 in the 3-phosphoglycerate-synthesizing AMP salvage pathway20, including the characteristic archaeal-like type II/III RuBisCO21.Fig. 3: Metabolic and cell features inferred from the genes encoded in the Vampirococcus genome.The diagram shows the host cell surface (bottom) with two stacking Vampirococcus cells attached to its surface (as in Fig. 1h).Full size imageThe genomes of Absconditabacteria/SR1 and Vampirococcus encode several electron carrier proteins (e.g., ferredoxin, cytochrome b5, several Fe-S cluster proteins) and a membrane F1FO-type ATP synthase. However, they apparently lack any standard electron transport chain and, therefore, they seem to be non-respiring19,20,22. The electron carrier proteins may be related with the oxidative stress response and/or the reoxidation of reduced ferredoxin or NADH20. In the absence of any obvious mechanism to generate proton motive force (PMF), the presence of the membrane ATP synthase is also intriguing. It has been speculated either that CPR bacteria might tightly adhere to their hosts and scavenge protons from them or that the membrane ATP synthase might work in the opposite direction as an ATPase, consuming ATP generated by substrate level phosphorylation to extrude protons and drive antiporters9. However, in the case of Vampirococcus the direct transport of protons from the host is unlikely since, as observed in the TEM sections (Fig. 1h), it seems that there is no direct contact with the host cell membrane. In fact, parasite and host cell membranes are separated by a relatively large space of ~100 nm, which would be largely conducive to proton diffusion and inefficient transfer between cells. Alternatively, since the Chiprana lake has high Na+ concentration (1.6 g l−1), it might be possible that the ATP synthase uses Na+ instead of protons. However, the Na+-binding domain of the subunit c of typical Na+-dependent ATP synthases exhibited several differences with that of Vampirococcus (Supplementary Fig. 6). Similar differences have been considered indicative of the use of protons instead of Na+ in other organisms23. Although the proton/cation antiporters (e.g., for Na+, K+, or Ca2+) encoded by Vampirococcus and the other Absconditabacteria/SR1 may serve to produce some PMF, it is improbable that this mechanism represents a major energy transducing system as cells would accumulate cations and disrupt their ionic balance; these antiporters are most likely involved in cation homeostasis.These observations prompted us to investigate other ways that these cells might use to generate PMF usable by their ATP synthase. We found a protein (Vamp_33_45) with an atypical tripartite domain structure. The N-terminal region, containing 8 transmembrane helices, showed similarity with several flavocytochromes capable of moving electrons and/or protons across the plasma membrane (e.g., 24). The central part of the protein was a rubredoxin-like nonheme iron-binding domain likely able to transport electrons. Finally, the C-terminal region, containing an NAD-binding motif, was similar to ferredoxin reductases involved in electron transfer25. This unusual Vampirococcus 3-domain protein is well conserved in the other Absconditabacteria/SR1 genomes sequenced so far, suggesting it plays an important function in this CPR phylum. Its architecture suggests that it can transport electrons and/or protons across the membrane using ferredoxin as electron donor and makes it a strong candidate to participate in a putative new PMF-generating system. Alternatively, this protein could play a similar role to that of some oxidoreductases in the strict anaerobic archaeon Thermococcus onnurineus, including a thioredoxin reductase, which couple reactive oxygen species detoxification with NAD(P) + regeneration from NAD(P)H to maintain the intracellular redox balance and enhance O2-mediated growth despite the absence of heme-based or cytochrome-type proteins26.Although our Vampirococcus genome sequence appears to be complete, genes encoding enzymes involved in the biosynthesis of essential cell building blocks such as amino acids, nucleotides and nucleosides, cofactors, vitamins, and lipids are almost completely absent (Fig. 2b). Therefore, the classical bacterial metabolic pathways for their synthesis27 do not operate in Vampirococcus. Such simplified metabolic potential, comparable to that of intracellular parasitic bacteria such as Mycoplasma28, implies that Vampirococcus must acquire these molecules from an external source and supports the predatory nature of the interaction with its photosynthetic host. An intriguing aspect of this interaction concerns the transfer of substrates from the host to Vampirococcus, especially considering that, despite examination of several serial ultrathin sections, the cell membranes of these two partners do not appear to be in direct contact (Fig. 1h). Vampirococcus encodes several virulence factors, including divergent forms of hemolysin and hemolysin translocator (Vamp_11_169 and Vamp_9_166, respectively), a phage holin (Vamp_5_129), and a membrane-bound lytic murein transglycosylase (Vamp_144_2). These proteins are likely involved in the host cell wall and membrane disruption leading to cell content release. Hemolysin has also been found in Saccharibacteria (formerly candidate phylum TM7), the only CPR phylum for which an epibiotic parasitic lifestyle has been demonstrated so far13,14. Recent coupled lipidomic-metagenomic analyses have shown that CPR bacteria that lack complete lipid biosynthesis are able to recycle membrane lipids from other bacteria29. In Vampirococcus, also devoid of phospholipid synthesis, a phospholipase gene (Vamp_34_196) predicted to be secreted and that has homologs involved in host phospholipid degradation in several parasitic bacteria30, may not only help disrupting the host membrane but also to generate a local source of host phospholipids that it can use to build its own cell membrane. Two Vampirococcus peptidoglycan hydrolases (Vamp_68-56_103 and Vamp_145_30), also predicted to be secreted, most probably contribute to degrade the host cell wall. The Vampirococcus genome also encodes two murein DD-endopeptidases (Vamp_311_38 and Vamp_41_33). As in other predatory bacteria, such as Bdellovibrio, one probably acts to degrade the prey cell wall whereas the other is involved in self-wall remodeling31. Despite their high sequence divergence, we could align both Vampirococcus sequences with those of Bdellovibrio and other bacteria (Supplementary Fig. 7). Both sequences conserved the characteristic active site serine residue of DD endopeptidases and, in contrast with the Bdellovibrio “predatory” enzymes, also the regulatory domain III. The deletion of this regulatory domain has been associated with the capacity of the Bdellovibrio “predatory” DD endopeptidases to act promiscuously on a wide variety of peptidoglycan substrates31. This difference most likely reflects that, whereas Bdellovibrio is able to prey on very diverse bacteria, Vampirococcus is a specialized predator of Chromatiaceae that has evolved specialized enzymes to degrade the wall of its particular prey. The Vampirococcus enzymes could also be aligned with the region where the ankyrin-repeat-containing self-protective regulatory inhibitor Bd3460 binds the Bdellovibrio “predatory” DD endopeptidases32, although only partially for the C-terminal part of Vamp_41_33, like in the self-wall Bdellovibrio enzyme Bd3244 (Supplementary Fig. 7). In that sense, the Vamp_311_38 enzyme seems more similar to the “predatory” Bdellovibrio ones. Interestingly, the “predatory” endopeptidase Bd3459 and the regulatory inhibitor Bd3460 are contiguous in the genomes of Bdellovibrio and other periplasmic predators but not in epibiotic predators32. Vampirococcus confirms this pattern since, although it possesses several ankyrin-repeat-containing proteins, none of them is encoded adjacent to the DD endopeptidase genes.Vampirococcus also possesses a number of genes encoding transporters, most of them involved in the transport of inorganic molecules (Fig. 3). One notable exception is the competence-related integral membrane protein ComEC (Vamp_67_106)33 which, together with ComEA (Vamp_21_186) and type IV pili (see below), probably plays a role in the uptake of host DNA that, once transported into the epibiont, can be degraded by various restriction endonucleases (five genes encoding them are present) and recycled to provide the nucleotides necessary for growth (Supplementary Fig. 8). These proteins are widespread in other CPR bacteria where they may have a similar function34. Vampirococcus also encodes an ABC-type oligopeptide transporter (Vamp_40_40) and a DctA-like C4-dicarboxylate transporter (Vamp_41_97), known to catalyze proton-coupled symport of several Krebs cycle dicarboxylates (succinate, fumarate, malate, and oxaloacetate)35. The first, coupled with the numerous peptidases present in Vampirococcus, most likely is a source of amino acids. By contrast, the role of DctA is unclear since Vampirococcus does not have a Krebs cycle.In sharp contrast with its simplified central metabolism, Vampirococcus possesses genes related to the construction of an elaborate cell surface, which seems to be a common theme in many CPR bacteria9,12. They include genes involved in peptidoglycan synthesis, several glycosyltransferases, a Sec secretion system, and a rich repertoire of type IV pilus proteins. The retractable type IV pili are presumably involved in the tight attachment of Vampirococcus to its host and in DNA uptake in cooperation with the ComEC protein. Other proteins probably play a role in the specific recognition and fixation to the host, including several very large proteins. In fact, the Vampirococcus membrane proteome is enriched in giant proteins. The ten longest predicted proteins (between 1392 and 4163 aa, see Supplementary Table 1) are inferred to have a membrane localization and are probably responsible of the conspicuous fibrous aspect of its cell surface (Fig. 1i). Most of these proteins possess domains known to be involved in the interaction with other molecules, including protein-protein (WD40, TRP, and PKD domains) and protein–lipid (saposin domain) interactions and cell adhesion (DUF11, integrin, and fibronectin domains). Two other large membrane proteins (Vamp_6_203, 2368 aa, and Vamp_19_245, 1895 aa) may play a defensive role as they contain alpha-2-macroglobulin protease-inhibiting domains that can protect against proteases released by the host. Several other smaller proteins complete the membrane proteome of Vampirococcus, some of them also likely involved in recognition and attachment to the host thanks to a variety of protein domains, such as VWA (Vamp_41_85) and flotillin (Vamp_11_100). We did not detect genes coding for flagellar components, confirming the absence of flagella observed under the microscope (Fig. 1).New CRISPR-Cas systems and other defense mechanisms in Vampirococcus
Although most CPR phyla are devoid of CRISPR-Cas36, some have been found to contain new systems with original effector enzymes such as CasY37. In contrast with most available Absconditabacteria genomes, Vampirococcus possesses two CRISPR-Cas loci (Fig. 4a and Supplementary Fig. 9). The first is a class II type V system that contains genes coding for Cas1, Cas2, Cas4, and Cpf1 proteins associated to 34 spacer sequences of 26–32 bp. Proteins similar to those of this system are encoded not only in genomes of close relatives of the Absconditabacteria (Gracilibacteria and Peregrinibacteria) but in many other CPR phyla. These sequences form monophyletic groups in phylogenetic analyses (e.g., Cas1, see Fig. 4b), which suggests that this type V system is probably ancestral in these CPR. The second system found in Vampirococcus belongs to the class I type III and contains genes coding for Cas1, Cas2, Csm3, and Cas10/Csm1 proteins associated to a cluster of 20 longer (35–46 bp) spacers. In contrast with the previous CPR-like system, the proteins of this second system did not show strong similarity with any CPR homolog but with sequences from other bacterial phyla, suggesting that they have been acquired by HGT. Phylogenetic analysis confirmed this and supported that Vampirococcus gained this CRISPR-Cas system from different distant bacterial donors (Supplementary Fig. 10). Interestingly, these two CRISPR-Cas systems encode a number of proteins that may represent new effectors. A clear candidate is the large protein Vamp_48_93 (1158 aa), located between Cpf1 and Cas1 in the type V system (Fig. 4a), which contains a DNA polymerase III PolC motif. Very similar sequences can be found in a few other CPR (some Roizmanbacteria, Gracilibacteria, and Portnoybacteria) and in some unrelated bacteria (Supplementary Fig. 11). As in Vampirococcus, the gene coding for this protein is contiguous to genes encoding different Cas proteins in several of these bacteria, including Roizmanbacteria, Omnitrophica, and the deltaproteobacterium Smithella sp. (Supplementary Fig. 11). This gene association, as well as the very distant similarity between this protein and Cpf1 CRISPR-associated proteins of bacterial type V systems, supports that it is a new effector in type V CRISPR-Cas systems. Additional putative new CRISPR-associated proteins likely exist also in the Vampirococcus type III system (Fig. 4a). Three proteins encoded by contiguous genes (Vamp_21_116, Vamp_21_127, and Vamp_21_128) exhibit very distant similarity with type III-A CRISPR-associated Repeat Associated Mysterious Proteins (RAMP) Csm4, Csm5, and Csm6 sequences, respectively, and most probably represent new RAMP subfamilies. To date, Absconditabacteria38 and Saccharibacteria39 are the only CPR phyla for which phages have been identified. Because of its proximity to Absconditabacteria, Vampirococcus is probably infected by similar phages, so that the function of its CRISPR-Cas systems may be related to the protection against these genetic parasites. Nevertheless, we did not find any similarity between the Vampirococcus spacers and known phage sequences, suggesting that it is infected by unknown phages. Alternatively, considering that Vampirococcus -as most likely many other CPR bacteria- seems to obtain nucleotides required for growth by uptaking host DNA, an appealing possibility is that the CRISPR-Cas systems participate in the degradation of the imported host DNA.Fig. 4: CRISPR-Cas systems in Vampirococcus.a Genes in the two systems encoded in the Vampirococcus lugosii genome, elements common to the two systems are highlighted in blue. b Maximum likelihood phylogenetic tree of the Cas1 protein encoded in the class II type V system, numbers at branches indicate bootstrap support.Full size imageAlthough CPR bacteria have been hypothesized to be largely depleted of classical defense mechanisms40, we found that Vampirococcus, in addition to the two CRISP-Cas loci, is endowed with various other protection mechanisms. These include an AbiEii-AbiEi Type IV toxin-antitoxin system, also present in other CPR bacteria, which may offer additional protection against phage infection41 and several restriction-modification systems, with three type I, one type II and one type III restriction enzymes and eight DNA methylases. In addition to a defensive role, these enzymes may also participate in the degradation of the host DNA. As in its sister-groups Absconditabacteria and Gracilibacteria5,19,20,42,43, Vampirococcus has repurposed the UGA stop codon to code for glycine. The primary function of this recoding remains unknown but it has been speculated that it creates a genetic incompatibility, whereby these bacteria would be “evolutionarily isolated” from their environmental neighbors, preventing their potential competitors from acquiring their genomic innovations by HGT19. However, the opposite might be argued as well, since the UGA codon reassignment can protect Vampirococcus from foreign DNA expression upon uptake by leading to aberrant protein synthesis via read-through of the UGA stop with Gly insertion. This can be important for these CPR bacteria because they are not only impacted by phages38 but they most likely depend on host DNA import and degradation to fulfill their nucleotide requirements. In that sense, it is interesting to note that the Vampirococcus ComEA protein likely involved in DNA transport44 is encoded within the class I type III CRISPR-Cas system (Fig. 4a). More

Invertebrate composition effects on primary productionTo evaluate the effects of fiddler crabs, marsh crabs, and mussels on benthic algae and cordgrass production, the dietary sources for fiddler and marsh crabs, respectively27,28, we measured benthic diatom biomass and cordgrass stem density every 4–6 weeks and quantified cordgrass biomass and grazing damage at the conclusion of the experiment in August 2017. Diatom biomass was enhanced in enclosures with mussels and/or marsh crabs relative to enclosures with only fiddler crabs or no invertebrates, and relative to all ambient plots (F36, 200 = 1.5; P = 0.04; Tukey’s HSD, all P  More

1.Hartman, K. et al. Cropping practices manipulate abundance patterns of root and soil microbiome members paving the way to smart farming. Microbiome 6, 14 (2018).Article 

Google Scholar 
2.Berzsenyi, Z., Győrffy, B. & Lap, D. Effect of crop rotation and fertilisation on maize and wheat yields and yield stability in a long-term experiment. Eur. J. Agron. 13, 225–244. https://doi.org/10.1016/S1161-0301(00)00076-9 (2000).Article 

Google Scholar 
3.Körschens, M. The importance of long-term field experiments for soil science and environmental research: a review. Plant Soil Environ. 52, 1–8 (2006).
Google Scholar 
4.Zuber, S. M., Behnke, G., Nafziger, E. & Villamil, M. B. Crop rotation and tillage effects on soil physical and chemical properties in Illinois. Agron. J. 107, 971–978. https://doi.org/10.2134/agronj14.0465 (2015).CAS 
Article 

Google Scholar 
5.Stott, D. E. Recommended Soil Health Indicators and Associated Laboratory Procedures. National Soil Health Specialist, Soil Health Division, U.S. Department of Agriculture (USDA), Natural Resources Conservation Service (NRCS), Washington, D.C. (2019).6.Barrios, E. Soil biota, ecosystem services and land productivity. Ecol. Econ. 64, 269–285. https://doi.org/10.1016/j.ecolecon.2007.03.004 (2007).Article 

Google Scholar 
7.Garbeva, P., van Veen, J. A. & van Elsas, J. D. Microbial diversity in soil: selection of microbial populations by plant and soil type and implications for disease suppressiveness. Annu. Rev. Phytopathol. 42, 243–270. https://doi.org/10.1146/annurev.phyto.42.012604.135455 (2004).CAS 
Article 
PubMed 

Google Scholar 
8.Hatfield, J., Prueger, J. & Kustas, W. Spatial and temporal variation of energy and carbon fluxes in central Iowa. Agron. J. 99, 285–296. https://doi.org/10.2134/agronj2005.0116S (2007).CAS 
Article 

Google Scholar 
9.Lehman, R. M. et al. Understanding and enhancing soil biological health: the solution for reversing soil degradation. Sustainability 7, 988–1027. https://doi.org/10.3390/su7010988 (2015).Article 

Google Scholar 
10.Weil, R. R., Islam, K. R., Stine, M. A., Gruver, J. B. & Samson-Liebig, S. E. Estimating active carbon for soil quality assessment: a simplified method for laboratory and field use. Am. J. Alter. Agric. 18, 3–17. https://doi.org/10.1079/AJAA200228 (2003).Article 

Google Scholar 
11.Idowu, O. et al. Farmer-oriented assessment of soil quality using field, laboratory, and VNIR spectroscopy methods. Int. J. Plant-Soil Relatsh. 307, 243–253. https://doi.org/10.1007/s11104-007-9521-0 (2008).CAS 
Article 

Google Scholar 
12.Gunapala, N. & Scow, K. M. Dynamics of soil microbial biomass and activity in conventional and organic farming systems. Soil Biol. Biochem. 30, 805–816. https://doi.org/10.1016/S0038-0717(97)00162-4 (1998).CAS 
Article 

Google Scholar 
13.Plaza-Bonilla, D., Álvaro-Fuentes, J. & Cantero-Martínez, C. Identifying soil organic carbon fractions sensitive to agricultural management practices. Soil Tillage Res. 139, 19–22. https://doi.org/10.1016/j.still.2014.01.006 (2014).Article 

Google Scholar 
14.Mirsky, S., Lanyon, L. & Needelman, B. Evaluating soil management using particulate and chemically labile soil organic matter fractions. Soil Sci. Soc. Am. J. 72, 180–185. https://doi.org/10.2136/sssaj2005.0279 (2008).ADS 
CAS 
Article 

Google Scholar 
15.Wright, S. F. & Upadhyaya, A. Extraction of an abundant and unusual protein from soil and comparison with hyphal protein of arbuscular mycorrhizal fungi. Soil Sci. 161, 575–586 (1996).ADS 
CAS 
Article 

Google Scholar 
16.Jan, M. T., Roberts, P., Tonheim, S. K. & Jones, D. L. Protein breakdown represents a major bottleneck in nitrogen cycling in grassland soils. Soil Biol. Biochem. 41, 2272–2282. https://doi.org/10.1016/j.soilbio.2009.08.013 (2009).CAS 
Article 

Google Scholar 
17.Nannipieri, P. & Eldor, P. The chemical and functional characterization of soil N and its biotic components. Soil Biol. Biochem. 41, 2357–2369. https://doi.org/10.1016/j.soilbio.2009.07.013 (2009).CAS 
Article 

Google Scholar 
18.Weintraub, M. N. & Schimel, J. P. Seasonal protein dynamics in Alaskan arctic tundra soils. Soil Biol. Biochem. 37, 1469–1475. https://doi.org/10.1016/j.soilbio.2005.01.005 (2005).CAS 
Article 

Google Scholar 
19.Ros, G. H., Temminghoff, E. J. M. & Hoffland, E. Nitrogen mineralization: a review and meta-analysis of the predictive value of soil tests. Eur. J. Soil Sci. 62, 162–173. https://doi.org/10.1111/j.1365-2389.2010.01318.x (2011).CAS 
Article 

Google Scholar 
20.Ros, G. H., Hanegraaf, M. C., Hoffland, E. & van Riemsdijk, W. H. Predicting soil N mineralization: relevance of organic matter fractions and soil properties. Soil Biol. Biochem. 43, 1714–1722. https://doi.org/10.1016/j.soilbio.2011.04.017 (2011).CAS 
Article 

Google Scholar 
21.Chang, E.-H., Chung, R.-S. & Tsai, Y.-H. Effect of different application rates of organic fertilizer on soil enzyme activity and microbial population. Soil Sci. Plant Nutr. 53, 132–140. https://doi.org/10.1111/j.1747-0765.2007.00122.x (2007).ADS 
CAS 
Article 

Google Scholar 
22.Liebig, M., Carpenter-Boggs, L., Johnson, J. M. F., Wright, S. & Barbour, N. Cropping system effects on soil biological characteristics in the Great Plains. Renew. Agric. Food Syst. 21, 36–48. https://doi.org/10.1079/RAF2005124 (2006).Article 

Google Scholar 
23.Wright, S. & Upadhyaya, A. Comparison of N-linked oligosaccharides of glomalin from arbuscular mycorrhizal fungi and soils by capillary electrophoresis. Soil Biol. Biochem. 30, 1853–1857 (1998).CAS 
Article 

Google Scholar 
24.Wright, S. F. & Upadhyaya, A. A survey of soils for aggregate stability and glomalin, a glycoprotein produced by hyphae of arbuscular mycorrhizal fungi. Plant Soil 198, 97–107 (1998).CAS 
Article 

Google Scholar 
25.Lovelock, C. E., Wright, S. F., Clark, D. A. & Ruess, R. W. Soil stocks of glomalin produced by arbuscular mycorrhizal fungi across a tropical rain forest landscape. J. Ecol. 92, 278–287. https://doi.org/10.1111/j.0022-0477.2004.00855.x (2004).CAS 
Article 

Google Scholar 
26.Emran, M., Gispert, M. & Pardini, G. Patterns of soil organic carbon, glomalin and structural stability in abandoned Mediterranean terraced lands. Eur. J. Soil Sci. 63, 637–649. https://doi.org/10.1111/j.1365-2389.2012.01493.x (2012).CAS 
Article 

Google Scholar 
27.Nichols, K. A. & Millar, J. Glomalin and soil aggregation under six management systems in the Northern Great Plains, USA. Open J. Soil Sci. 03(08), 5. https://doi.org/10.4236/ojss.2013.38043 (2013).CAS 
Article 

Google Scholar 
28.Rillig, M., Ramsey, P., Morris, S. & Paul, E. Glomalin, an arbuscular-mycorrhizal fungal soil protein, responds to land-use change. Int. J. Plant-Soil Relatsh. 253, 293–299. https://doi.org/10.1023/A:1024807820579 (2003).CAS 
Article 

Google Scholar 
29.Klose, S. & Tabatabai, M. A. Response of phosphomonoesterases in soils to chloroform fumigation. J. Plant Nutr. Soil Sci. 165, 429–434. https://doi.org/10.1002/1522-2624(200208)165:4%3c429::AID-JPLN429%3e3.0.CO;2-S (2002).CAS 
Article 

Google Scholar 
30.Wang, X.-C. & Lu, Q. Beta-glucosidase activity in paddy soils of the Taihu Lake Region, China. Pedosphere 16, 118–124. https://doi.org/10.1016/S1002-0160(06)60033-7 (2006).CAS 
Article 

Google Scholar 
31.Wilson, D. B. Microbial diversity of cellulose hydrolysis. Curr. Opin. Microbiol. 14, 259–263. https://doi.org/10.1016/j.mib.2011.04.004 (2011).CAS 
Article 
PubMed 

Google Scholar 
32.Shewale, J. G. β-Glucosidase: Its role in cellulase synthesis and hydrolysis of cellulose. Int. J. Biochem. 14, 435–443. https://doi.org/10.1016/0020-711X(82)90109-4 (1982).CAS 
Article 
PubMed 

Google Scholar 
33.Acosta-Martínez, V., Reicher, Z., Bischoff, M. & Turco, R. F. The role of tree leaf mulch and nitrogen fertilizer on turfgrass soil quality. Biol. Fertil. Soils 29, 55–61. https://doi.org/10.1007/s003740050524 (1999).Article 

Google Scholar 
34.Krogh, K. et al. Characterization and kinetic analysis of a thermostable GH3 β-glucosidase from Penicillium brasilianum. Appl. Microbiol. Biotechnol. 86, 143–154. https://doi.org/10.1007/s00253-009-2181-7 (2010).CAS 
Article 
PubMed 

Google Scholar 
35.Chen, M. et al. Isolation and characterization of a β-glucosidase from Penicillium decumbens and improving hydrolysis of corncob residue by using it as cellulase supplementation. Enzyme Microb. Technol. 46, 444–449. https://doi.org/10.1016/j.enzmictec.2010.01.008 (2010).CAS 
Article 
PubMed 

Google Scholar 
36.Günata, Z. & Vallier, M.-J. Production of a highly glucose-tolerant extracellular β-glucosidase by three Aspergillus strains. Biotechnol. Lett. 21, 219–223. https://doi.org/10.1023/A:1005407710806 (1999).Article 

Google Scholar 
37.Riou, C., Salmon, J.-M., Vallier, M.-J., Gunata, Z. & Barre, P. Purification, characterization, and substrate specificity of a novel highly glucose-tolerant beta -glucosidase from Aspergillus oryzae. Appl. Environ. Microbiol. 64, 3607 (1998).CAS 
Article 

Google Scholar 
38.Tsukada, T., Igarashi, K., Yoshida, M. & Samejima, M. Molecular cloning and characterization of two intracellular β-glucosidases belonging to glycoside hydrolase family 1 from the basidiomycete Phanerochaete chrysosporium. Appl. Microbiol. Biotechnol. 73, 807–814. https://doi.org/10.1007/s00253-006-0526-z (2006).CAS 
Article 
PubMed 

Google Scholar 
39.Yang, S., Wang, L., Yan, Q., Jiang, Z. & Li, L. Hydrolysis of soybean isoflavone glycosides by a thermostable β-glucosidase from Paecilomyces thermophila. Food Chem. 115, 1247–1252. https://doi.org/10.1016/j.foodchem.2009.01.038 (2009).CAS 
Article 

Google Scholar 
40.Arévalo Villena, M., Úbeda Iranzo, J. F., Gundllapalli, S. B., Cordero Otero, R. R. & Briones Pérez, A. I. Characterization of an exocellular β-glucosidase from Debaryomyces pseudopolymorphus. Enzyme Microb. Technol. 39, 229–234. https://doi.org/10.1016/j.enzmictec.2005.10.018 (2006).CAS 
Article 

Google Scholar 
41.Amouri, B. & Gargouri, A. Characterization of a novel β-glucosidase from a Stachybotrys strain. Biochem. Eng. J. 32, 191–197. https://doi.org/10.1016/j.bej.2006.09.022 (2006).CAS 
Article 

Google Scholar 
42.Okamoto, K., Sugita, Y., Nishikori, N., Nitta, Y. & Yanase, H. Characterization of two acidic β-glucosidases and ethanol fermentation in the brown rot fungus Fomitopsis palustris. Enzyme Microb. Technol. 48, 359–364. https://doi.org/10.1016/j.enzmictec.2010.12.012 (2011).CAS 
Article 
PubMed 

Google Scholar 
43.Singhania, R. R., Patel, A. K., Sukumaran, R. K., Larroche, C. & Pandey, A. Role and significance of beta-glucosidases in the hydrolysis of cellulose for bioethanol production. Biores. Technol. 127, 500–507. https://doi.org/10.1016/j.biortech.2012.09.012 (2013).CAS 
Article 

Google Scholar 
44.Okamoto, K., Nakano, H., Yatake, T., Kiso, T. & Kitahata, S. Purification and some properties of a β-glucosidase from Flavobacterium johnsonae. Biosci. Biotechnol. Biochem. 64, 333–340. https://doi.org/10.1271/bbb.64.333 (2000).CAS 
Article 
PubMed 

Google Scholar 
45.Spano, G. et al. A β-glucosidase gene isolated from wine Lactobacillus plantarum is regulated by abiotic stresses. J. Appl. Microbiol. 98, 855–861. https://doi.org/10.1111/j.1365-2672.2004.02521.x (2005).CAS 
Article 
PubMed 

Google Scholar 
46.Mendes, R. et al. Deciphering the rhizosphere microbiome for disease-suppressive bacteria. Science 332, 1097–1100. https://doi.org/10.1126/science.1203980 (2011).ADS 
CAS 
Article 
PubMed 

Google Scholar 
47.Crookston, R., Kurle, J., Copeland, P. J., Ford, J. H. & Lueschen, W. E. Rotational cropping sequence affects yield of corn and soybean. Agron. J. 83, 108–113 (1991).Article 

Google Scholar 
48.Meese, B. G., Carter, P. R., Oplinger, E. S. & Pendleton, J. W. Corn/soybean rotation effect as influenced by tillage, nitrogen, and hybrid/cultivar. J. Prod. Agric. 4, 74–80 (1991).Article 

Google Scholar 
49.Kelley, K. W., Long, J. H. & Todd, T. C. Long-term crop rotations affect soybean yield, seed weight, and soil chemical properties. Field Crop Res. 83, 41–50. https://doi.org/10.1016/S0378-4290(03)00055-8 (2003).Article 

Google Scholar 
50.Mourtzinis, S. et al. Corn and soybean yield response to tillage, rotation, and nematicide seed treatment. Crop Sci. 57, 1704–1712. https://doi.org/10.2135/cropsci2016.09.0792 (2017).Article 

Google Scholar 
51.Farmaha, B. S. et al. Rotation impact on on-farm yield and input-use efficiency in high-yield irrigated maize-soybean systems. Agron. J. 108, 2313–2321. https://doi.org/10.2134/agronj2016.01.0046 (2016).Article 

Google Scholar 
52.Crookston, R. K. & Kurle, J. E. Corn residue effect on the yield of corn and soybean grown in rotation. Agron. J. 81, 229–232. https://doi.org/10.2134/agronj1989.00021962008100020018x (1989).Article 

Google Scholar 
53.Whiting, K. R. & Crookston, R. K. Host-specific pathogens do not account for the corn-soybean rotation effect. Crop Sci. 33, 539–543. https://doi.org/10.2135/cropsci1993.0011183X003300030024x (1993).Article 

Google Scholar 
54.Copeland, P. J., Allmaras, R. R., Crookston, R. K. & Nelson, W. W. Corn-soybean rotation effects on soil water depletion. Agron. J. 85, 203–210. https://doi.org/10.2134/agronj1993.00021962008500020008x (1993).Article 

Google Scholar 
55.Li, J. et al. Soil-plant indices help explain legume response to crop rotation in a semiarid environment. Front. Plant Sci. https://doi.org/10.3389/fpls.2018.01488 (2018).Article 
PubMed 
PubMed Central 

Google Scholar 
56.Nickel, S. E., Crookston, R. K. & Russelle, M. P. Root growth and distribution are affected by corn-soybean cropping sequence. Agron. J. 87, 895–902. https://doi.org/10.2134/agronj1995.00021962008700050020x (1995).Article 

Google Scholar 
57.Bennett, A. J., Bending, G. D., Chandler, D., Hilton, S. & Mills, P. Meeting the demand for crop production: the challenge of yield decline in crops grown in short rotations. Biol. Rev. 87, 52–71. https://doi.org/10.1111/j.1469-185X.2011.00184.x (2012).Article 
PubMed 

Google Scholar 
58.Johnson, N., Copeland, P. J., Crookston, R. & Pfleger, F. L. Mycorrhizae: possible explanation for yield decline with continuous corn and soybean. Agron. J. 84, 387–390 (1992).Article 

Google Scholar 
59.Chen, S., Porter, P. M., Reese, C. D. & Stienstra, W. C. Crop sequence effects on soybean cyst nematode and soybean and corn yields this research was supported by Minnesota soybean producers check-off funding through Minnesota research and promotion council and Minnesota agric. exp. stn.. Crop Sci. 41, 1843–1849. https://doi.org/10.2135/cropsci2001.1843 (2001).Article 

Google Scholar 
60.Grabau, Z. J. & Chen, S. Determining the role of plant-parasitic nematodes in the corn-soybean crop rotation yield effect using nematicide application: II. Soybean. Agron. J. 108, 1168–1179. https://doi.org/10.2134/agronj2015.0432 (2016).CAS 
Article 

Google Scholar 
61.Hoss, M., Behnke, G., Davis, A., Nafziger, E. & Villamil, M. B. Short corn rotations do not improve soil quality. Compared with corn monocultures. Agron. J. 110, 1274–1288. https://doi.org/10.2134/agronj2017.11.0633 (2018).CAS 
Article 

Google Scholar 
62.Plaza, C., Courtier-Murias, D., Fernández, J. M., Polo, A. & Simpson, A. J. Physical, chemical, and biochemical mechanisms of soil organic matter stabilization under conservation tillage systems: a central role for microbes and microbial by-products in C sequestration. Soil Biol. Biochem. 57, 124–134. https://doi.org/10.1016/j.soilbio.2012.07.026 (2013).CAS 
Article 

Google Scholar 
63.Tardy, V. et al. Shifts in microbial diversity through land use intensity as drivers of carbon mineralization in soil. Soil Biol. Biochem. 90, 204–213. https://doi.org/10.1016/j.soilbio.2015.08.010 (2015).CAS 
Article 

Google Scholar 
64.Tiedje, J. M., Asuming-Brempong, S., Nüsslein, K., Marsh, T. L. & Flynn, S. J. Opening the black box of soil microbial diversity. Appl. Soil. Ecol. 13, 109–122. https://doi.org/10.1016/S0929-1393(99)00026-8 (1999).Article 

Google Scholar 
65.Hussain, S., Ghaffar, A. & Aslam, M. Biological-control of macrophomina-phaseolina charcoal rot of sunflower and mung bean. J. Phytopathol. 130, 157–160. https://doi.org/10.1111/j.1439-0434.1990.tb01163.x (1990).Article 

Google Scholar 
66.Khan, A. N. et al. Molecular identification and genetic characterization of Macrophomina phaseolina strains causing pathogenicity on sunflower and chickpea. Front. Microbiol. 8, 1309. https://doi.org/10.3389/fmicb.2017.01309 (2017).Article 
PubMed 
PubMed Central 

Google Scholar 
67.Ramezani, M. et al. Soybean charcoal rot disease fungus Macrophomina phaseolina in Mississippi produces the phytotoxin (-)-botryodiplodin but no detectable phaseolinone. J. Nat. Prod. 70, 128–129. https://doi.org/10.1021/np060480t (2007).CAS 
Article 
PubMed 

Google Scholar 
68.Smith, L.J., Datnoff, L.E., Pernezny, K. & Schlub, R.L. Phylogenetic and pathogenic characterization of Corynespora cassiicola isolates. In II International Symposium on Tomato Diseases 808, 51–56 (2007)69.Deon, M. et al. Characterization of a cassiicolin-encoding gene from Corynespora cassiicola, pathogen of rubber tree (Hevea brasiliensis). Plant Sci. 185–186, 227–237. https://doi.org/10.1016/j.plantsci.2011.10.017 (2012).CAS 
Article 
PubMed 

Google Scholar 
70.Videira, S. I. R. et al. Mycosphaerellaceae: chaos or clarity?. Stud. Mycol. 87, 257–421 (2017).CAS 
Article 

Google Scholar 
71.Wijayawardene, N. N. et al. Outline of ascomycota: 2017. Fungal Divers. 88, 167–263. https://doi.org/10.1007/s13225-018-0394-8 (2018).Article 

Google Scholar 
72.Olofsson, J., Ericson, L., Torp, M., Stark, S. & Baxter, R. Carbon balance of Arctic tundra under increased snow cover mediated by a plant pathogen. Nat. Clim. Change 1, 220–223. https://doi.org/10.1038/Nclimate1142 (2011).ADS 
CAS 
Article 

Google Scholar 
73.Wells, L. D. & McManus, P. S. A photographic diagnostic guide for identification of the principal cranberry fruit rot pathogens. Plant Health Prog. https://doi.org/10.1094/php-2013-0729-01-dg (2013).Article 

Google Scholar 
74.Yeager, C. M. et al. Polysaccharide degradation capability of actinomycetales soil isolates from a semiarid grassland of the Colorado Plateau. Appl. Environ. Microbiol. https://doi.org/10.1128/aem.03020-16 (2017).Article 
PubMed 
PubMed Central 

Google Scholar 
75.Loria, R., Bukhalid, R. A., Fry, B. A. & King, R. R. Plant pathogenicity in the genus Streptomyces. Plant Dis. 81, 836–846. https://doi.org/10.1094/Pdis.1997.81.8.836 (1997).Article 
PubMed 

Google Scholar 
76.Li, Y., Liu, J., Diaz-Cruz, G., Cheng, Z. & Bignell, D. R. D. Virulence mechanisms of plant-pathogenic Streptomyces species: an updated review. Microbiology 165, 1025–1040. https://doi.org/10.1099/mic.0.000818 (2019).CAS 
Article 
PubMed 

Google Scholar 
77.Abdalla, M. H. Solubilization of rock phosphates by rhizobium and bradyrhizobium. Folia Microbiol. 39, 53–56. https://doi.org/10.1007/Bf02814530 (1994).CAS 
Article 

Google Scholar 
78.Bargaz, A., Lyamlouli, K., Chtouki, M., Zeroual, Y. & Dhiba, D. Soil microbial resources for improving fertilizers efficiency in an integrated plant nutrient management system. Front Microbiol. https://doi.org/10.3389/fmicb.2018.01606 (2018).Article 
PubMed 
PubMed Central 

Google Scholar 
79.Moebius-Clune, B. N. Comprehensive Assessment of Soil Health: The Cornell Framework Manual (Cornell University, 2016).
Google Scholar 
80.Deng, S. P. & Tabatabai, M. A. Cellulase activity of soils. Soil Biol. Biochem. 26, 1347–1354. https://doi.org/10.1016/0038-0717(94)90216-X (1994).CAS 
Article 

Google Scholar 
81.Riesenfeld, C. S., Goodman, R. M. & Handelsman, J. Uncultured soil bacteria are a reservoir of new antibiotic resistance genes. Environ. Microbiol. 6, 981–989 (2004).CAS 
Article 

Google Scholar 
82.Gohl, D. et al. Systematic improvement of amplicon marker gene methods for increased accuracy in microbiome studies. Nat. Biotechnol. 34, 942–949. https://doi.org/10.1038/nbt.3601 (2016).CAS 
Article 
PubMed 

Google Scholar 
83.R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing (Vienna, Austria, 2020).84.Bates, D., Machler, M., Bolker, B. M. & Walker, S. C. Fitting linear mixed-effects models using lme4. J. Stat. Softw. 67, 1–48 (2015).Article 

Google Scholar 
85.Bolyen, E. et al. Qiime 2: reproducible, interactive, scalable, and extensible microbiome data science. Report No. 2167–9843, (PeerJ Preprints, 2018).86.Callahan, B. J. et al. DADA2: high-resolution sample inference from Illumina amplicon data. Nat. Methods 13, 581–583. https://doi.org/10.1038/nmeth.3869 (2016).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
87.Mandal, S. et al. Analysis of composition of microbiomes: a novel method for studying microbial composition. Microb. Ecol. Health Dis. https://doi.org/10.3402/mehd.v26.27663 (2015).Article 
PubMed 
PubMed Central 

Google Scholar 
88.Morton, J. T. et al. Balance trees reveal microbial niche differentiation. mSystems https://doi.org/10.1128/mSystems.00162-16 (2017).Article 
PubMed 
PubMed Central 

Google Scholar 
89.Nilsson, R. H. et al. The UNITE database for molecular identification of fungi: handling dark taxa and parallel taxonomic classifications. Nucleic Acids Res. 47, D259–D264. https://doi.org/10.1093/nar/gky1022 (2019).CAS 
Article 
PubMed 

Google Scholar 
90.Quast, C. et al. The SILVA ribosomal RNA gene database project: improved data processing and web-based tools. Nucleic Acids Res. 41, D590-596. https://doi.org/10.1093/nar/gks1219 (2013).CAS 
Article 
PubMed 

Google Scholar 
91.Pedregosa, F. et al. Scikit-learn: machine learning in python. J. Mach. Learn. Res. 12, 2825–2830 (2011).MathSciNet 
MATH 

Google Scholar 
92.Kruskal, W. H. & Wallis, W. A. Use of ranks in one-criterion variance analysis. J. Am. Stat. Assoc. 47, 583–621. https://doi.org/10.1080/01621459.1952.10483441 (1952).Article 
MATH 

Google Scholar  More

1.Hewitt, G. M. Quaternary phylogeography: the roots of hybrid zones. Genetica 139, 617–638 (2011).PubMed 
Article 

Google Scholar 
2.Gillespie, R. G. & Roderick, G. K. Evolution: geology and climate drive diversification. Nature 509, 297–298 (2014).ADS 
CAS 
PubMed 
Article 

Google Scholar 
3.Hewitt, G. M. The genetic legacy of the quaternary ice ages. Nature 405, 907–913 (2000).ADS 
CAS 
PubMed 
Article 

Google Scholar 
4.Gómez, A. & Lunt, D. H. Refugia within refugia: patterns of phylogeographic concordance in the Iberian Peninsula. In Phylogeography of Southern European Refugia (eds Weiss, S. & Ferrand, N.) 155–188 (Springer, 2007).
Google Scholar 
5.Abellán, P. & Svenning, J. C. Refugia within refugia—patterns in endemism and genetic divergence are linked to Late Quaternary climate stability in the Iberian Peninsula. Biol. J. Linn. Soc. 113, 13–28 (2014).Article 

Google Scholar 
6.Avise, J. C. Phylogeography: The History and Formation of Species (Harvard University Press, 2000).
Google Scholar 
7.Juan, C., Emerson, B. C., Oromí, P. & Hewitt, G. M. Colonization and diversification: towards a phylogeographic synthesis for the Canary Islands. Trends Ecol. Evol. 15, 104–109 (2000).CAS 
PubMed 
Article 

Google Scholar 
8.Chung, M. Y., López-Pujol, J. & Chung, M. G. The role of the Baekdudaegan (Korean Peninsula) as a major glacial refugium for plant species: a priority for conservation. Biol. Conserv. 206, 236–248 (2017).Article 

Google Scholar 
9.Chung, M. Y. et al. The Korean baekdudaegan mountains: a glacial refugium and a biodiversity hotspot that needs to be conserved. Front. Genet. 9, 489 (2018).PubMed 
PubMed Central 
Article 

Google Scholar 
10.AmphibiaWeb. https://amphibiaweb.org/. (Accessed: 15th July 2020).11.Borzée, A. & Min, M.-S. Disentangling the impacts of speciation, sympatry and the island effect on the morphology of seven Hynobius sp. salamanders. Animals 11, 187 (2021).PubMed 
Article 

Google Scholar 
12.Baek, H.-J., Lee, M.-Y., Lee, H. & Min, M.-S. Mitochondrial DNA data unveil highly divergent populations within the Genus Hynobius (Caudata: Hynobiidae) in South Korea. Mol. Cells 31, 105–112 (2011).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
13.Suk, H. Y. et al. Phylogenetic structure and ancestry of Korean clawed salamander, Onychodactylus koreanus (Caudata: Hynobiidae). Mitochondrial DNA Part A DNA Mapp. Seq. Anal. 29, 650–658 (2018).CAS 
Article 

Google Scholar 
14.Min, M.-S. et al. Discovery of the first Asian plethodontid salamander. Nature 435, 87–90 (2005).ADS 
CAS 
PubMed 
Article 

Google Scholar 
15.Wiens, J. J., Engstrom, T. N. & Chippindale, P. T. Rapid diversification, incomplete isolation, and the “speciation clock” in North American salamanders (genus Plethodon): testing the hybrid swarm hypothesis of rapid radiation. Evolution 60, 2585–2603 (2006).CAS 
PubMed 

Google Scholar 
16.Kozak, K. H., Mendyk, R. W. & Wiens, J. J. Can parallel diversification occur in sympatry? Repeated patterns of body-size evolution in coexisting clades of North American salamanders. Evolution 63, 1769–1784 (2009).PubMed 
Article 

Google Scholar 
17.Zhang, P. & Wake, D. B. Higher-level salamander relationships and divergence dates inferred from complete mitochondrial genomes. Mol. Phylogenet. Evol. 53, 492–508 (2009).CAS 
PubMed 
Article 

Google Scholar 
18.Vieites, D. R., Román, S. N., Wake, M. H. & Wake, D. B. A multigenic perspective on phylogenetic relationships in the largest family of salamanders, the Plethodontidae. Mol. Phylogenet. Evol. 59, 623–635 (2011).PubMed 
Article 

Google Scholar 
19.Shen, X. X. et al. Enlarged multilocus data set provides surprisingly younger time of origin for the Plethodontidae, the largest family of salamanders. Syst. Biol. 65, 66–81 (2016).PubMed 
Article 

Google Scholar 
20.Wake, D. B. Persistent plethodontid themes: species, phylogenies, and biogeography. Herpetologica 73, 242–251 (2017).Article 

Google Scholar 
21.Wake, D. B. The enigmatic history of the European, Asian and American plethodontid salamanders. Amphib-reptil 34, 323–336 (2013).Article 

Google Scholar 
22.Vieites, D. R., Min, M.-S. & Wake, D. B. Rapid diversification and dispersal during periods of global warming by plethodontid salamanders. Proc. Natl. Acad. Sci. U.S.A. 104, 19903–19907 (2007).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
23.IUCN SSC Amphibian Specialist Group. Karsenia koreana. The IUCN Red List of Threatened Species 2019. (2019). https://www.iucnredlist.org/species/61903/110101886. (Accessed: 15th July 2020).24.Borzée, A. et al. Climate change-based models predict range shifts in the distribution of the only Asian plethodontid salamander: Karsenia koreana. Sci. Rep. 9, 1–9 (2019).Article 
CAS 

Google Scholar 
25.Sessions, S. K. et al. Cytogenetic analysis of the Asian Plethodontid salamander, Karsenia koreana: evidence for karyotypic conservation, chromosome repatterning, and genome size evolution. Chromosom. Res. 16, 563–574 (2008).CAS 
Article 

Google Scholar 
26.Buckley, D., Wake, M. H. & Wake, D. B. Comparative skull osteology of Karsenia koreana (Amphibia, Caudata, Plethodontidae). J. Morphol. 271, 533–558 (2010).PubMed 
Article 

Google Scholar 
27.Sever, D. M., Pinsoneault, A. D., Mackenzie, B. W., Siegel, D. S. & Staub, N. L. A description of the skin glands and cloacal morphology of the plethodontid salamander Karsenia koreana. Copeia 104, 816–823 (2016).Article 

Google Scholar 
28.Moon, K. Y. & Park, D. Report of Karsenia koreana eggs oviposited within a semi-natural terrarium constructed at natural habitat. Korean J. Herpetol. 7, 1–5 (2016).
Google Scholar 
29.Song, J.-Y. et al. Life history of a unique Asian plethodontid salamander, Karsenia koreana. Zool. Sci. 34, 122–128 (2017).Article 

Google Scholar 
30.Jung, J.-H., Lee, E.-J., Lee, W.-S. & Park, C.-D. Habitat suitability models of Korean crevice salamander (Karsenia koreana) at forested area in Daejeon metropolitan city, Republic of Korea. J. For. Res. 24, 349–355 (2019).CAS 
Article 

Google Scholar 
31.Rozen, S. & Skaletsky, H. Primer3 on the WWW for general users and for biologist programmers. In Bioinformatics Methods and Protocols: Methods in Molecular Biology Vol. 132 (eds Misener, S. & Krawetz, S. A.) 365–386 (Humana Press, 2000).
Google Scholar 
32.Su, X. Z., Wu, Y., Sifri, C. D. & Wellems, T. E. Reduced extension temperatures required for PCR amplification of extremely A+T-rich DNA. Nucleic Acids Res. 24, 1574–1575 (1996).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
33.Smit, A. F. A., Hubley, R. & Green, P. RepeatModeler Open-1.0. 2008–2015. (2014). http://www.repeatmasker.org.34.Stieneke, D. L. & Eujayl, I. L. Imperfect SSR Finder Version 1.0. United States Department of Agriculture. (2019). https://data.nal.usda.gov/dataset/imperfect-ssr-finder.35.Thompson, J. D., Higgins, D. G. & Gibson, T. J. CLUSTAL W: improving the sensitivity of progressive multiple sequence alignment through sequence weighting, position-specific gap penalties and weight matrix choice. Nucleic Acids Res. 22, 4673–4680 (1994).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
36.Kumar, S., Stecher, G., Li, M., Knyaz, C. & Tamura, K. MEGA X: molecular evolutionary genetics analysis across computing platforms. Mol. Biol. Evol. 35, 1547–1549 (2018).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
37.Librado, P. & Rozas, J. DnaSP v5: a software for comprehensive analysis of DNA polymorphism data. Bioinformatics 25, 1451–1452 (2009).CAS 
Article 

Google Scholar 
38.Nei, M. Molecular Evolutionary Genetics (Columbia University Press, 1987).Book 

Google Scholar 
39.Tajima, F. Evolutionary relationship of DNA sequences in finite populations. Genetics 105, 437–460 (1983).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
40.Excoffier, L., Smouse, P. & Quattro, J. Analysis of molecular variance inferred from metric distances among DNA haplotypes: application to human mitochondrial DNA restriction data. Genetics 131, 479–491 (1992).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
41.Excoffier, L. & Lischer, H. E. L. Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and Windows. Mol. Ecol. Resour. 10, 564–567 (2010).PubMed 
Article 

Google Scholar 
42.Bandelt, H.-J., Forster, P. & Röhl, A. Median-joining networks for inferring intraspecific phylogenies. Mol. Biol. Evol. 16, 37–48 (1999).CAS 
PubMed 
Article 

Google Scholar 
43.Wake, D. B. Taxonomy of salamanders of the family Plethodontidae (Amphibia: Caudata). Zootaxa 3484, 75–82 (2012).Article 

Google Scholar 
44.Lanfear, R., Calcott, B., Ho, S. Y. & Guidon, S. PartitionFinder: combined selection of partitioning schemes and substitution models for phylogenetic analyses. Mol. Biol. Evol. 29, 1695–1701 (2012).CAS 
PubMed 
Article 

Google Scholar 
45.Lanfear, R., Frandsen, P. B., Wright, A. M., Senfeld, T. & Calcott, B. PartitionFinder 2: new methods for selecting partitioned models of evolution for molecular and morphological phylogenetic analyses. Mol. Biol. Evol. 34, 772–773 (2017).CAS 
PubMed 

Google Scholar 
46.Stamatakis, A. RAxML version 8: a tool for phylogenetic analysis and post-analysis of large phylogenies. Bioinformatics 30, 1312–1313 (2014).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
47.Miller, M. A., Pfeiffer, W. & Schwartz, T. Creating the CIPRES Science Gateway for inference of large phylogenetic trees. In Proceedings of the Gateway Computing Environment Workshop (GCE) 1–8 (2010).48.Ronquist, F. & Huelsenbeck, J. P. Bayesian phylogenetic inference under mixed models. Bioinformatics 19, 1572–1574 (2003).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
49.Tavaré, S. Some probabilistic and statistical problems in the analysis of DNA sequences. In Some Mathematical Questions in Biology: DNA Sequence Analysis (ed. Miura, R. M.) 57–86 (American Mathematical Society, 1986).
Google Scholar 
50.Rambaut, A., Drummond, A. J., Xie, D., Baele, G. & Suchard, M. A. Posterior summarization in Bayesian phylogenetics using Tracer 1.7. Syst. Biol. 67, 901–904 (2018).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
51.Bouckaert, R. et al. BEAST 2.5: an advanced software platform for Bayesian evolutionary analysis. PLoS Comput. Biol. 15, e1006650 (2019).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
52.Brown, R. P. & Yang, Z. Rate variation and estimation of divergence times using strict and relaxed clocks. BMC Evol. Biol. 11, 271 (2011).PubMed 
PubMed Central 
Article 

Google Scholar 
53.Skilling, J. Nested sampling for general Bayesian computation. Bayesian Anal. 1, 833–860 (2006).MathSciNet 
MATH 
Article 

Google Scholar 
54.Russel, P. M., Brewer, B. J., Klaere, S. & Bouckaert, R. Model selection and parameter inference in phylogenetics using nested sampling. Syst. Biol. 68, 219–233 (2018).Article 

Google Scholar 
55.Heled, J. & Drummond, A. J. Bayesian inference of population size history from multiple loci. BMC Evol. Biol. 8, 289 (2008).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
56.Hasegawa, M., Kishino, H. & Yano, T. Dating of the human-ape splitting by a molecular clock of mitochondrial DNA. J. Mol. Evol. 22, 160–174 (1985).CAS 
PubMed 
Article 

Google Scholar 
57.R Core Team. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. (2013). http://www.r-project.org/.58.Tajima, F. The effect of change in population size on DNA polymorphism. Genetics 123, 597–601 (1989).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
59.Fu, Y.-X. Statistical tests of neutrality of mutations against population growth, hitchhiking and background selection. Genetics 147, 915–9125 (1997).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
60.van Oosterhout, C., Hutchinson, B., Wills, D. & Shipley, P. MICRO-CHECKER: software for identifying and correcting genotyping errors in microsatellite data. Mol. Ecol. Notes 4, 535–538 (2004).Article 
CAS 

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

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

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

Google Scholar 
64.Garza, J. C. & Williamson, E. G. Detection of reduction in population size using data from microsatellite loci. Mol. Ecol. 10, 305–318 (2001).CAS 
PubMed 
Article 

Google Scholar 
65.Harley, E. H. AGARst: a programme for calculating allele frequencies, Gst and Rst from microsatellite data, version 2. University of Cape Town, Cape Town, South Africa. (2001).66.Wright, S. The interpretation of population structure by F-statistics with special regard to systems of mating. Evolution 19, 395–420 (1965).Article 

Google Scholar 
67.Slatkin, M. A measure of population subdivision based on microsatellite allele frequencies. Genetics 139, 457–462 (1995).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
68.Goudet, J. Fstat (ver. 2.9.4), a program to estimate and test population genetics parameters. (2003). https://www2.unil.ch/popgen/softwares/fstat.htm69.Rousset, F. Genetic differentiation and estimation of gene flow from F-statistics under isolation by distance. Genetics 145, 1219–1228 (1997).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
70.Manni, F., Guérard, E. & Heyer, E. Geographic patterns of (genetic, morphologic, linguistic) variation: how barriers can be detected by “Monmonier’s algorithm”. Hum. Biol. 76, 173–190 (2004).71.Monmonier, M. Maximum-difference barriers: an alternative numerical regionalization method. Geogr. Anal. 3, 245–261 (1973).
Google Scholar 
72.Nei, M., Tajima, F. & Tateno, Y. Accuracy of estimated phylogenetic trees from molecular data. J. Mol. Evol. 19, 153–170 (1983).ADS 
CAS 
PubMed 
Article 

Google Scholar 
73.Dieringer, D. & Schlötterer, C. Microsatellite analyser (MSA): a platform independent analysis tool for large microsatellite data sets. Mol. Ecol. Notes 3, 167–169 (2003).CAS 
Article 

Google Scholar 
74.Kitada, S., Nakamichi, R. & Kishino, H. The empirical Bayes estimators of fine-scale population structure in high gene flow species. Mol. Ecol. Resour. 17, 1210–1222 (2017).CAS 
PubMed 
Article 

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

Google Scholar 
76.Earl, D. A. & Vonholdt, B. M. STRUCTURE HARVESTER: a website and program for visualizing STRUCTURE output and implementing the Evanno method. Conserv. Genet. Resour. 4, 359–361 (2012).Article 

Google Scholar 
77.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 (2005).CAS 
PubMed 
Article 

Google Scholar 
78.Cornuet, J.-M. et al. DIYABC v2.0: a software to make approximate Bayesian computation inferences about population history using single nucleotide polymorphism, DNA sequence and microsatellite data. Bioinformatics 30, 1187–1189 (2014).CAS 
PubMed 
Article 

Google Scholar 
79.Penny, D. Relativity for molecular clocks. Nature 436, 183–184 (2005).ADS 
CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
80.Dakin, E. E. & Avise, J. C. Microsatellite null alleles in parentage analysis. Heredity 93, 504–509 (2004).CAS 
PubMed 
Article 

Google Scholar 
81.Oromi, N. et al. Genetic structure of lake and stream populations in a Pyrenean amphibian (Calotriton asper) reveals evolutionary significant units associated with paedomorphosis. J. Zool. Syst. Evol. Res. 57, 418–430 (2019).Article 

Google Scholar 
82.Janes, J. K. et al. The K = 2 conundrum. Mol. Ecol. 26, 3594–3602 (2017).PubMed 
Article 

Google Scholar 
83.Chiari, Y. et al. Phylogeography of Sardinian cave salamanders (genus Hydromantes) is mainly determined by geomorphology. PLoS ONE 7, e32332 (2012).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
84.Cimmaruta, R., Lucente, D. & Nascetti, G. Persistence, isolation and diversification of a naturally fragmented species in local refugia: the case of Hydromantes strinatii. PLoS ONE 10, e0131298 (2015).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
85.Kuchta, S. R., Haughey, M., Wynn, A. H., Jacobs, J. F. & Highton, R. Ancient river systems and phylogeographical structure in the spring salamander, Gyrinophilus porphyriticus. J. Biogeogr. 43, 639–652 (2016).Article 

Google Scholar 
86.Pan, T. et al. Long-term sky islands generate highly divergent lineages of a narrowly distributed stream salamander (Pachyhynobius shangchengensis) in mid-latitude mountains of East Asia. BMC Evol. Biol. 19, 1–15 (2019).PubMed 
PubMed Central 
Article 

Google Scholar 
87.Suk, H. Y. et al. Genetic and phylogenetic structure of Hynobius quelpaertensis, an endangered endemic salamander species on the Korean Peninsula. Genes Genom. 42, 165–178 (2020).CAS 
Article 

Google Scholar 
88.Cameron, A. C., Anderson, J. J. & Page, R. B. Assessment of intra and interregional genetic variation in the Eastern Red-backed Salamander, Plethodon cinereus, via analysis of novel microsatellite markers. PLoS ONE 12, e0186866 (2017).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
89.Yoshikawa, N. & Nagata, N. Eighteen SSR Markers for the Japanese clawed salamander, Onychodactylus japonicus, and cross-amplification in its congeners. Curr. Herpetol. 36, 153–158 (2017).90.Moritz, C. Applications of mitochondrial DNA analysis in conservation: a critical review. Mol. Ecol. 3, 401–411 (1994).CAS 
Article 

Google Scholar 
91.Estoup, A., Jarne, P. & Cornuet, J. M. Homoplasy and mutation model at microsatellite loci and their consequences for population genetics analysis. Mol. Ecol. 11, 1591–1604 (2002).CAS 
PubMed 
Article 

Google Scholar 
92.Toews, D. P. L. & Brelsford, A. The biogeography of mitochondrial and nuclear discordance in animals. Mol. Ecol. 21, 3907–3930 (2012).CAS 
PubMed 
Article 

Google Scholar 
93.Karl, S. A., Toonen, R. J., Grant, W. S. & Bowen, B. W. Common misconceptions in molecular ecology: echoes of the modern synthesis. Mol. Ecol. 21, 4171–4189 (2012).CAS 
PubMed 
Article 

Google Scholar 
94.Cook, B. D., Bunn, S. E. & Hughes, J. M. Molecular genetic and stable isotope signatures reveal complementary patterns of population connectivity in the regionally vulnerable southern pygmy perch (Nannoperca australis). Biol. Conserv. 138, 60–72 (2007).Article 

Google Scholar 
95.Hewitt, G. M. Some genetic consequences of ice ages, and their role in divergence and speciation. Biol. J. Linn. Soc. 58, 247–276 (1996).Article 

Google Scholar 
96.Avise, J. C., Walker, D. & Johns, G. C. Speciation durations and Pleistocene effects on vertebrate phylogeography. Proc. R. Soc. Lond. Ser. B Biol. Sci. 265, 1707–1712 (1998).CAS 
Article 

Google Scholar 
97.Alexandrino, J., Froufe, E., Arntzen, J. W. & Ferrand, N. Genetic subdivision, glacial refugia and postglacial recolonization in the golden-striped salamander, Chioglossa lusitanica (Amphibia: Urodela). Mol. Ecol. 9, 771–781 (2000).CAS 
PubMed 
Article 

Google Scholar 
98.Alexandrino, J., Arntzen, J. W. & Ferrand, N. Nested clade analysis and the genetic evidence for population expansion in the phylogeography of the golden-striped salamander, Chioglossa lusitanica (Amphibia: Urodela). Heredity 88, 66–74 (2002).CAS 
PubMed 
Article 

Google Scholar 
99.Rovito, S. M. Lineage divergence and speciation in the Web-toed Salamanders (Plethodontidae: Hydromantes) of the Sierra Nevada, California . Mol. Ecol. 19, 4554–4571 (2010).PubMed 
Article 

Google Scholar 
100.Shafer, A. B. A., Cullingham, C. I., Côté, S. D. & Coltman, D. W. Of glaciers and refugia: a decade of study sheds new light on the phylogeography of northwestern North America. Mol. Ecol. 19, 4589–4621 (2010).PubMed 
Article 

Google Scholar 
101.Zhang, R.-Z. Geological events and mammalian distribution in China. Acta Zool. Sin. 48, 141–153 (2002).
Google Scholar 
102.Matsui, M., Tominaga, A., Liu, W. Z. & Tanaka-Ueno, T. Reduced genetic variation in the Japanese giant salamander, Andrias japonicus (Amphibia: Caudata). Mol. Phylogenet. Evol. 49, 318–326 (2008).CAS 
PubMed 
Article 

Google Scholar 
103.Matsui, M. et al. Phylogenetic relationships of two Salamandrella species as revealed by mitochondrial DNA and allozyme variation (Amphibia: Caudata: Hynobiidae). Mol. Phylogenet. Evol. 48, 84–93 (2008).CAS 
PubMed 
Article 

Google Scholar 
104.Malyarchuk, B., Derenko, M. & Denisova, G. Phylogeny and genetic history of the Siberian salamander (Salamandrella keyserlingii, Dybowski, 1870) inferred from complete mitochondrial genomes. Mol. Phylogenet. Evol. 67, 348–357 (2013).CAS 
PubMed 
Article 

Google Scholar 
105.Honda, A. et al. Late Pleistocene climate change and population dynamics of Japanese Myodes voles inferred from mitochondrial cytochrome b sequences. J. Mammal. 100, 1156–1168 (2019).PubMed 
PubMed Central 
Article 

Google Scholar 
106.Moritz, C. Defining “evolutionarily significant units” for conservation. Trends Ecol. Evol. 9, 373–375 (1994).CAS 
PubMed 
Article 

Google Scholar  More

1.Cremer, S., Armitage, S. A. O. & Schmid-Hempel, P. Social immunity. Curr. Biol. 17, R693–R702 (2007).CAS 
Article 

Google Scholar 
2.Page, P. et al. Social apoptosis in honey bee superorganisms. Sci. Rep. 6, 1–6 (2016).Article 

Google Scholar 
3.Winston, M. L. The Biology of the Honey Bee. (Harvard University Press, 1991).4.Beye, M. et al. Exceptionally high levels of recombination across the honey bee genome. Genome Res. 16, 1339–1344 (2006).CAS 
Article 

Google Scholar 
5.Kent, C. F., Minaei, S., Harpur, B. A. & Zayed, A. Recombination is associated with the evolution of genome structure and worker behavior in honey bees. Proc. Natl. Acad. Sci. 109, 18012–18017 (2012).ADS 
CAS 
Article 

Google Scholar 
6.Rosenkranz, P., Aumeier, P. & Ziegelmann, B. Biology and control of Varroa destructor. J. Invertebr. Pathol. 103, S96–S119 (2010).Article 

Google Scholar 
7.Traynor, K. S. et al. Varroa destructor: a complex Parasite, Crippling Honey Bees Worldwide. Trends Parasitol. 36, 592–606 (2020).CAS 
Article 

Google Scholar 
8.Neumann, P. & Carreck, N. L. Honey bee colony losses. J. Apic. Res. 49, 1–6 (2010).Article 

Google Scholar 
9.Thompson, C. E., Biesmeijer, J. C., Allnutt, T. R., Pietravalle, S. & Budge, G. E. Parasite Pressures on Feral Honey Bees (Apis mellifera sp.). PLoS ONE 9, (2014).10.Camazine, S. Differential Reproduction of the Mite, Varroa jacobsoni (Mesostigmata: Varroidae), on Africanized and European Honey Bees (Hymenoptera: Apidae). Ann. Entomol. Soc. Am. 79, 801–803 (1986).Article 

Google Scholar 
11.Corrêa-Marques, M.-H. & De Jong, D. Uncapping of worker bee brood, a component of the hygienic behavior of Africanized honey bees against the mite Varroa jacobsoni Oudemans. Apidologie 29, 283–289 (1998).Article 

Google Scholar 
12.Allsopp, M. H. Analysis of Varroa destructor infestation of southern African honeybee populations. (University of Pretoria, 2007).13.Locke, B., Le Conte, Y., Crauser, D. & Fries, I. Host adaptations reduce the reproductive success of Varroa destructor in two distinct European honey bee populations. Ecol. Evol. 2, 1144–1150 (2012).Article 

Google Scholar 
14.Oddie, M. A. Y., Dahle, B. & Neumann, P. Norwegian honey bees surviving Varroa destructor mite infestations by means of natural selection. PeerJ 5, e3956 (2017).Article 

Google Scholar 
15.Locke, B. Natural Varroa mite-surviving Apis mellifera honeybee populations. Apidologie 47, 467–482 (2016).Article 

Google Scholar 
16.Villegas, A. J. & Villa, J. D. Uncapping of pupal cells by European bees in the United States as responses to Varroa destructor and Galleria mellonella. J. Apic. Res. 45, 203–206 (2006).Article 

Google Scholar 
17.Spivak, M. & Gilliam, M. Facultative expression of hygienic behaviour of honey bees in relation to disease resistance. J. Apic. Res. 32, 147–157 (1993).Article 

Google Scholar 
18.Le Conte, Y., Arnold, G. & Desenfant, P. Influence of brood temperature and hygrometry variations on the development of the honey bee Ectoparasite Varroa jacobsoni (Mesostigmata: Varroidae). Environ. Entomol. 19, 1780–1785 (1990).Article 

Google Scholar 
19.Kraus, B. & Velthuis, H. H. W. High humidity in the honey bee (Apis mellifera L.) Brood nest limits reproduction of the parasitic mite varroa jacobsoni oud. Naturwissenschaften 84, 217–218 (1997).ADS 
CAS 
Article 

Google Scholar 
20.Harris, J. W., Danka, R. G. & Villa, J. D. Changes in infestation, cell cap condition, and reproductive status of varroa destructor (Mesostigmata: Varroidae) in brood exposed to honey bees with varroa sensitive hygiene. Ann. Entomol. Soc. Am. 105, 512–518 (2012).Article 

Google Scholar 
21.Oddie, M. A. Y. et al. Rapid parallel evolution overcomes global honey bee parasite. Sci. Rep. 8, 1–9 (2018).ADS 
CAS 
Article 

Google Scholar 
22.Martin, S. J. et al. Varroa destructor reproduction and cell re-capping in mite-resistant Apis mellifera populations. Apidologie 51, 369–381 (2020).CAS 
Article 

Google Scholar 
23.Martin, S. J. Ontogenesis of the mite Varroa jacobsoni Oud. in worker brood of the honeybee Apis mellifera L under natural conditions. Exp. Appl. Acarol. 18, 87–100 (1994).Article 

Google Scholar 
24.Donzé, G., Herrmann, M., Bachofen, B. & Guerin, P. R. M. Effect of mating frequency and brood cell infestation rate on the reproductive success of the honeybee parasite Varroa jacobsoni. Ecol. Entomol. 21, 17–26 (1996).Article 

Google Scholar 
25.Harris, J. W., Danka, R. G. & Villa, J. D. Honey bees (Hymenoptera: Apidae) with the Trait of varroa sensitive hygiene remove brood with all reproductive stages of varroa mites (Mesostigmata: Varroidae). Ann. Entomol. Soc. Am. 103, 146–152 (2010).Article 

Google Scholar 
26.Harris, J. W. & Harbo, J. R. Low sperm counts and reduced fecundity of mites in colonies of honey bees (Hymenoptera: Apidae) resistant to varroa jacobsoni (mesostigmata: Varroidae). J. Econ. Entomol. 92, 83–90 (1999).Article 

Google Scholar 
27.Peck, D. T. & Seeley, T. D. Mite bombs or robber lures? The roles of drifting and robbing in Varroa destructor transmission from collapsing honey bee colonies to their neighbors. PLoS ONE 14, (2019).28.Arathi, H. S., Ho, G. & Spivak, M. Inefficient task partitioning among nonhygienic honeybees, Apis mellifera L, and implications for disease transmission. Anim. Behav. 72, 431–438 (2006).Article 

Google Scholar 
29.Kirrane, M. J. et al. Asynchronous Development of Honey Bee Host and Varroa destructor (Mesostigmata: Varroidae) Influences Reproductive Potential of Mites. J. Econ. Entomol. 104, 1146–1152 (2011).Article 

Google Scholar 
30.Locke, B. & Fries, I. Characteristics of honey bee colonies (Apis mellifera) in Sweden surviving Varroa destructor infestation. Apidologie 42, 533–542 (2011).Article 

Google Scholar 
31.Chantawannakul, P., Ramsey, S., vanEngelsdorp, D., Khongphinitbunjong, K. & Phokasem, P. Tropilaelaps mite: an emerging threat to European honey bee. Curr. Opin. Insect Sci. 26, 69–75 (2018).Article 

Google Scholar 
32.Bates, D., Maechler, M., Bolker, B. & Walker, S. Fitting linear mixed-effects models using lme4. J. Stat. Softw. 67, 1–48 (2015).Article 

Google Scholar 
33.R Core team. R: A language and environment for statistical computing. (R Foundation for Statistical Computing, 2019), https://www.R-project.org. More

1.Ferrarini, A., Dai, J., Bai, Y. & Alatalo, J. M. Redefining the climate niche of plant species: A novel approach for realistic predictions of species distribution under climate change. Sci. Total Environ. 671, 1086–1093 (2019).ADS 
CAS 
Article 

Google Scholar 
2.Ferrarini, A., Alsafran, M. H. S. A., Dai, J. & Alatalo, J. M. Improving niche projections of plant species under climate change: Silene acaulis on the British Isles as a case study. Clim. Dyn. 52, 1413–1423 (2019).Article 

Google Scholar 
3.Walther, G.-R. et al. Ecological responses to recent climate change. Nature 416, 389–395 (2002).ADS 
CAS 
Article 

Google Scholar 
4.Thuiller, W., Lavorel, S., Araujo, M. B., Sykes, M. T. & Prentice, I. C. Climate change threats to plant diversity in Europe. Proc. Natl. Acad. Sci. 102, 8245–8250 (2005).ADS 
CAS 
PubMed 
Article 

Google Scholar 
5.Mousavi Kouhi, S. M. & Erfanian, M. B. Predicting the present and future distribution of medusahead and barbed goatgrass in Iran. Ecopersia 8, 41–46 (2020).
Google Scholar 
6.Alavi, S. J., Ahmadi, K., Hosseini, S. M., Tabari, M. & Nouri, Z. The response of English yew (Taxus baccata L.) to climate change in the Caspian Hyrcanian Mixed Forest ecoregion. Reg. Environ. Change 19, 1495–1506 (2019).Article 

Google Scholar 
7.Huntley, B., Berry, P. M., Cramer, W. & McDonald, A. P. Special paper: Modelling present and potential future ranges of some European higher plants using climate response surfaces. J. Biogeogr. 22, 967 (1995).Article 

Google Scholar 
8.Pearson, R. G. & Dawson, T. P. Predicting the impacts of climate change on the distribution of species: Are bioclimate envelope models useful?: Evaluating bioclimate envelope models. Glob. Ecol. Biogeogr. 12, 361–371 (2003).Article 

Google Scholar 
9.Hällfors, M. H. et al. Assessing the need and potential of assisted migration using species distribution models. Biol. Conserv. 196, 60–68 (2016).Article 

Google Scholar 
10.Kamakhina, G. L. Kopetdagh-Khorassan Flora: Regional Features of Central Kopetdagh. In Biogeography and Ecology of Turkmenistan (eds. Fet, V. & Atamuradov, K. I.) Vol. 72 129–148 (Springer Netherlands, 1994).11.Memariani, F., Zarrinpour, V. & Akhani, H. A review of plant diversity, vegetation, and phytogeography of the Khorassan-Kopet Dagh floristic province in the Irano-Turanian region (northeastern Iran–southern Turkmenistan). Phytotaxa 249, 8 (2016).Article 

Google Scholar 
12.Fet, V. Biogeographic Position of the Khorassan-Kopetdagh. In Biogeography and Ecology of Turkmenistan (eds. Fet, V. & Atamuradov, K. I.) Vol. 72 197–204 (Springer Netherlands, 1994).13.Memariani, F. Khorassan-Kopet Dagh mountains. In Plant Biogeography and Vegetation of High Mountains of Central and South-West Asia (ed. Noroozi, J.) (Springer, 2020). https://datadryad.org/stash/dataset/doi:10.5061/dryad.4sb638314.Behroozian, M., Ejtehadi, H., Peterson, A. T., Memariani, F. & Mesdaghi, M. Climate change influences on the potential distribution of Dianthus polylepis Bien. ex Boiss. (Caryophyllaceae), an endemic species in the Irano-Turanian region. PLoS ONE 15, e0237527 (2020).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
15.Erfanian, M. B. et al. Data from: Plant community responses to environmentally-friendly piste management in northeast Iran. Dryad Dataset. https://datadryad.org/stash/dataset/doi:10.5061/dryad.4sb6383 (2019).16.Jamzad, Z. Flora of Iran vol. 76 Lamiaceae. (Research Institute of Forests & Rangelands, 2012).17.Sagharyan, M., Ganjeali, A. & Cheniany, M. Investigating the effect of antioxidant compounds and various concentrations of BAP and NAA on the improvement of in vitro stem and root formation of Nepeta binaloudensis Jamzad. NBR 6, 198–205 (2019).Article 

Google Scholar 
18.Nadjafi, F., Koocheki, A., Moghaddam, P. R. & Rastgoo, M. Ethnopharmacology of Nepeta binaludensis Jamzad a highly threatened medicinal plant of Iran. J. Med. Plants 8, 29–35 (2009).
Google Scholar 
19.Nadjafi, F., Koocheki, A., Honermeier, B. & Asili, J. Autecology, ethnomedicinal and phytochemical studies of Nepeta binaludensis Jamzad a highly endangered medicinal plant of Iran. J. Essent. Oil Bear. Plants 12, 97–110 (2009).CAS 
Article 

Google Scholar 
20.Memariani, F., Akhani, H. & Joharchi, M. R. Endemic plants of Khorassan-Kopet Dagh floristic province in Irano-Turanian region: Diversity, distribution patterns and conservation status. Phytotaxa 249, 31 (2016).Article 

Google Scholar 
21.Salmaki, Y. & Joharchi, M. R. Phlomoides binaludensis (Phlomideae, Lamioideae, Lamiaceae), a new species from northeastern Iran. Phytotaxa 172, 265 (2014).Article 

Google Scholar 
22.Pahlevani, A. H., Liede-Schumann, S. & Akhani, H. Seed and capsule morphology of Iranian perennial species of Euphorbia (Euphorbiaceae) and its phylogenetic application: Perennial Species of Euphorbia in Iran. Bot. J. Linn. Soc. 177, 335–377 (2015).Article 

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

Google Scholar 
24.Djamali, M. et al. Application of the global bioclimatic classification to Iran: Implications for understanding the modern vegetation and biogeography. Ecol. Mediterr. 37, 91–114 (2011).Article 

Google Scholar 
25.Farashi, A., Shariati, M. & Hosseini, M. Identifying biodiversity hotspots for threatened mammal species in Iran. Mamm. Biol. 87, 71–88 (2017).Article 

Google Scholar 
26.Hosseinzadeh, M. S., Fois, M., Zangi, B. & Kazemi, S. M. Predicting past, current and future habitat suitability and geographic distribution of the Iranian endemic species Microgecko latifi (Sauria: Gekkonidae). J. Arid Environ. 183, 104283 (2020).ADS 
Article 

Google Scholar 
27.Noroozi, J. et al. Endemic diversity and distribution of the Iranian vascular flora across phytogeographical regions, biodiversity hotspots and areas of endemism. Sci. Rep. 9, 12991 (2019).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
28.Erfanian, M. B., Ejtehadi, H., Vaezi, J. & Moazzeni, H. Plant community responses to multiple disturbances in an arid region of northeast Iran. Land Degrad. Dev. 30, 1554–1563 (2019).Article 

Google Scholar 
29.Erfanian, M. B. et al. Plant community responses to environmentally friendly piste management in northeast Iran. Ecol. Evol. 9, 8193–8200 (2019).PubMed 
PubMed Central 
Article 

Google Scholar 
30.Memariani, F. et al. Plant diversity of the Khorassan-Kopet Dagh Floristic Province (Irano-Turanian Region). (Magnolia Press, 2016)31.Memariani, F., Joharchi, M. R., Ejtehadi, H. & Emadzade, K. A contribution to the flora and vegetation of Binalood mountain range, NE Iran: Floristic and chorological studies in Fereizi region. Ferdowsi Univ. Int. J. Biol. Sci. J. Cell Mol. Res. 1, 1–17 (2009).
Google Scholar 
32.Memariani, F. & Joharchi, M. R. Iris ferdowsii (Iridaceae), a new species of section Regelia from northeast of Iran. Phytotaxa 291, 192 (2017).Article 

Google Scholar 
33.Thuiller, W., Georges, D., Engler, R. & Breiner, F. biomod2: Ensemble Platform for Species Distribution Modeling. R Package. https://cran.r-project.org/package=biomod2 (2019).34.R Core Team. R: A Language and Environment for Statistical Computing. (R Foundation for Statistical Computing, 2020).35.Aiello-Lammens, M. E., Boria, R. A., Radosavljevic, A., Vilela, B. & Anderson, R. P. spThin: An R package for spatial thinning of species occurrence records for use in ecological niche models. Ecography 38, 541–545 (2015).Article 

Google Scholar 
36.Fick, S. E. & Hijmans, R. J. WorldClim 2: New 1-km spatial resolution climate surfaces for global land areas. Int. J. Climatol. 37, 4302–4315 (2017).Article 

Google Scholar 
37.Ahmadi, M., Dadashi Roudbari, A. A., Akbari Azirani, T. & Karami, J. The performance of the HadGEM2-ES model in the evaluation of seasonal temperature anomaly of Iran under RCP scenarios. J. Earth Space Phys. 45, 625–644 (2019).
Google Scholar 
38.Dray, S. & Dufour, A.-B. The ade4 Package: Implementing the duality diagram for ecologists. J. Stat. Softw. 22, 1–20 (2007).Article 

Google Scholar 
39.Guisan, A., Thuiller, W. & Zimmermann, N. E. Habitat Suitability and Distribution Models: With Applications in R. (Cambridge University Press, 2017).40.Naimi, B., Hamm, N. A. S., Groen, T. A., Skidmore, A. K. & Toxopeus, A. G. Where is positional uncertainty a problem for species distribution modelling. Ecography 37, 191–203 (2014).Article 

Google Scholar 
41.Menard, S. W. Applied Logistic Regression Analysis (Sage Publications, Thousand Oaks, 2002).Book 

Google Scholar 
42.Landis, J. R. & Koch, G. G. The measurement of observer agreement for categorical data. Biometrics 33, 159 (1977).CAS 
MATH 
Article 

Google Scholar 
43.Araujo, M. & New, M. Ensemble forecasting of species distributions. Trends Ecol. Evol. 22, 42–47 (2007).PubMed 
Article 

Google Scholar 
44.Breiner, F. T., Guisan, A., Bergamini, A. & Nobis, M. P. Overcoming limitations of modelling rare species by using ensembles of small models. Methods Ecol. Evol. 6, 1210–1218 (2015).Article 

Google Scholar 
45.Kaky, E., Nolan, V., Alatawi, A. & Gilbert, F. A comparison between Ensemble and MaxEnt species distribution modelling approaches for conservation: A case study with Egyptian medicinal plants. Ecol. Inform. 60, 101150 (2020).Article 

Google Scholar 
46.Hao, T., Elith, J., Lahoz-Monfort, J. J. & Guillera-Arroita, G. Testing whether ensemble modelling is advantageous for maximising predictive performance of species distribution models. Ecography 43, 549–558 (2020).Article 

Google Scholar 
47.Abdelaal, M., Fois, M., Fenu, G. & Bacchetta, G. Using MaxEnt modeling to predict the potential distribution of the endemic plant Rosa arabica Crép, Egypt. Ecol. Inform. 50, 68–75 (2019).Article 

Google Scholar 
48.Thuiller, W. et al. Endemic species and ecosystem sensitivity to climate change in Namibia. Glob. Change Biol. 12, 759–776 (2006).ADS 
Article 

Google Scholar 
49.Chitale, V. S., Behera, M. D. & Roy, P. S. Future of endemic flora of biodiversity hotspots in India. PLoS ONE 9, e115264 (2014).ADS 
PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
50.Fois, M., Bacchetta, G., Cogoni, D. & Fenu, G. Current and future effectiveness of the Natura 2000 network for protecting plant species in Sardinia: A nice and complex strategy in its raw state?. J. Environ. Plan. Manag. 61, 332–347 (2018).Article 

Google Scholar 
51.Mamet, S. D., Brown, C. D., Trant, A. J. & Laroque, C. P. Shifting global Larix distributions: Northern expansion and southern retraction as species respond to changing climate. J. Biogeogr. 46, 30–44 (2019).Article 

Google Scholar 
52.Thuiller, W., Lavorel, S. & Araújo, M. B. Niche properties and geographical extent as predictors of species sensitivity to climate change: Predicting species sensitivity to climate change. Glob. Ecol. Biogeogr. 14, 347–357 (2005).Article 

Google Scholar 
53.Hosseini, S. S., Ejtehadi, H. & Memariani, F. The first report Nepeta binaloudensis Jamzad in Hezar masjed mountains of Khorasan Razavi province. In Proceedings of the 9th National Congress and 7th International Congrees of Bilogy of Iran (2016).54.Dullinger, S. et al. Extinction debt of high-mountain plants under twenty-first-century climate change. Nat. Clim. Change 2, 619–622 (2012).ADS 
Article 

Google Scholar 
55.Wiens, J. J. Climate-related local extinctions are already widespread among plant and animal species. PLoS Biol. 14, e2001104 (2016).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
56.Casazza, G. et al. Climate change hastens the urgency of conservation for range-restricted plant species in the central-northern Mediterranean region. Biol. Conserv. 179, 129–138 (2014).Article 

Google Scholar 
57.Zhang, M.-G. et al. Major declines of woody plant species ranges under climate change in Yunnan, China. Divers. Distrib. 20, 405–415 (2014).CAS 
Article 

Google Scholar 
58.Sanjerehei, M. M. & Rundel, P. W. The impact of climate change on habitat suitability for Artemisia sieberi and Artemisia aucheri (Asteraceae)—A modeling approach. Pol. J. Ecol. 65, 97–109 (2017).Article 

Google Scholar 
59.Abolmaali, S.M.-R., Tarkesh, M. & Bashari, H. MaxEnt modeling for predicting suitable habitats and identifying the effects of climate change on a threatened species, Daphne mucronata, in central Iran. Ecol. Inform. 43, 116–123 (2018).Article 

Google Scholar 
60.Di Musciano, M. et al. Dispersal ability of threatened species affects future distributions. Plant Ecol. 221, 265–281 (2020).Article 

Google Scholar 
61.Fois, M., Cuena-Lombraña, A., Fenu, G., Cogoni, D. & Bacchetta, G. The reliability of conservation status assessments at regional level: Past, present and future perspectives on Gentiana lutea L. ssp. lutea in Sardinia. J. Nat. Conserv. 33, 1–9 (2016).Article 

Google Scholar  More

Plant materialsIn April 2019, 2-year-old Chinese fir seedlings were collected from the Yalin Center of the Chinese Academy of Forestry in good condition and free from pests and diseases. (The use of Chinese fir seedlings in the experiment complies with national regulations).MediumPikovskava (PVK) solid medium: glucose 10 g, Ca3(PO4) 25 g, CaCO35g, (NH4)2SO40.5 g, NaCl 0.2 g, MgSO4 7H2O 0.1 g, KCl 0.1 g, MnSO4 0.002 g, FeSO4 7H2O 0.002 g, agar 18 g, Distilled water 1000 mL, pH 7.0.PVK liquid medium: PVK solid medium without agar.Luria-Bertan (LB) medium: tryptone 10 g, yeast extract powder 5 g, NaCl 10 g, agar 18 g, distilled water 1000 mL, pH7.0.LB liquid medium:LB solid medium without agar.Isolation and purification of endophytes from Chinese fir seedlingsThe roots, stems, and leaves of the selected Chinese fir seedlings were washed away with running water to remove the surface soil, and then washed with running water for 24 h to 36 h, and the surface moisture was absorbed with sterile filter paper. Weigh 1 g of roots, stems and leaves in a petri dish, and then carry out surface disinfection in a sterile operating table with 75% alcohol (C2H5OH) for 30 s, 5% sodium hypochlorite (NaClO) 10 min, wash the sterile water 7 times, and use sterile filter paper to absorb the water. The material after surface sterilization is cut into 2 mm × 2 mm with sterile surgical scissors, placed in a sterilized mortar and grated with a small amount of sterile quartz sand and ground into a homogenate, then diluted with sterile water to 10–1, 10–2 and 10–3, pipette to draw 150 µL of sample grinding fluid, spread on the medium, set the sterile water of the last rinse of Chinese fir tissue as a blank control, incubate with other plates under the same conditions, and verify whether the surface disinfection. Cultivate in a 28 ℃ incubator according to the characteristics of the colony phenotype, and use the streak separation method to further purify and isolate the strains until the isolation of the colony morphology is uniform for each isolate.Note: 75% C2H5OH: 75 mL absolute ethanol + 25 mL sterile deionized water; 5% NaClO: sodium hypochlorite solution with 10% available chlorine: sterile deionized water = 1:1.Characterization of PGP traitsDetermination of phosphorus solubilizing abilityThe strain was inoculated into PVK liquid medium and cultured at 28 ℃ for 180 days/min on a reciprocating shaker for 7 days, and then the pH value of the medium was measured. The culture solution was centrifuged at 8000 rpm for 15 min to remove bacterial cells. Take the supernatant and use the molybdenum antimony scandium colorimetric method23 to determine the soluble phosphorus content in the culture broth.Determination of nitrogenase activityAn aliquot of 200 µL fresh culture was inoculated to 20 mL of nutrient broth and incubated overnight at 30℃. Bacterial growth was collected by centrifugation and was washed twice using sterile water, and resuspended by liquid limited nitrogen culture medium (OD600 = 0.2). The 3 mL suspension was transferred to a 25 mL sterilized serum vial and 2.4 mL acetylene gas (99.9999%) was driven into the serum bottle, and then incubated at 30 °C for 12 h. The ethylene content and the protein of bacterial suspension were determined as You et al.24.1-Aminocyclopropane-1-carboxylate (ACC) deaminase activity determinationACC deaminase activity was determined by the method of Glick et al.25 using N-free medium (Nfb)26 for bacteria and minimal medium (MM)27 for actinomycetes containing 0.3 m mol L−1 ACC (Sigma, USA) as a sole nitrogen source. MM with 0.1% (w/v) NH4(SO4)2 was used as a positive control and cultivation without ACC was used as a negative control. After incubation at 28 ℃ for 7 days for non-actinomycete bacteria and 14 days for actinomycetes, colony growth on Nfb or MM with addition of ACC indicated ACC deaminase activity.Indole-3-acetic acid (IAA) productionIAA production was measured by colorimetric assay27. Bacterial isolates were cultured for 3 days in TY broth (without L-tryptophan or supplemented with 500 μg/mL of l-tryptophan) in the dark at 28 °C. Cells were removed from the culture medium by centrifugation at 13,000×g for 10 min; then, 1 mL of the supernatant was mixed vigorously with 2 mL of Salkowski’s reagent (4.5 g of FeCl3 per L in 10.8 M H2SO4). Samples were incubated at room temperature for 30 min and the IAA production was estimated from the optical density at 600 nm (OD600) by comparison with a standard curve prepared from known concentrations of IAA.Siderophore productionSiderophore production was examined by using chrome azurol S (CAS) agar28. Isolate was inoculated onto CAS agar, cultured at 28 °C for 2 days, and the positive strain was indicated by an orange halo around the bacterial colony. Determine the ratio (D/d) of orange aperture diameter (D) to colony diameter (d) to determine the iron-producing carrier capacity of the strain.Physiological and biochemical tests of phosphate-solubilizing bacteriaThe conventional physiological and biochemical identification of PSB is carried out according to the methods in the “Common Bacterial System Identification Manual”, which mainly includes Gram stain, glucose hydrolysis test, lactose hydrolysis test, methyl red test, Voges-Proskauer (VP) test, hydrogen sulfide production test, gelatin liquefaction Test, citrate utilization test, malonate utilization test, denitrification test.16 SrRNA gene sequencingTaking the screened multifunctional PSB as the object, the bacterial genomic DNA extraction kit of Beijing Bomed Biotechnology Co., Ltd. was used to extract the DNA of the strain, using the DNA as a template, and using the bacterial universal primer 27F (5′-AGAGTTTGATCCTGGCTCAG-3′) and 1492R (5′-GGTTACCTTGTTACGACTT-3′) PCR amplification, the amplification system is as follows: DNA template 1 uL, primer 27F 0.5 µL, 1492R 0.5 µL, 2 × TaqMix 12.5 µL, ddH2O10.5 µL. The PCR procedure is as follows: 93℃ for 3 min, 93℃ for 30 s, 56 ℃ for 30 s, 72℃ for 2 min, 32 cycles; 72 ℃ for 7 min. The amplified products were sequenced bidirectionally by BGI. After splicing the measured 16SrDNA sequences in ContigExpress, search in GenBank, EzTaxon, BIGSdb databases respectively, select the model strains with high homology. The phylogenetic tree was constructed by the neighbor-joining method using MEGA version 7.0 with the Kimura 2-parameter model29, the robustness of the tree was evaluated by performing bootstrap analyses based on 1000 replications30.Evaluation of plant growth promotion by individual inoculationPreparation of inoculumThe single colonies were picked out and incubated in LB broth at 180 rpm at 28 ℃ for 12 h. The above solution was inserted into 200 mL of LB broth at 1% inoculation, incubated at 180 rpm at 28℃ for 48 h and the cell pellet was resuspended in sterile distilled water and made up to a final concentration of 3 × 108 CFU/mL.Experimental seedlings and soilChinese fir seedlings were provided by the Experimental Center of Subtropical Forestry, Chinese Academy of Forestry (117° 67′ E, 27° 82′ N), Jiangxi Province, China. The use of Chinese fir seedlings in the experiment complies with national regulations. Five months old seedlings with vigorous and apparently disease and pest free were used. The height and root collar diameters of seedlings were 8.7 cm and 1.36 mm, respectively. The seedling container was made of non-woven fabric, the specification was 4.5 cm × 8.0 cm (Diameter × Height). The soil was consisted of nursery medium and loess at a ratio of 9:1, was thoroughly mixed and homogenized with 3 kg slow-release fertilizer per cubic. The slow-release fertilizer is produced by American Abbes (180 g kg−1 total N, 80 g kg−1 available P, and 80 g kg−1 total K, the fertilizer effect period is 9 months). The soil exhibited the following properties: 6.34 g kg−1 total N, 0.80 g kg−1 total P, 2.50 g kg−1 total K, and a pH-value of 6.00.Test designThe SSP2, JRP22 and HRP2, which were confirmed to have the characteristics of promoting plant growth, so pot experiments were conducted. To determine the effectiveness of phosphorus-solubilizing bacteria in plant growth of Chinese fir, a pot culture experiment was conducted between August and November 2019 in an open-sided greenhouse in Experimental Center of Subtropical Forestry, Chinese Academy of Forestry, Jiangxi Province, China. The experiment was carried out in three-factor orthogonal design with five replications for each treatment. The orthogonal experimental design of experiment is provided in Table 1, and strain, dilution ratio, inoculation method contained 3 levels. The pots with water was used as control (CK). For the irrigation of the rhizosphere (IR) treatments, 30 mL of diluted inoculum was added to the soil in the vicinity of the roots of Chinese fir. For the foliar spray (FS) treatments, 30 mL diluted bacterial cell suspension was inoculated in the leaves of Chinese fir by using a syringe. For the rhizosphere + foliar spray (IS) treatments, 15 mL of diluted inoculum was inoculated in the rhizosphere of seedlings, 15 mL was added to the leaves of seedlings by foliar spray. Each treatment contained sixteen seedlings for a total of nine treatments. A total of 3 inoculations were given in the middle of each month. The plant height and stem diameter were recorded before the first inoculation. The plants were harvested after 90 days (16 plantlets/replicate/treatment, i.e., a total of 80 plantlets per treatment) and the root biomass, stem biomass, leaf biomass, plant height and stem diameter were measured.Table 1 L9(34) Orthogonal design of experiment.Full size tableDetermination of growth indicatorsDuring the test, before each inoculation, the height of the seedlings was measured with a ruler and the ground diameter of the seedlings was measured with a vernier caliper. After the experiment, 10 plants of Chinese fir seedlings in average growth were randomly selected from each treatment, a total of 30 plants were washed with clean water to remove surface impurities, the filter paper was dried and the roots, stems, and leaves were put into paper bags respectively at 105 °C. After being degraded for 0.5 h, dried at 70 °C to a constant weight, weighed and recorded the biomass of each part.Determination of leaf and soil nutrient contentTotal N content of leaf and soil was measured by a 2300 Kjeltec Analyzer Unit (FOSS, Höganäs, Sweden). Total P and total K, total Mg and total Fe of leaf, soil TP, TK, AP and AK were extracted according to literature31 and were determined by ICP (Kleve, Germany).Determination of soil enzyme activitiesActivities of the soil urease, cellulase, sucrase, dehydrogenase and acid phosphatase were determined by spectrophotometry. Firstly, 0.05 g of soil was added to 450 mL of phosphate buffer solution (PBS, 0.1 mol L−1, pH 7.4). Then the solution was mixed by shaking, and centrifuged at 2000 rpm at 4 ℃ for 10 min and supernatant was collected with a new centrifugal tube. The supernatant and reagents were added according to the kit instructions (Shanghai Enzyme-linked Biotechnology Co., Ltd., Shanghai, China). Absorbance at 450 nm was measured on a SpectraMax Paradigm Multi-Mode detection platform (Molecular devices, San Jose, CA, USA).Data processing and analysisAll statistical analyses were performed using SPSS24. Data are presented in terms of means (± SE; standard error). Statistical differences were tested by one-factor ANOVA to evaluate the differences in the nutrient content of soil and plant growth status. In MEGA7.0, the Neighbor-Joining method was used to construct the phylogenetic tree, and the Bootstrap value was 1000. More