Philippa Kaur

Site descriptionWe define “kelp forest” as rocky-reef habitat within the 5–20 m depth range that supports dense stands of giant kelp, Macrocystis pyrifera. For this study, we considered the Santa Barbara Channel (SBC) to include the mainland region between Point Conception (−120.476° longitude, 34.455° latitude) and Point Mugu (−119.065° longitude, 34.079° latitude), as well the northern and southern sides of the four northern Channel Islands (Fig. 1). Although the SBC is a subset of the Southern California Bight, its strong west-east gradient in cold to warm temperature means the study system includes many of the kelp-forest species throughout California31. This means the SBC kelp-forest food web is a large “metaweb”, characterizing kelp forest meta-communities, rather than a site-specific web. In other words, the system includes cold water and warm water species that might not necessarily co-occur at a single site. However, there are site-specific food webs embedded in the metaweb at particular locations where a subset of species occur.Data sourcesOur goal was to assemble the food web using both published and novel empirical observations. To this end, we first used published data sets and species’ range boundaries to create free-living species lists. The initial list of fishes, algae, and free-living invertebrates was assembled from the Channel Islands National Park Kelp Forest Monitoring program (CINP KFM, annual reports available at https://irma.nps.gov/DataStore/SavedSearch/Profile/1508, accessed March 6, 2017, or visit https://www.nps.gov/im/medn/kelp-forest-communities.htm to contact David Kushner or Joshua Sprague) and the SBC Long Term Ecological Research program’s ongoing kelp-forest community timeseries (SBC LTER, https://sbclter.msi.ucsb.edu/data/catalog/, accessed March 12, 2017). We added to these lists using primary literature, technical reports (e.g., NOAA, USFW), direct observations, expert opinion, crowd-sourced observations (e.g., eBird.org), guidebooks, and grey literature. We sampled the local kelp forest zooplankton and the algae-associated small-invertebrate community, because these organisms were not well represented in surveys (see below).We created initial lists of parasite species using published literature and host-parasite databases. A systematic review was conducted to collect parasite records for each free-living species. We searched the Natural History Museum (NHM) of London host-parasite database (https://www.nhm.ac.uk/research-curation/scientific-resources/taxonomy-systematics/host-parasites/database/search.jsp), the FishPest database32, WoRMs (http://www.marinespecies.org/aphia.php?p = search), BIOSIS citation index (http://webofscience.com), and Google Scholar™(https://scholar.google.com/) (Genus + species + parasit*, expanded to Genus + parasit* if no records were found). For each host species, we recorded the number of records found in BIOSIS and NHM as an estimate of study effort. Although parasites are often reported at the host and parasite species level, we were often able to infer parasite and host life stages based on knowledge about life cycles. We added to these lists by sampling local fish and invertebrates, with a focus on hosts that were common in the system and not well-studied (see below). As for any food-web study, we were most interested in including common or important parasites, rather than rarities.Published diet observations (including in grey literature), direct observations, and inference were used to determine trophic links (see below).Free-living species sampling methodsCertain groups of free-living species were under-represented in published survey data, so we conducted sampling to assess species diversity in the following areas.Zooplankton towsWe conducted vertical zooplankton tows within kelp forests at two island locations (on the same date) and two mainland locations (repeated tows, four dates at one site, three of those dates at a second site, including one nighttime sampling date), for eight site by date samples30. While the vessel was at anchor within a kelp forest, a 30 cm diameter, 200 micron plankton net was dropped to the bottom and pulled to the surface at a rate of 0.33 m per second. Care was taken not to scrape the net against kelp plants. The collection jar attached to the net was kept vertical with a small lead weight to ensure that the net did not collect organisms on the way down to the bottom. The depth and time of collection were recorded30. We held collection jars on ice while in the field, then preserved specimens in 95% ethanol when we returned to the lab (within a few hours of collection). All organisms were counted and identified to species when possible, but some groups were identified to Order or Family, and then cross-checked with lists of known local species. If this was not possible, specimens were assigned to morphospecies, indicating they appeared to be a unique species based on morphology. Representative specimens from each species or morphospecies were photographed and measured.Giant kelp holdfastsGiant kelp holdfasts were sampled for free-living invertebrates. In the field, holdfast circumference and two slant height measures were taken, as well as basal stipe circumference. A subsample of approximately 25% of the holdfast was collected by SCUBA in a large plastic zip bag, and frozen until processing (n = 7). The samples were processed for organisms  > 500 microns, and holdfast tissue was weighed after organisms and debris were removed. Organisms were counted, identified to species or morphospecies when possible, and measured30. Some groups were identified to Family, and then matched to lists of known local species.Taxon-specific methods: gastropodsSmall gastropods are a diverse but overlooked group that lives in benthic turf algae. Algal clumps were collected haphazardly by either laying down a 7 × 7 cm quadrat and collecting all algae within the quadrat, or by collecting clumps of a particular alga and weighing at the lab. All gastropods were removed by hand under a stereomicroscope, counted, identified to species or morphospecies, measured, and photographed30.Parasitological collectionsWe collected fish and invertebrates and dissected them for parasites, with the goal of identifying the most common parasites in the food web. We targeted host groups that are known to transmit trophically-transmitted parasites in other systems. We collected most organisms from mainland sites, and sampled opportunistically at sites on Anacapa, Santa Cruz, and Santa Rosa islands30 (Fig. 2). A list of all species dissected and sample sizes is provided30.Fish collectionsWe prioritized collecting the most common and abundant fish species based on survey data from 2000–2014 (SBC LTER), as well as personal observation, expert opinion, and amount of parasite data in the literature. Other species (lower abundance or higher past study effort) were collected opportunistically. Fish were collected primarily by spear on SCUBA. Specific size classes were not targeted and the spear tips used were appropriate for the focal species. Small benthic fish were collected using dip nets. All fish were collected under UCSB IACUC protocol 549.2. Fish were either stored on ice and processed within 24 hours of collection or frozen until processing.Invertebrate collectionsInvertebrates are necessary intermediate hosts in many parasite life cycles, but relatively few parasite life cycles have been described in marine environments. We targeted invertebrate species that were abundant and potentially important as intermediate hosts for parasites. We did not collect sessile colonial taxa, such as hydroids, gorgonians, sponges, and tunicates, as they were not expected to be hosts for trophically transmitted parasites (but these hosts do merit further study). Most sampled invertebrates were gastropods and small crustaceans, as they host trophically-transmitted parasites in other food webs. Bivalves, large crustaceans, echinoderms, and polychaetes were also dissected. Large invertebrates were collected by hand or using a rock chisel and scraper when appropriate. Small invertebrates were sampled by collecting benthic substrates in plastic or fine mesh bags and removing organisms in the lab. Invertebrates were held live in flow-through seawater until the time of dissection or frozen until processing.Parasitological assessmentFor each host dissection, the exterior and all internal soft tissues were examined for parasite life stages. For larger species, entire host organs were usually searched by pressing soft tissues thin between two glass plates (“squashed”) and examining with a stereomicroscope. However, to increase sample size, bilaterally symmetric organs (e.g. gills) were examined from one randomly determined side, and large organs (e.g. muscle, liver) were subsampled in larger fishes. Small crustaceans and soft-bodied invertebrates were squashed whole. We identified gut contents where feasible to improve host diet data and inform parasite life cycles. We recorded host mass, length (or other species-appropriate measurement), collection method, and host condition at time of dissection (e.g. frozen, fresh). We counted and identified all parasites to the lowest possible taxonomic level and assigned a morphospecies code when species-level identification was not possible. Only a few putative parasites were excluded from additional analysis because they had no identifying features. Dissection data30 includes species not included in the full food web (see below for discussion of justifications for node inclusion).Node list assemblyNodes in the web included free-living species that used the water column and benthic zones within kelp forests as feeding habitat (including transient kelp-forest visitors but excluding rare and vagrant species) and parasites of those free-living species. Species was the preferred taxonomic unit, and life stages were included as separate nodes if that life stage was present in the system and had distinct trophic interactions from the adult stage. The fully-resolved free-living food web was constructed with life stage (e.g., larva, adult) nested within species (or morpho-species) (excepting benthic diatoms, planktonic diatoms, dinoflagellates, foraminifera, free-living nematodes, bacteria, free-living ciliates, copepod nauplii, filamentous algae, and invertebrate eggs, which were aggregate nodes). As various forms of detritus are important to energy flow in kelp forests, detritus was broken into four categories based on the typical feeding modes of detritivores and main sources of detritus: carrion, drift macroalgae, small mixed origin (such as would be consumed by a deposit or suspension feeder, with the recognition that this alone is a complex system deserving further resolution) and dissolved organic material. The “drift macroalgae” component was especially important to distinguish, as certain herbivores (sea urchins) are known to prefer drift algae as food but will turn to feeding on live algae when drift algae are sparse. This is a very distinct type of interaction from suspension feeders, which consume small particles of detritus that may be largely bacteria. “Parasites” are consumers which fit the seven types of parasitism defined by Lafferty and Kuris33. Commensal organisms were also recorded. We limited the parasite species list to metazoan species that use kelp-forest species as hosts for at least one stage in their life cycle. Bacterial, viral, fungal, and protozoan pathogens that are important in kelp-forest food webs merit inclusion in further work.We assigned each node a justification code (see below), confidence level, literature reference, and locality of the reference. Additional node metadata includes site on host (ecto-vs. endoparasite), taxonomic information, and life cycle information30 (see below). The node list contains columns with a species ID, and a species-by-stage ID. To work with the life-stage resolution, select the species-by-stage ID as the node identifier in analyses. To work with the species version, select the species ID as the node identifier in analyses. This will collapse all of the interactions to the species, so all of the trophic interactions are preserved and linked to the species node. Network analysis packages in R (such as Cheddar34) will automatically remove duplicate links if they are generated in this process.Life stages as nodesSpecies were partitioned into life-stage nodes (e.g., larva, juvenile, adult) if a species changed its trophic position from one stage to the other and multiple stages were present in the system. Whether or not a distinct life stage resided in the kelp forest was indicated by various data sources (e.g. dissections, published records), or inferred from species life history or trophic interactions. For example, amphipods brood offspring and have crawl-away juveniles. These juveniles remain in the kelp forest (rather than having a pelagic phase), and due to their small size are subject to different predators than adults (e.g. adults are eaten by fishes, while juveniles are eaten by hydroids). This was justification for juvenile amphipods being a distinct node from adult amphipods. On the other hand, many species have planktonic larvae that develop outside of the kelp forest, so only the adult stages were included at the species level. Larval stages of parasites were included if there was no feasible alternative for the focal host to become infected. We assumed that kelp-forest resident hosts became infected through life-cycle stages found within the kelp-forest food web, but that transient hosts could have acquired some parasites outside the kelp forest (e.g., if intermediate hosts were not known from the kelp forest). Likewise, presence of larval parasites in dissections was evidence for including adult stages. For some species, there was insufficient data on life history to infer additional stages. Metadata in the node list indicates whether parasites have additional life stages inside the kelp forest, outside, or unknown. When comparing this food web with others (which rarely separate species into life stages), using our data it is easy to collapse life-stage nodes into species nodes.Justifications for node inclusionWe used multiple lines of evidence to justify whether or not to include a node in the food web. Free-living species were included if they were known from the SBC (see site description above) and were indicated by the data sources described above (e.g. reports, surveys, published papers, guidebooks, expert opinion, etc.). Species lists from regional guidebooks included non-kelp-forest species, so these lists were compared with species lists from long-term monitoring surveys. Following the methods of Lafferty et al. 2006, we excluded most rare species (0.5, we assumed that an unobserved link actually occurred unless otherwise contradicted by species life history. We also then noted the probability of a false positive link (1 – ({widehat{F}}_{{rm{ij}}})). We further identified those few host and parasite species that generated substantial error in the network. To keep the overall error rate to More

Single plasmid, single-population modelsTo understand the dynamics of parasitic plasmids in complex ecologies, we first need to understand their behavior in simple scenarios. In this section, we analyze the dynamics of plasmids spreading by different HGT mechanisms in single populations. We begin by modeling competition between plasmid-free cells and cells containing a conjugative plasmid. A nutrient, with concentration (C), is supplied to the system at rate (S). Cells grow at a rate proportional to (C) with proportionality constant (alpha) for plasmid-free cells or ((1 ,-, {Delta})alpha) for plasmid-containing cells. Since we are interested in parasitic plasmids, we assume that ({Delta} in (0,1)). Cells of both types die at a rate (delta). When a plasmid-containing cell divides there is a loss probability, (p_ell), for one of the daughter cells to contain no plasmids. As long as a daughter cell contains at least one plasmid, the original plasmid copy number (the number of copies of the plasmid maintained per cell) is regenerated (as depicted in Fig. 1A). Plasmids can spread horizontally by conjugation, as illustrated in Fig. 1B, wherein a plasmid-free cell and a plasmid-containing cell interact to produce two plasmid-containing cells. We model the rate of conjugation by a mass-action term with rate (gamma _{mathrm{c}}). The equations governing the dynamics of conjugation are therefore:$$ frac{{drho }}{{dt}} ,=, alpha Crho ,-, gamma _{mathrm{c}}rho rho _{mathrm{p}} ,+, p_ell (1 ,-, {Delta})alpha Crho _{mathrm{p}} ,-, delta rho ,\ frac{{drho _{mathrm{p}}}}{{dt}} ,=, (1 ,-, {Delta})alpha Crho _{mathrm{p}} ,+, gamma _{mathrm{c}}rho rho _{mathrm{p}} ,-, p_ell (1 ,-, {Delta})alpha Crho _{mathrm{p}} ,-, delta rho _{mathrm{p}},\ frac{{dC}}{{dt}} ,=, S ,-, alpha Crho ,-, (1 ,-, {Delta})alpha Crho _{mathrm{p}}.$$
(1)
Fig. 1: Different modeled mechanisms of plasmid transfer lead to distinct ecological phase diagrams, but all such mechanisms leave individual populations susceptible to runaway plasmid invasion.A At each division, plasmids are randomly segregated between daughter cells. Original plasmid copy number is regenerated if at least one plasmid remains in a daughter cell. B Schematic of plasmid transfer mechanisms. Left: spread of plasmids by plasmid-containing cells conjugating with plasmid-free cells. Right: spread of plasmids by extracellular plasmids infecting plasmid-free cells via transformation. C Phase diagram for conjugative plasmids as a function of plasmid cost, ({Delta}), and (gamma _{mathrm{c}}); (delta ,=, 0.1), (S ,=, 1), (p_ell ,=, 0), and (alpha ,=, 1) (see Eq. 4). D Phase diagram for transformative plasmids as a function of ({Delta}) and (gamma _{mathrm{t}}). Parameters as in C with (delta _{mathrm{p}} ,=, 0.3) and (n_{{mathrm{eff}}} ,=, 0.6) (see Eq. 9). See “Methods” for details. E In model multiplasmid cells, plasmid types segregate independently. If at least one plasmid of a given type remains in a daughter cell, the full copy number of that plasmid type is regenerated. F Fitness cost as a function of number of unique plasmid types in a cell for multiplicative case ({Delta}_{{mathrm{tot}}} ,=, 1 ,-, (1 ,-, {Delta})^m) with ({Delta} ,=, 0.05). G Steady-state distribution of number of plasmid types per cell at different conjugation rates, measured relative to (gamma _{mathrm{c}}^ ast) (the critical conjugation rate necessary for invasion of a single plasmid into a plasmid-free population, see Eq. 4). Results for eight unique plasmid types with (delta ,=, 1), ({Delta} ,=, 0.1), (alpha ,=, 1), (S ,=, 1), and (p_ell ,=, 0.05).Full size imageIn this model, what are the conditions for a parasitic conjugative plasmid to be able to invade a plasmid-free population? Invasibility implies that the equilibrium containing only plasmid-free cells is locally unstable, which occurs when$$qquadqquadqquadgamma _{mathrm{c}}rho ^ ast , > , delta {Delta} ,+, delta p_ell (1 ,-, {Delta}),$$
(2)
where (rho ^ ast ,=, S/delta) is the steady-state abundance of the plasmid-free cells at the plasmid-free equilibrium. This invasibility condition has an intuitive physical interpretation: to invade, the rate of conjugation must overcome losses due to reduced host growth rate as well as plasmid loss during division. This condition is similar to those found in previous studies [15].Given the condition for plasmid invasion in Eq. 4, what is the optimal behavior for a parasitic conjugative plasmid? The left-hand-side of the expression is linear in the plasmid-free population, meaning that it is more difficult for a plasmid to invade smaller populations. To favor invasion, the plasmid can minimize the right-hand-side of the equation. For a plasmid that relies on random segregation upon cell division, both the plasmid cost ({Delta}) and the loss probability (p_ell) are functions of plasmid copy number, (n_{mathrm{p}}), a property controlled by the plasmid itself. If the primary cost of a plasmid is its replication and its gene products, plasmid cost will scale with copy number such that ({Delta} ,=, {Delta}_{mathrm{p}}n_{mathrm{p}}), where ({Delta}_{mathrm{p}}) is the cost of an individual plasmid copy. The loss probability will be (p_ell ,=, 2^{1 ,-, n_{mathrm{p}}}), i.e., the probability that a daughter cell receives zero plasmids from random segregation. The right-hand-side of the invasion condition Eq. 4 is therefore (delta ({Delta}_{mathrm{p}}n_{mathrm{p}} ,+, 2^{1 ,-, n_{mathrm{p}}}(1 ,-, {Delta}_{mathrm{p}}n_{mathrm{p}}))), which has a minimum at finite (n_{mathrm{p}}). The minimum in the invasion boundary at finite (n_{mathrm{p}}) indicates that in our framework optimal conjugative plasmids have a moderate copy number.What kinds of ecological dynamics does our model for a conjugative parasitic plasmid exhibit? To answer this question, we characterize the stability of the system’s equilibria (see SI Appendix 1 for details). For conjugative plasmids with the optimal copy number, the dominant form of loss will be from reduced host fitness (see SI Fig. S1), and thus we characterize the case of negligible loss rate (p_ell ,=, 0) (we consider the case of finite loss rates in SI Fig. S2 and find similar results). In Fig. 1C we show the phase diagram of possible ecological outcomes as a function of plasmid cost ({Delta}) and conjugation rate (gamma _{mathrm{c}}). For high values of plasmid cost and low values of conjugation rate, the plasmid is unable to invade and the plasmid-free equilibrium is the only stable state. As plasmid cost decreases or conjugation rate increases, plasmids are able to invade and there is a state of stable coexistence between plasmid-free and plasmid-containing cells. The range of conjugation rates permitting coexistence is larger for costlier plasmids. Once the plasmid cost is sufficiently low or the conjugation rate is sufficiently high, the unique stable state consists only of plasmid-containing cells (note that for finite values of loss rate (p_ell), this plasmid-only state will contain a small fraction of plasmid-free cells due to plasmid loss).Conjugation is the best studied mechanism of plasmid transmission, but plasmids can instead be transmitted by transformation, whereby plasmid-free cells are infected by free-floating plasmids, as illustrated in Fig. 1B. We therefore consider a model for plasmid-spread via transformation in which cell death results in release of free-floating plasmids which can then infect cells by mass action at rate (gamma _{mathrm{t}}). For every cell death, (n_{{mathrm{eff}}}) free-floating plasmids are released and these plasmids decay at a rate (delta _{mathrm{p}}). The dynamics of transformative plasmids are therefore:$$ frac{{drho }}{{dt}} ,=, alpha Crho – gamma _{mathrm{t}}rho P ,+, p_ell (1 ,-, {Delta})alpha Crho _{mathrm{p}} ,-, delta rho ,\ frac{{drho _{mathrm{p}}}}{{dt}} ,=, (1 ,-, {Delta})alpha Crho _{mathrm{p}} ,+, gamma _{mathrm{t}}rho P ,-, p_ell (1 ,-, {Delta})alpha Crho _{mathrm{p}} ,-, delta rho _{mathrm{p}},\ frac{{dC}}{{dt}} ,=, S ,-, alpha Crho ,-, (1 ,-, {Delta})alpha Crho _{mathrm{p}},\ frac{{dP}}{{dt}} ,=, n_{{mathrm{eff}}}delta rho _{mathrm{p}} ,-, gamma _{mathrm{t}}rho P ,-, delta _{mathrm{p}}P.$$
(3)
What is the condition for transformative plasmid invasion? The plasmid-free equilibrium is unstable if$$qquadqquadquadgamma _{mathrm{t}}rho ^ ast , > , delta _{mathrm{p}}left( {frac{{{Delta} ,+, p_ell (1 ,-, {Delta})}}{{n_{{mathrm{eff}}} ,-, {Delta} ,-, p_ell (1 ,-, {Delta})}}} right).$$
(4)
The left-hand-side of Eq. 9 is similar to the conjugative plasmid invasion condition, with the conjugation rate (gamma _{mathrm{c}}) replaced by the transformation rate (gamma _{mathrm{t}}). The numerator of the right-hand-side is also similar, with the cell death rate (delta) replaced with the plasmid decay rate (delta _{mathrm{p}}). The primary difference is in the denominator, which is the difference between the number of plasmids released on cell death, (n_{{mathrm{eff}}}), and the total replication deficit of plasmid-containing cells. If this denominator is negative, the inequality reverses and the plasmid-free equilibrium is always stable.The invasion condition in Eq. 9 determines the optimal (n_{mathrm{p}}) of transformative plasmids: if each plasmid within a cell has a fixed probability of remaining viable after cell death, (p_{mathrm{v}}), then (n_{{mathrm{eff}}}) will scale linearly with (n_{mathrm{p}}) such that (n_{{mathrm{eff}}} ,=, p_{mathrm{v}}n_{mathrm{p}}). If the denominator of Eq. 9 is positive, the optimal copy number will be (n_{mathrm{p}} ,=, 1/{Delta}_{mathrm{p}}), the point at which the host’s growth rate is driven to zero and the plasmid relies entirely on horizontal transfer to survive. These results are substantially different than in the case of conjugation: instead of restricting itself to a limited portion of the host’s metabolic budget, a transformative parasite maximizes its spread by using as much of the host’s resources as possible. This is reminiscent of the behavior of phages—suggesting a possible evolutionary link between parasitic plasmids and phages.As in the conjugation case, we now explore the ecological outcomes possible with transformative plasmids. We again consider the case of negligible loss rate (p_ell ,=, 0) and characterize the stability of the equilibria (see SI Appendix 1 for details). For (n_{{mathrm{eff}}} , > , 1), the system has similar ecological outcomes to the conjugative case, with the system transitioning through no-plasmid, coexistence, and plasmid-only equilibria as ({Delta}) decreases and (gamma _{mathrm{t}}) increases. Interestingly, when (n_{{mathrm{eff}}} , , 0.} end{array}$$
(5)
Fig. 2: Competition between populations may prevent runaway plasmid invasion.A Illustration of multiple populations, each occupying an isolated “deme”. During each epoch, populations compete for demes, with plasmid invasion occurring randomly (see Eq. 11 for details). In the example shown, in the first epoch, the population with two plasmids is replaced by the population with zero plasmids. In the second epoch, the population with magenta plasmids is invaded by the green plasmid. B Multiplasmid fitness costs for different types of epistasis. With no epistasis, fitness burden is multiplicative as in Fig. 1F. With positive epistasis, fitness burden increases sub-multiplicatively (pictured: ({Delta}_{{mathrm{tot}}} ,=, {Delta}) for (m , > , 0)). For negative epistasis, fitness burden increases super-multiplicatively (pictured: ({Delta}_{{mathrm{tot}}} ,=, 1 ,-, (1 ,-, {Delta})^{m^{3/2}})). C Steady-state distributions of number of plasmid types per cell in the Wright–Fisher model (see SI Appendix 3). Parameters ({Delta} ,=, 0.01) and plasmid invasion probability for each time period (q ,=, 0.005).Full size imageA population’s fitness is dependent on the number of unique plasmid types it contains. Thus far, we have considered a simple multiplicative model. However, it has been demonstrated that plasmid–plasmid interactions can modulate plasmid properties. For example, one study found that the presence of a plasmid can reduce the fitness cost of an invading plasmid [12]. To account for this epistasis between plasmids, we also consider fitness costs that increase sub-multiplicatively (positive epistasis) or super-multiplicatively (negative epistasis). We show examples of positive epistasis, negative epistasis, and no epistasis in Fig. 2B.What is the distribution of unique plasmid types across populations in our model with HGT barriers? We derive the stationary distribution of this model for the three different epistasis functions in Fig. 2B and plot them in Fig. 2C (see SI Appendix 3 for details). For the case of no epistasis, the stationary distribution is Poisson-like. Positive epistasis favors carriage of multiple plasmids and results in an exponential-like distribution with a long tail. Negative epistasis has the opposite effect: it penalizes carriage of multiple plasmids and results in a sub-Poissonian distribution with a reduced tail. Importantly, in all cases the runaway invasion of plasmids is stopped. While there is nothing stopping individual populations from being overrun by invading plasmids, these populations are more likely to be out-competed by populations with fewer plasmids. Thus, the single-population “tragedy of the commons” is counteracted at a higher level of selection.Analysis of natural genomesHow does our predicted distribution of unique plasmid types per cell compare to that in natural genomes? To make this comparison, we downloaded all complete bacterial genomes from NCBI (a total of 17,725 genomes) and analyzed their plasmid content. In Fig. 3A, we show the overall distribution of unique plasmid types per genome and corresponding model fits for both positive and no epistasis cases (see “Methods” for fitting details). The natural distribution is exponential-like and is well-fit by a model with positive epistasis. The model fit with no epistasis has too short a tail to be able to fit the data, and this problem becomes even more severe for negative epistasis. Thus, interestingly, we find that the distribution of unique plasmid types in real-world genomes is consistent with parasitic plasmids that ameliorate each other’s fitness costs. The degree of positive epistasis suggested by the data is quite strong—the distribution is nearly a pure exponential. In our model, this corresponds to the case in which the cost of all plasmids beyond the first is zero, such that for (m , > , 1) the parameters controlling both population replication and plasmid invasion are independent of plasmid number. This means that the ratio between consecutive elements of the distribution is constant, yielding an exponential tail. In order to determine whether our conclusions are influenced by oversampling of clinically relevant species, we excluded 91 genera known to be clinically relevant or human-associated and repeated our analysis. The remaining dataset contains nearly 5000 genomes and still shows clear exponential behavior (see SI Fig. S4). We also analyzed whether the presence of engineered strains within the NCBI database influences our results. We found that there are only a small number of these engineered strains and that removing them had negligible impact on our results (see SI Fig. S5).Fig. 3: Comparison of distributions of number of unique plasmid types per cell in natural genomes to Wright–Fisher model.A Distribution of number of plasmid types per cell in 17,725 complete NCBI genomes. Positive epistasis distribution fit with the fitness function ({Delta}_{{mathrm{tot}}} ,=, {Delta}) for (m , > , 0) (best-fit parameters: ({Delta} ,=, 9.8 ,times, 10^{ – 3}), (q ,=, 5.4 ,times, 10^{ – 3})), no epistasis distribution fit with ({Delta}_{{mathrm{tot}}} ,=, 1 ,-, (1 ,-, {Delta})^m) (best-fit parameters: ({Delta} ,=, 3.9 ,times, 10^{ – 3}), (q ,=, 1.4 ,times, 10^{ – 2})). B Distribution of number of plasmid types per cell in 1153 complete Escherichia genomes, with a positive epistasis fit using the fitness function ({Delta}_{{mathrm{tot}}} ,=, 1 ,-, (1 ,-, {Delta})^{m^a}) (best-fit parameters: ({Delta} ,=, 8.3 ,times, 10^{ – 3}), (q ,=, 8 ,times, 10^{ – 3}), (a ,=, 0.33)). C Distribution of number of plasmid types per cell in 576 complete Klebsiella genomes, with a positive epistasis fit using the fitness function as in (B) (best-fit parameters: ({Delta} ,=, 7 ,times, 10^{ – 3}), (q ,=, 9.7 ,times, 10^{ – 3}), (a ,=, 0.43)). Note that in certain limits of our models, only the ratio of (q) and ({Delta}) can be properly estimated, effectively reducing them to single parameter (see SI Appendix 3). D Distribution of number of plasmid types per cell in genomes containing and not containing cas genes. Genomes are considered cas containing if at least one chromosome or plasmid within the genome contains a cas gene. See “Methods” for details.Full size imageCan our model capture variation within smaller, related groups of genomes? In Fig. 3B we show the distribution of unique plasmid types per cell within the genus Escherichia. As can be seen, the data is very well fit by a model of parasitic plasmids with positive epistasis. However, our model was not able to capture some of the within-genus distributions we encountered. A notable exception is the distribution of unique plasmid types per cell in the genus Klebsiella, shown in Fig. 3C. In this genus, there is a substantial discontinuity between the zero-plasmid class and the rest of the distribution. While our simple Wright–Fisher model with some positive epistasis can capture the tail of the distribution, it then fails to capture the first few classes. Despite such exceptions, we find that the positive epistasis model is generally able to capture the overall trends in plasmid distributions over the bulk of natural genomes (see SI Fig. S6).It should be noted that our current model of constant plasmid invasion probability and strong positive epistasis is not the only Wright–Fisher model that can produce an exponential distribution matching the data. We analyzed a more general form of the Wright–Fisher model in which the invasion probability and total fitness cost are arbitrary functions of unique plasmid number (see SI Appendix). We find that the general condition to yield an exponential is that the plasmid invasion probability and total fitness cost must be comparable regardless of the number of plasmids in the cell. These results indicate that even if there is no epistasis in fitness cost, an exponential can still result if there is positive epistasis in the invasion probability (i.e., if existing plasmids make it more likely for a new plasmid to successfully invade).HGT barriers are not the only mechanism that can plausibly limit runaway plasmid invasion. Cells also have specialized systems to defend against foreign DNA, notably the CRISPR-Cas system [32]. To explore whether CRISPR-Cas is responsible for limiting plasmid invasion in natural genomes, we searched for cas genes within the NCBI complete bacterial genomes using HMMER (see “Methods” for details). We expect that if CRISPR-Cas plays a major role in limiting the spread of plasmids, the distribution of unique plasmid types per cell would be shifted towards lower plasmid numbers in genomes containing cas genes versus those lacking cas genes. In Fig. 3D, we show the distribution of unique plasmid types per genome in genomes containing at least one cas gene and those not containing any cas genes. The distributions are very similar, with no large differences between them. These results suggest that CRISPR-Cas is not a major mechanism limiting the spread of plasmids in bacteria. There are additional defense systems that may also influence plasmid carriage. However, a prior bioinformatics study found results similar to ours for restriction-modification (RM) systems, another defense system that protects against foreign DNA; the study examined the distribution of RM systems in bacterial genomes and found almost no relationship between the number of RM systems a genome encodes and the presence of plasmids (in one subset of data the authors actually found a positive relation) [33]. More

Here we present five current shortcomings identified in China’s biodiversity protection framework.Varying threat-assessment quality and uniform treatment of speciesIn this section, we highlight how the threat classifications of the Catalogue of Wildlife under Special State Conservation can lead to sentences that are not commensurate with the species’ threat level. In recent amendments to the catalogue, insect species occur in the highest protection classes (3 species out of 234 in Class I and 72 species out of 746 in Class II; Fig. 2) with similar sentencing standards as for large mammals and birds. For instance, killing more than six individuals of Class I protected insects is treated equally to killing one giant panda, with a punishment of at least ten years’ imprisonment according to the Judicial Interpretation of Several Questions Concerning the Application of Law in the Trial of Criminal Cases of Destruction of Wildlife Resources.Fig. 2: Example species with the highest protection status but considerably different life histories.a,b, Mammals such as the giant panda (a) and insects such as the butterfly T. aureus (b) both occur in the highest protection category in the Catalogue of Wildlife under Special State Conservation. Credit: Juping Zeng (b).Full size imageIn June 2002, 10 poachers captured 263 adults of the butterfly Teinopalpus aureus, meant to be sold on the black market. As T. aureus is listed in Class I of the Catalogue of Wildlife under Special State Conservation, based on the assumption of being rare, the punishment was 5 to 13 years’ imprisonment20. However, recent observations indicate both a wider distribution range21,22 and larger population sizes than initially assumed23. Further, the reproduction rate of insects is generally much higher than that of mammals, which usually makes insects more resistant to the removal of specimens. This case raised some controversy about the scientific basis for classification and the financial profit that can be made with insects compared with mammals24. On the black market, T. aureus can be sold for 700 Chinese yuan per male (~US$100; US$1 = 6.9932 yuan, 21 July 2020; gross domestic product (GDP) per capita: 30,808 yuan in 2010, 54,139 yuan in 2016) and 3,500 yuan per female (~US$500; personal communication with collectors in 2011), while a pair of giant pandas is usually rented to abroad zoos for about 7 million yuan (~US$1 million) per year25.In 2015, a college student and a farmer took 16 fledglings of the Eurasian hobby (Falco subbuteo), a Class II protected species, and were sentenced to 10.5 and 10 years’ imprisonment and fines of 10,000 and 5,000 yuan, respectively26. However, ecological studies indicate that the distribution range, population density and reproduction rate of F. subbuteo in China seem sufficient for sustaining viable populations27, highlighting the potential of overly harsh punishment when classification lacks scientific basis.In contrast to valuation according to (black) market prices, wild species also provide higher-level socioeconomic benefits28. For instance, the value of insect pollination services in China was estimated to be 886.5 billion yuan (US$131 billion) in 201529. In comparison, the ecosystem services related to the giant panda were estimated at between 18 billion and 48 billion yuan per year (US$2.6–6.9 billion) in 2010, but they seem more indirect via regulating, provisioning and cultural services provided by the panda reserves30. However, pollination services are provided by multiple species within a highly flexible network31,32 and the impact of removing a particular amount of specimens is hard to assess, whereas large mammals, such as the giant panda, are irreplaceable in ecosystems and their roles as umbrella species. Thus, differences between insects and mammals are striking not only in terms of direct financial profit but also in terms of ecological and socioeconomic damage, and therefore it is questionable that they are both listed in the highest protection class with the same stringent punishment.Lack of quantitative sentencing standards for herbaceous plants, fungi and algaeHere, we discuss how limited scientific knowledge for particular species groups can lead to legal uncertainties and consequently to limited protection or overly harsh punishment. The Regulations of the People’s Republic of China on the Protection of Wild Plants identify the legal responsibilities for the protection of wild plants (excluding trees), but have not yet reached the status of a law and thus are without judicial interpretation of the Supreme People’s Court and respective sentencing standards. Instead, stipulations of ‘seriousness’ are used with regard to the sentences used for trees, defined in the Judicial Interpretation of Several Questions Concerning the Specific Application of Law in the Trial of Criminal Cases of Destruction of Forest Resources (Box 1), and respective sentencing standards, defined in the Criminal Law of the People’s Republic of China, are applied (up to seven years’ imprisonment). With this analogy, an offender was sentenced to three years in prison in 2016 (suspended sentence) and a fine of 1,000 yuan for digging out three stems of Cymbidium faberi33, an orchid listed in Appendix II of the Convention on International Trade in Endangered Species of Wild Fauna and Flora (CITES34; Fig. 3d) but with high market value. Some uncertainty in the legal position regarding herbaceous plants is expressed by another case in the same year, in which an offender was sentenced to one year of imprisonment (fine of 5,000 yuan) for digging out 55 stems of C. faberi35, and the later revocation of the sentences given that C. faberi is not listed in the Catalogue of Wild Plants under Special State Conservation36.Fig. 3: Example species with changing threat status.a–d, Wildlife protection laws need to be adaptive to reflect the recovery of formerly threatened species, such as the snow leopard (Panthera uncia; a) or the kiang (E. kiang; b), or the increasing endangerment of initially non-threatened species, such as the butterfly Bhutanitis lidderdalii (c) or the orchid C. faberi (d). Credit: Zhi Lu (a, b); Lixin Zhu (c); Yu Ren (d).Full size imageSimilar to the non-discrimination of large mammals and insects, we find such an approach also questionable for precious trees and other plants. Such analogies might become almost impossible when applied to algae such as Nostoc flagelliforme, an important water and soil conservation and high-priced food algae but under Class I protection37. The main reason for the lack of quantitative sentencing standards for these organisms is limited evidence. Therefore, we think it is necessary to raise the Regulations of the People’s Republic of China on the Protection of Wild Plants to become law with respective judicial interpretations and to establish comprehensive scientific assessments targeting herbaceous plants, fungi and algae to provide a solid basis for the development of sentencing standards.Lack of legislative flexibility to reflect dynamic changes in status and taxonomyWe identified a lack of regular updates of the Catalogues of Wildlife and Wild Plants under Special State Conservation needed to address the dynamic changes in taxonomy and threat status. Since its promulgation, the Wildlife Protection Law of the People’s Republic of China has been revised four times and the Regulations of the People’s Republic of China on the Protection of Wild Plants was amended once in 200138, but the Catalogues of Wildlife and Wild Plants under Special State Conservation have basically remained unchanged for the past 32 and 20 years, respectively, with the exception of a recent amendment of the Catalogue of Wildlife in February 2021 and a pending amendment of the Catalogue of Wild Plants (Box 1). Taxonomies change dynamically, which can lead to considerable incongruences among scientifically accepted species names and those in the respective protection lists39. Until this recent amendment, there has been a mismatch in the names of 25 threatened species as listed under CITES compared with the Catalogue of Wildlife under Special State Conservation, putting them at particular risk because their protection status might be questioned, for example, when species such as the Himalayan goral (Naemorhedus goral), or even genera such as the leaf monkeys (Presbytis spp.), have been split into different units with different names that are not listed in the respective catalogues40. Although the Catalogue of Wildlife under Special State Conservation has been updated very recently, it is still recommended that such updates are done regularly and in a coordinated manner, not only in China but across all CITES signatory nations40.Additional legislative flexibility is also needed when formerly endangered species have recovered11, while others have become endangered16,41 (Fig. 3). Recently, several mammals such as the giant panda, snow leopards or the kiang (Equus kiang)11,42 have considerably recovered and their threat status has been reduced by the International Union for Conservation of Nature (IUCN)11. Although the Chinese government does not follow such a downgrade because of precautionary reasons, we think that the sentencing threshold for such species should be adapted in the Judicial Interpretation of Several Questions Concerning the Application of Law in the Trial of Criminal Cases of Destruction of Wildlife Resources. On the other hand, species whose endangerment has increased since the promulgation of the Catalogues of Wildlife and Wild Plants under Special State Conservation, such as the narrow-ridged finless porpoise43, many birds44, snakes45, turtles46, frogs40, butterflies47 or herbaceous (medicinal) plants2, have long been with low or no protection until the recent amendment. Cultivation can also increase endangerment of wild species by hybridization between the cultivars and the wild populations (for example, rice, wheat, soybean and cotton)48.Outdated punishment standards based on economic profitsSimilar to the lack of flexibility covering species’ taxonomic and threat status, here we highlight that punishment standards are outdated and regular updates are required to reflect economic developments and guarantee balanced sentencing. For instance, according to the Judicial Interpretation of Several Questions Concerning the Application of Law in the Trial of Criminal Cases of Destruction of Wildlife Resources, the illegal purchase, transport and sale of precious and endangered wildlife products will be considered as a ‘serious crime’ if the financial profit is more than 100,000 yuan and as an ‘extremely serious crime’ if the profit is 200,000 yuan or more. The sentencing standard was developed in the year 2000, but with the rapid development of China’s economy, nationwide per capita income has increased more than fourfold from 6,279 yuan in 2000 to 28,228 yuan in 201849. To reflect economic developments, the penalty standards need to be adjusted to comply with the principle of balanced sentencing. In comparison, the Chinese standards for corruption and bribery have been increased from 4,886 yuan in 1997 to currently 30,715 yuan for crimes involving a ‘relatively large amount’, which might serve as a guideline for adapting the sentencing standards for wildlife protection50.Potential for excessive punishment because of non-discrimination between organized and individual wildlife crimeIn this section, we highlight that ignoring the motivational, educational and economic backgrounds of offenders is against the principle of proportionality and may lead to inappropriate deterrence strategies. China’s laws are very strict with quite harsh penalty sentencing; for example, 10.5 years’ imprisonment and a fine of 10,000 yuan for a student taking birds26, 12 years and a fine of 10,000 yuan for a farmer killing a giant panda51 or 13 years and a fine of 2,000 yuan for a farmer taking butterflies20, all cases representing ‘extremely serious crimes’ with a minimum sentencing standard of 10 years’ imprisonment (no maximum defined). Even in comparison with other criminal fields in China and internationally, these standards seem very stringent. For instance, sentences of more than 10 years’ imprisonment apply to larceny only if the value of the stolen goods is larger than 500,000 yuan, or to the theft of first-class cultural relics (all valued in the millions; Criminal Law of the People’s Republic of China, Article 264). Also in comparison, the United Nations Convention Against Transnational Organized Crime52 defines much lower sentencing standards, with at least four years’ imprisonment for a ‘serious crime’. In contrast to China, the wildlife protection laws of Western and many other developing countries prioritize monetary fines over imprisonment. Under European wildlife law53, for example, hunting or destroying Class I protected species is generally punishable by a fine and will be sentenced with fixed-term imprisonment only if the case is ‘extremely serious’. In the United States, the maximum imprisonment is a year, with fines of up to US$50,000 (340,000 yuan)54; in the UK, 6 months and fines of up to £20,000 (177,000 yuan)55,56; in India, 3–7 years and a minimum fine of 25,000 rupees (2,300 yuan)57; or in Brazil, 3 months to a year plus fines58.The wildlife protection laws of such countries may provide useful examples for China, but to adhere to the principle of proportionality, motivational, educational and economic backgrounds, in particular a differentiation between organized wildlife crimes and individual violations needs to be considered. Individual and organized crimes are currently not differentiated in the Criminal Law of the People’s Republic of China. Historically, wildlife crime was considered a local activity performed by single individuals. However, at present criminal networks are highly involved59 and resulting economic damage from environmental crime has been estimated to range between US$91 billion and US$259 billion globally60, with the profits of illegal wildlife trade ranging between US$7 billion and US$23 billion61, which is of similar orders to human trafficking, and arms and drug dealing62. In China, the consumption of illegal wildlife products has increased with growing economic wealth63, while China has also been identified as one of the major exporters of such products64. Key players in both cases are organized crime groups65,66, causing severe ecological damage while making enormous financial profits67. In such cases, high fines might be simply factored in as part of the ‘business model’. Thus, the current focus on severe jail sentences seems appropriate, and the level is comparable to other Southeast Asian countries (Indonesia: 10 years; Singapore: 2 years; Thailand: 7 years; Vietnam: 15 years)68,69.In contrast to organized wildlife crime, we also noticed that many cases of harvesting or poaching protected wildlife happened in remote and less-developed regions, conducted by individuals seeking to earn some extra income but without good knowledge of the protection laws20,51. The resulting ecological damage and profits gained are much lower compared with cases of organized wildlife crime, and thus applying the same harsh punishments, as shown in our earlier examples, is clearly against the principle of proportionality. Moreover, it has been shown that the mentality of different types of offender and how they perceive different punishments (imprisonment, fines or both) need to be considered for designing appropriate deterrence strategies for different offence categories, suggesting that imprisonment as the main policy instrument is inappropriate70. Imprisonment is not necessarily a deterrent for every offender, especially when the price of time in prison falls relative to the price of time outside71. Consequently, a penalty that eliminates any financial gain should eliminate the incentive to engage in such conduct72. A shift in focus from imprisonment to fines, at best coupled with local or regional GDP per capita and in combination with raising public awareness, might not only increase proportionality and effectiveness of environmental laws but also comply with other international standards, where, for example, the Council of Europe’s Recommendation (92)17, concerning consistency in sentencing, paragraph B5(2), states that “custodial sentences should be regarded as a sanction of last resort, and should therefore be imposed only in cases where, taking due account of other relevant circumstances, the seriousness of the offence would make any other sentence clearly inadequate”. More

1.Tapia, J., Audry, S., Townley, B. & Duprey, J. L. Geochemical background, baseline and origin of contaminants from sediments in the mining-impacted Altiplano and Eastern Cordillera of Oruro, Bolivia. Geochemistry 12, 3–20. https://doi.org/10.1144/1467-7873/10-RA-049 (2012).CAS 
Article 

Google Scholar 
2.Sarricolea, E. P. & Romero, H. Variabilidad y cambios climáticos observados y esperados en el Altiplano del norte de Chile. Revista de Geografía Norte Grande. 62, 169–183 (2015).Article 

Google Scholar 
3.Garreaud, R., Vuille, M. & Clement, C. A. The climate of the Altiplano: observed current conditions and mechanisms of past changes. Palaeogeogr. Palaeoclimatol. Palaeoecol. 194(1–3), 5–22. https://doi.org/10.1016/S0031-0182(03)00269-4 (2003).Article 

Google Scholar 
4.Vuille, M. & Keiming, F. Interannual variability of summertime convective cloudiness and precipitation in the central andes derived from ISCCP-B3 data. J. Clim. 17(17), 3334–3348. https://doi.org/10.1175/15200442(2004)017%3c3334:IVOSCC%3e2.0.CO;2 (2004).ADS 
Article 

Google Scholar 
5.Cerveny, R. Present climates of South America. In Climates of the Southern Continents: Present, Past and Future (ed. Hobbs, J. E.) 107–135 (Wiley, Chichester, 1998).
Google Scholar 
6.Garreaud, R. & Aceituno, P. Interannual rainfall variability over the South American Altiplano. J. Clim. 14(12), 2779–2789. https://doi.org/10.1175/1520-0442(2001)014%3c2779:IRVOTS%3e2.0.CO;2 (2001).ADS 
Article 

Google Scholar 
7.Coronel, J., Declerck, S. & Brendonck, L. High-altitude peatland temporary pools in Bolivia house a high cladoceran diversity. Wetlands 27(4), 1166–1174. https://doi.org/10.1672/0277-5212(2007)27[1166:HPTPIB]2.0.CO;2 (2007).Article 

Google Scholar 
8.Dorador, C., Vila, I., Witzel, K. P. & Imhoff, J. F. Bacterial and archaeal diversity in high altitude wetlands of the Chilean Altiplano. Fundam. Appl. Limnol. 182(2), 135–159. https://doi.org/10.1127/1863-9135/2013/0393 (2013).CAS 
Article 

Google Scholar 
9.Garcia, E. & Otto, M. Caracterización ecohidrológica de humedales alto andinos usando imágenes de satélite multitemporales en la cabecera de cuenca del río Santa, Ancash, Perú. Ecología Aplicada 14(2), 115–125 (2013).
Google Scholar 
10.Buytaert, W., Camacho, F. C. & Tobón, C. Potential impacts of climate change on the environmental services of humid tropical alpine regions. Glob. Ecol. Biogeogr. 20(1), 19–33. https://doi.org/10.1111/j.1466-8238.2010.00585.x (2011).Article 

Google Scholar 
11.Hribljan, J. A. et al. Carbon storage and long-term rate of accumulation in high-altitude Andean peatlands of Bolivia. Mires Peat 15(12), 1–14 (2015).
Google Scholar 
12.Yager, K. et al. Dimensiones socioecológicas del cambio del paisaje pastoral andino: puente entre el conocimiento ecológico tradicional y el análisis de imágenes satelitales en Sajama Parque Nacional, Bolivia. Cambio ambiental regional 17, 27–37 (2019).
Google Scholar 
13.Urrutia, R. & Vuille, M. Climate change projections for the tropical Andes using a regional climate model: temperature and precipitation simulations for the end of the 21st century. J. Geophys. Res. 114, D02108. https://doi.org/10.1029/2008JD011021 (2009).ADS 
Article 

Google Scholar 
14.Buytaert, W. et al. Uncertainties in climate change projections and regional downscaling in the tropical Andes: implications for water resources management. Hydrol. Earth. Syst. Sci. 14, 1247–1258. https://doi.org/10.5194/hess-14-1247-2010 (2010).ADS 
Article 

Google Scholar 
15.Rabatel, A. et al. Current state of glaciers in the tropical Andes: a multi-century perspective on glacier evolution and climate change. Cryosphere 7, 81–102. https://doi.org/10.5194/tc-7-81-2013 (2013).ADS 
Article 

Google Scholar 
16.Prieto, G. et al. A mass sacrifice of children and camelids at the Huanchaquito-Las Llamas site, Moche Valley, Peru . PLoS ONE 14(3), e0211691. https://doi.org/10.1371/journal.pone.0211691 (2019).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
17.Babidge, S., Kalazich, F., Prieto, M. & Yager, K. That’s the problem with that lake; it changes sides’: mapping extraction and ecological exhaustion in the Atacama. J. Political Ecol. 26(1), 738–760. https://doi.org/10.2458/v26i1.23169 (2019).Article 

Google Scholar 
18.Prieto, M., Fragkou, M. & Calderón, M. Water policy and management in Chile. In The Wiley Encyclopedia of Water: Science, Technology, and Society (ed. Strickland, C.) 2589–2600 (Wiley, New York, 2020).
Google Scholar 
19.Fritz, S. C., Baker, P. A., Tapia, P., Spanbauer, T. & Westover, K. Evolution of the Lake Titicaca basin and its diatom flora over the last 370,000 years. Palaeogeogr. Palaeoclim. Palaeoecol. 317–318, 93–103 (2012).ADS 
Article 

Google Scholar 
20.Cohen, S. C. Scientific drilling and biological evolution in ancient lakes: lessons learned and recommendations for the future. Hydrobiologia 682(1), 3–25. https://doi.org/10.1007/s10750-010-0546-7 (2012).CAS 
Article 

Google Scholar 
21.Tapia, P. M., Fritz, S. C., Baker, P. A., Seltzer, G. A. & Dunbar, R. B. A Late Quaternary diatom record of tropical climatic history from Lake Titicaca (Peru and Bolivia). Palaeogeogr. Palaeoclimatol. Palaeoecol. 194, 1–3. https://doi.org/10.1016/S0031-0182(03)00275-X (2003).Article 

Google Scholar 
22.Vining, B. R., Steinman, B. A., Abbott, M. B. & Woods, A. Paleoclimatic and archaeological evidence from Lake Suches for highland Andean refugia during the arid middle-Holocene. The Holocene 29(2), 328–344. https://doi.org/10.1177/0959683618810405 (2019).ADS 
Article 

Google Scholar 
23.Fritz, S. C., Baker, P. A., Tapia, P. & Garland, J. Spatial and temporal variation in cores from Lake Titicaca, Bolivia/Peru during the last 13,000 years. Quat. Int. 158(1), 23–29. https://doi.org/10.1016/j.quaint.2006.05.014 (2006).Article 

Google Scholar 
24.Hernández, A. et al. Biogeochemical processes controlling oxygen and carbon isotopes of diatom silica in Late Glacial to Holocene lacustrine rhythmites. Palaeogeogr. Palaeoclimatol. Palaeoecol. 299(3–4), 413–425. https://doi.org/10.1016/j.palaeo.2010.11.020 (2012).Article 

Google Scholar 
25.Placzek, C. et al. Climate in the dry central Andes over Geologic, millennial, and interannual timescales. Ann. Mo. Bot. Gard. 96(3), 386–397. https://doi.org/10.3417/2008019 (2009).Article 

Google Scholar 
26.Cerda, M. et al. A new 20th century lake sedimentary record from the Atacama Desert/Chile reveals persistent PDO (Pacific Decadal Oscillation) impact. J. S. Am. Earth Sci. 95, 102302. https://doi.org/10.1016/j.jsames.2019.102302 (2019).Article 

Google Scholar 
27.Aránguiz-Acuña, A. et al. Aquatic community structure as sentinel of recent environmental changes unraveled from lake sedimentary records from the Atacama Desert, Chile . PLoS ONE 15(2), e0229453. https://doi.org/10.1371/journal.pone.0229453 (2020).CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
28.Flores-Varas, A. et al. Ascotán and Carcote salt flats as sensors of humidity fluctuations and anthropic impacts in the transition zone of the Andean Altiplano. J. S. Am. Earth Sci. 105, 102934. https://doi.org/10.1016/j.jsames.2020.102934 (2021).CAS 
Article 

Google Scholar 
29.Maher, B. A. & Taylor, R. M. Formation of ultrafine-grained magnetite in soils. Nature 336, 368–370 (1988).ADS 
CAS 
Article 

Google Scholar 
30.Dearing, J. A. et al. Frequency-dependent susceptibility measurements of environmental materials. Geophys. J. Int. 124, 228–240. https://doi.org/10.1111/j.1365-246X.1996.tb06366.x (1996).ADS 
Article 

Google Scholar 
31.Evans, M. & Heller, F. Environmental magnetism: principles and applications of enviromagnetics. Int. Geophys. 86, 202 (2003).
Google Scholar 
32.Pizarro, H. et al. The origin of the magnetic record in Eocene-Miocene coarse-grained sediments deposited in hyper-arid/arid conditions: examples from the Desert. Palaeogeogr. Palaeoclimatol. Palaeoecol. 516, 322–335. https://doi.org/10.1016/j.palaeo.2018.12.009 (2019).Article 

Google Scholar 
33.Risacher, F., Alonso, H. & Salazar, C. Geoquímica de Aguas en Cuencas cerradas: I, II y III Regiones-Chile. Volumen III Estudio de Cuencas de la II Región. (Convenio de Cooperación DGA-UCN-IRD. S.I.T. Nº51, 1999).34.Tapia, R. et al. Glacial differences of Southern Ocean Intermediate Waters in the Central South Pacific. Quat. Sci. Rev. 208, 105–117. https://doi.org/10.1016/j.quascirev.2019.01.016 (2019).ADS 
Article 

Google Scholar 
35.Horne, D. J. Life-cycles of podocopid Ostracoda – a review (with particular reference to marine and brackish-water species). In Applications of Ostracoda. Proceedings of the Eighth International Symposium on Ostracoda (ed. Maddocks, R.) 581–590 (University of Houston, Texas, 1983).
Google Scholar 
36.Cohen, A. C. & Morin, J. G. Patterns of reproduction in ostracodes; a review. J. Crust. Biol. 10(2), 184–211. https://doi.org/10.2307/1548480 (1990).Article 

Google Scholar 
37.Mesquita-Joanes, F., Smith, A. J. & Viehberg, F. A. The ecology of ostracoda across levels of biological organization from individual to ecosystem. J. Quat. Sci. 17, 15–35. https://doi.org/10.1016/B978-0-444-53636-5.00002-0 (2012).Article 

Google Scholar 
38.McLay, C. L. The population biology of Cyprinotus carolinensis and Herpetocypris reptans (Crustacea, Ostracoda). Can. J. Zool. 56(5), 1170–1179. https://doi.org/10.1139/z78-161 (1978).ADS 
Article 

Google Scholar 
39.Hamouda, S. A., Sames, B., Mohammed, A. & Bensalah, M. First record of non-marine ostracods from the Paleogene “hamadian deposits” of Méridja area, west of Bechar (southwestern Algeria). Annales de Paléontologie 104(1), 27–44. https://doi.org/10.1016/j.annpal.2017.12.001 (2018).Article 

Google Scholar 
40.Bergue, C. T., Maranhao, M. D. S. A. S. & Fauht, G. Paleolimnological inferences based on Oligocene ostracods (Crustacea: Ostracoda) from Tremembé Formation. Southeast Brazil. An. Acad. Bras. Cienc. 87(3), 1531–1544. https://doi.org/10.1590/0001-3765201520140366 (2015).Article 
PubMed 

Google Scholar 
41.Sylvestre, F., Servant-Vildary, S. & Roux, M. Diatom-based ionic concentration and salinity models from the south Bolivian Altiplano (15–23°S). J. Paleolimnol. 25, 279–295 (2001).ADS 
Article 

Google Scholar 
42.Nunnery, J. A., Fritz, S. C., Baker, P. A. & Selenbien, W. Lake-level variability in Salar de Coipasa, Bolivia during the past ∼40,000 yr. Quat. Res. https://doi.org/10.1017/qua.2018.108 (2018).Article 

Google Scholar 
43.Herrera, C. et al. Investigaciones hidrogeológicas en la laguna Tuyajto perteneciente a la Reserva Nacional de los Flamencos (Atacama, Chile). Bol. Geol. Min. 130(4), 789–806. https://doi.org/10.21701/bolgeomin.130.4.011 (2019).Article 

Google Scholar 
44.Houston, J. Variability of precipitation in the Atacama Desert: its causes and hydrological impact. Q. J. R. Meteorol. Soc. 26(15), 2181–2198. https://doi.org/10.1002/joc.1359 (2006).Article 

Google Scholar 
45.Herrera, C. et al. Groundwater flow in a closed basin with a saline shallow lake in a volcanic area: Laguna Tuyajto, northern Chilean Altiplano of the Andes. Sci. Total Environ. 541, 303–318. https://doi.org/10.1016/j.scitotenv.2015.09.060 (2016).ADS 
CAS 
Article 
PubMed 

Google Scholar 
46.Munk, L. A., Boutt, D. F., Hynek, S. A. & Moran, B. J. Hydrogeochemical fluxes and processes contributing to the formation of lithium-enriched brines in a hyper-arid continental basin. Chem. Geol. 493, 37–57. https://doi.org/10.1016/j.chemgeo.2018.05.013 (2018).ADS 
CAS 
Article 

Google Scholar 
47.Godfrey, L. & Álvarez-Amado, F. Volcanic and Saline Lithium Inputs to the Salar de Atacama. Minerals 10(2), 201. https://doi.org/10.3390/min10020201 (2020).CAS 
Article 

Google Scholar 
48.Marazuela, M. A., Ayora, C., Vázquez-Suñé, E., Olivella, S. & García-Gil, A. Hydrogeological constraints for the genesis of the extreme lithium enrichment in the Salar de Atacama (NE Chile): A thermohaline flow modelling approach. Sci. Total Environ. 739, 139959. https://doi.org/10.1016/j.scitotenv.2020.139959 (2020).ADS 
CAS 
Article 
PubMed 

Google Scholar 
49.Bobst, A. L. et al. A 106 ka paleoclimate record from drill core of the Salar de Atacama, northern Chile. Palaeogeogr. Palaeoclimatol. Palaeoecol. 173(1–2), 21–42. https://doi.org/10.1016/S0031-0182(01)00308-X (2001).Article 

Google Scholar 
50.Baspineiro, C. F., Franco, J. & Flexer, V. Potential water recovery during lithium mining from high salinity brines. Sci. Total Environ. 720, 137523. https://doi.org/10.1016/j.scitotenv.2020.137523 (2020).ADS 
CAS 
Article 
PubMed 

Google Scholar 
51.Marazuela, M. A., Vázquez-Suñé, E., Ayora, C. & García-Gil, A. Towards more sustainable brine extraction in salt flats: Learning from the Salar de Atacama. Sci. Total Environ. 703, 135605. https://doi.org/10.1016/j.scitotenv.2019.135605 (2020).ADS 
CAS 
Article 
PubMed 

Google Scholar 
52.Babidge, S. Sustaining ignorance: the uncertainties of groundwater and its extraction in the Salar de Atacama, northern Chile. J. R. Anthropol. Inst. 25(1), 83–102. https://doi.org/10.1111/1467-9655.12965 (2018).Article 

Google Scholar 
53.Sonter, L. J., Dade, M. C., Watson, J. E. M. & Valenta, R. K. Renewable energy production will exacerbate mining threats to biodiversity. Nat. Commun. 11, 4174. https://doi.org/10.1038/s41467-020-17928-5 (2020).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
54.Stahl, A. T., Fremier, A. K. & Cosens, B. A. Mapping legal authority for terrestrial conservation corridors along streams. Conserv. Biol. 34(4), 943–955. https://doi.org/10.1111/cobi.13484 (2020).Article 
PubMed 
PubMed Central 

Google Scholar 
55.García, M., Prieto, M. & Kalazich, F. The protection of the mountain ecosystems of the Southern Central Andes: tensions between Aymara herding practices and conservation policies. Eco. Mont. 13(1), 22–30 (2021).
Google Scholar 
56.Vila, T. Geología de los depósitos salinos andinos, provincia de Antofagasta, Chile. Revista de Geología de Chile 2, 41–55 (1975).
Google Scholar 
57.CIREN. Catastro Agrícola, estudio de metodología para la realización y actualización de catastro agrícola regional en base a la utilización de tecnología geoespacial 1–35 (CCIRA, Atacama, 2013).
Google Scholar 
58.Villagrán, C., Kalin-Arroyo, M. T. & Marticorena, C. Efectos de la destización en la distribución de la flora andina de Chile. Rev. Chil. Hist. Nat. 56, 137–157 (1983).
Google Scholar 
59.CONAF. Actualización Plan de Manejo Participativo Reserva Nacional Los Flamencos, Region de Antofagasta (2008).60.Núñez, L., Grosjean, M. & Cartajena, I. Ocupaciones humanas y paleoambientes en la Puna de Atacama (Universidad Católica del Norte-Taraxacum, Antofagasta, 2005).
Google Scholar 
61.Los Ostracodos, M. P. VI. 4f. In El lago Titicaca, síntesis del conocimiento limnológico actual (eds Dejoux, C. & Iltis, A.) 345–352 (Orstom, New Caledonia, 1991).
Google Scholar 
62.Karanovic, I. Recent Freshwater Ostracods of the World, Crustacea, Ostracoda, Podocopida (Springer , Berlin, 2012).
Google Scholar 
63.Palacios-Fest, M. R., Cusminsky, G. C. & McGlue, M. M. Late Quaternary lacustrine ostracods (Ostracoda, Crustacea) and charophytes (Charophyta, Charales) from the Puna Plateau, Argentina. Micropaleontology 35, 66–78 (2016).
Google Scholar 
64.Brandão, S. N., Angel, M. V., Karanovic, I., Perrier, V. & Meidla, T. World Ostracoda Database. http://www.marinespecies.org/ostracoda/aphia.php?p=taxdetails&id=1091 on 2020–01–15 (2018).65.Fatela, F. & Taborda, R. Confidence limits of species proportions in microfossil assemblages. Mar. Micropaleontol. 45(2), 169–174. https://doi.org/10.1016/S0377-8398(02)00021-X (2002).ADS 
Article 

Google Scholar 
66.Díaz, C. P. & Maidana, N. I. Diatomeas de los Salares Atacama y Punta Negra, II Región-Chile (Centro de Ecología Aplicada , La Reina, 2005).
Google Scholar 
67.Diatoms of North America. The source for diatom identification and ecology. https://diatoms.org (2019).68.Hammer, Ø. & Harper, D. Paleontological Data Analysis (Blackwell Publishing , Hoboken, 2006).
Google Scholar 
69.Oksanen, J. et al. vegan: Community Ecology Package. Ordination methods, diversity analysis and other functions for community and vegetation ecologists. Version 2.5–1. URL https://CRAN.R-project.org/package=vegan (2019).70.McArdle, B. H. & Anderson, M. J. Fitting multivariate models to community data: a comment on distance-based redundancy analysis. Ecology 82(1), 290–297. https://doi.org/10.1890/0012-9658(2001)082[0290:FMMTCD]2.0.CO;2 (2001).Article 

Google Scholar  More

Study areaThe study was conducted along the Agi-gawa River, a tributary of the Kiso-gawa River system in central Japan (35°23 42″–35°26 49″N, 137°25 12″–137°28 01″E; Fig. 1), with the Agi-gawa Dam (110 km from the river mouth, 418 m a.s.l.). The Agi-gawa River is a 3rd to 4th-order river with a naturally sand-rich bed derived from weathered granite that characterizes the local geology36. The Agi-gawa Dam (35°25 32″N, 137°25 55″E) had begun operations in 1990; it is a 102 m high rockfill dam with a catchment area of 82 km2, a storage capacity of 4.8 × 107 m3, a mean depth of ~ 45 to 50 m at the dam site, and a hydraulic residence time of 71 days. Although three small sub-dams at the upstream end of the impoundment trap particulates, the sediment speed in the reservoir has been 1,000,000 m3 for 24 years. The dam serves multiple purposes, including flood control, industrial and urban water supply, and the maintenance of baseflow. Further information on the Agi-gawa Dam is available in Katano et al.37.Figure 1The study area shows six study reaches in three stream segments along the Agi-gawa River and Iinuma-gawa Stream, Gifu Prefecture, Japan. Gray circles denote reaches, which are numbered from upstream to downstream within each segment: UD1 and UD2 are upstream of the dam, DD1 and DD2 are downstream of the dam, whilst TR1 and TR2 are in the tributary. The two black circles denote the sediment replenished reaches (S1 and S2). The three small rectangles at the upstream ends of the impoundment are sub-dams, constructed to reduce the inputs of particulates to the impoundment. This map is based on the Digital Topographic Map 25,000 published by Geospatial Information Authority of Japan.Full size imageSediment replenishment and sampling sitesSediment replenishment was undertaken 0.8 and 1.8 km downstream of the Agi-gawa Dam (S1 and S2, Fig. 1) on February 16 and 27, 2005. A total of 1,200 m3 of sediment (D50 ≈ 0.6 mm; mainly sand) was mined from the upstream sub-dams and transported to S1 and S2. We estimate that this constituted 0.086% of the annual sedimentation in the Agi-gawa Dam (e.g., in 2007, replenished sediment per year × 100/sedimentation in the reservoir). The sediment (800 and 400 m3) was replenished at high-flow banks in both sites. The replenished sediment was gradually washed during the high flows at the end of June (visual observation by dam administrators) (Fig. 2). We confirmed that this replenished sediment remained on both banks in March, and no sediments remained on both banks in early July.Figure 2(a) Precipitation (mm·d) (b) mean inflow to the impoundment per day (m3·s-1); and (c) mean outflow from the Agi-gawa Dam per day (m3·s-1). The vertical broken line indicates the study period. Note that the y-axes for (b) and (c) have a logarithmic scale.Full size imageField sampling was conducted twice between March 15 and 18, 2005, prior to sediment flushing and between August 22 and 24, 2005, following sediment flushing [7 weeks after the end of the sediment drift out (Fig. 2)]. The later sampling date was scheduled to investigate the continuous effects (i.e., not immediate effects) of replenished sediment on the riverbed environment and macroinvertebrate assemblages before the replenished sediment had completely been transported further downstream from S1 and S2.Three study segments (length: 1–2 km each) were selected: (1) upstream of the dam and impounded area (UD); (2) downstream of the dam (DD); and (3) in the tributary (TR). These sites were along a 6.0 km stretch of the Agi-gawa River and a 1.0 km stretch of the Iinuma-gawa Stream (catchment area = 24 km2); the latter is a tributary that flows into the Agi-gawa River 2.7 km, downstream of the dam (Fig. 1, Table 1). Each segment contained two study reaches (six reaches in total), and each study reach was 160 m long with several pool–riffle sequences; all reaches were > 300 m apart. DD1 and DD2 were located immediately downstream of the sediment-displaced banks (S1 and S2; 100 m upstream of DD1 and DD2, respectively). Measurements at the two reaches within the same segment were completed on the same day, and the reaches were surveyed in an upstream direction. The dominant land use along the study area was paddy fields, with sparse riparian forest.Table 1 General characteristics of the three study segments and two seasons.Full size tableAlthough the most suitable reference site for DD is the DD prior to dam construction, we were unable to investigate the site prior to dam construction. Therefore, we treated the reference sites as sites that were less affected by the dam than DD on the present day. Katano et al.37 indicated that the difference between the TR and UD sites was smaller than that between DD and UD/TR sites in terms of biota and geology. However, UD was characterized by a wider channel and higher discharge than TR, due to differences in their catchment areas (Table 1). As we did not have a definitive reference, we treated both UD and TR as reference sites (see “Statistical analysis” section). Therefore, how DD in March and DD in August is different from UD and TR can be interpreted as the effect of sediment reduction.Physical environment and water qualitySix riffles were selected at each study reach, and a sampling location (50 × 50 cm quadrat) was established in the mid-channel area of each riffle. Prior to invertebrate sampling, physical environmental variables were measured.Substrate coarseness was measured by gently floating a Plexiglas observation box (50 × 50 × 10 cm deep) divided into four grid squares (25 × 25 cm) on the surface water such that the grid had projected onto the streambed. The size of the substrate material was coded based on the intermediate-axis length: 1 = sand (particles  16 mm) and sieved through a 0.25 mm mesh sieve. Sieved samples and substrate material smaller than pebbles were mixed in a container and preserved in 5% formalin in the field.The material in each container was later divided into two size fractions using 1-and 0.25 mm mesh sieves. To simplify the sorting process, all material retained in the 0.25 mm sieve was mixed and divided into 2n subsamples (maximum n = 32) using a splitter (Idea Co., Tokyo, Japan), following the method described by Vinson and Hawkins43. All macroinvertebrates in subsamples in the 1 mm sieve were counted and identified to the lowest taxonomic level possible, usually to genus or species level using the taxonomic keys of Kawamura and Ueno44, Merritt and Cummins34, Kathman and Brinkhurst45, Kawai and Tanida35, and Torii46.Macroinvertebrate taxa were also classified into five functional feeding groups (FFGs) according to Kawamura and Ueno44, Merritt and Cummins34, Kathman and Brinkhurst45, Kawai and Tanida (2005)35, and Torii46. FFGs were defined as collector-filterers, collector-gatherers, predators, scrapers, and shredders. If a species belonged to ≥ 2 FFGs, the number of individuals was apportioned across the FFGs. We also counted the number of burrowers (#burrowers), inorganic case-bearing caddisflies (#ICB), and net-spinners (#net spinners) of macroinvertebrate assemblages according Kawamura and Ueno44, Merritt and Cummins34, Kathman and Brinkhurst45, Kawai and Tanida35, and Torii46 (see Supplementary Table S1). This classification was carried out as such life-habit traits are important for surviving in a regulated river containing reduced quantities of sand and gravel on the riverbed37. The Chironomidae family was excluded in the life-habit analysis as they consist of various life forms. Once all invertebrates were removed, dry mass (mg m−2) and ash-free dry mass (AFDM, mg m−2) of benthic coarse particulate organic matter (BCPOM,  > 1 mm), and benthic fine particulate organic matter (BFPOM,  0.25 mm) were obtained by drying in an oven at 60 °C for 1 day and combusting in a muffle furnace at 550 °C for 4 h. BCPOM and BFPOM were calculated based on the difference between the dry mass and the AFDM.The total number of invertebrate individuals and the AFDM of BFPOM in each sample were estimated by multiplying by the corresponding 2n value. The number of taxa and density of invertebrates in each sample were calculated as the sum of the values in both size fractions. Additionally, we determined Shannon’s diversity index (H), Simpson’s evenness index, and the percentage of Ephemeroptera, Plecoptera, and Trichoptera (%EPT)47. A sample from UD2 in March had been lost and therefore could not be included in the analyses.Periphyton was sampled from cobbles adjacent to each sampling location. Periphyton was removed from a 5 × 5 cm area on the upper surface of each cobble with a toothbrush. Each sample was placed in a separate container with 250 mL of water. Within 24 h of sample collection, a subsample of the well-mixed content in each container was filtered using a glass-fiber filter (GF/C; Whatman Co., Maidstone, UK). Each filter was placed in a separate vial with 20 mL of 99.5% ethanol and stored in a dark refrigerator at 4 °C for 24 h. The extracted pigments were measured using a spectrophotometer (U-1800; Shimadzu Co., Kyoto, Japan), following the method of Lorenzen48.Analysis of case materials of an inorganic case-bearing caddisflyWe compared the particle size structure of replenished sediment, riverbed sediment, and case materials for case-bearing caddisfly. The replenished sediment was directly sampled in a 1 L polyethylene jar at the upstream replenished bank (S1) on March 16, 2005 (Fig. 1). Riverbed sediment was sampled at two stations; 100 m upstream of S1, and 100 m upstream of DD1 between August 22 and 24, 2005. At each station of the river, a metallic narrow cup (200 mL) with a lid was pushed into a vacancy between the cobbles, which had been randomly selected, and fine sediments (up to small gravel) in the vacancies were sampled by closing the lid underwater. Sampling was carried out three times (i.e., three different vacancies in the cobbles), and subsamples were pooled for measurement. The replenished and riverbed sediment was combusted at 550 °C for 2 h in a muffle furnace to remove organic contamination. Combusted samples were separated with eight sieves with a mesh size range of 0.075–9.5 mm (JIS A 1204). Each fraction was weighed, and the grain size accumulation curve of each type of sediment and its D50 were obtained.In a macroinvertebrate sample at DD1 between August 22 and 24, 2005, ten individuals from two case-bearing caddisfly larvae, Glossosoma sp. and Gumaga orientalis, which were prevalent at DD1 during this period (see Results), were randomly selected from the formalin-fixed sample. The case was carefully removed from the larvae and combusted as described above for the replenished and riverbed sediment. The number of case material grains was measured using a dissection microscope.Statistical analysesWe described results based on two main assumptions: (1) the DD in March is the dam-affected reach (cf. unregulated reaches UD and TR), and (2) the changes in DD from March to August were mainly a result of sediment replenishment. In the statistical analyses, the p criterion (⍺) was set at 0.05.To consider the effects of the segment, replicate reach, and season on variables, nested multivariate analysis of variance (MANOVA) was used to test whether any measured variables at the riffle scale differed between segments (UD, DD, and TR). Three segments and two replicate reaches were nested within each season (March and August) and segment (i.e., Season/Segment/Reach), whereby measurements within each reach were treated as subsamples. In the MANOVA, we also consider the interactions of the variances to interpret the interactions among the sampling segments and seasons to consider the independent effects on the factors.To perform MANOVA, we assumed that temporal variability was greater than spatial variability within each reach for variables measured over 24 h (e.g., water quality), and the opposite would hold true for variables measured only once (e.g., macroinvertebrates). Therefore, subsamples within each reach were either spatially or temporally replicated, depending on the variable type. Temporal replicates (four samples collected every 6 h) were treated as a repeated factor (time factor). A nested MANOVA was used for variables quantified once at each location (e.g., macroinvertebrates), and nested repeated-measures MANOVA (rm-MANOVA) were used for variables quantified over a 24 h period at each reach (e.g., water quality). When a significant difference was detected by MANOVA with non-significant interactions, each variable was tested separately with a nested ANOVA for variable groups once at each location or the nested rm-ANOVA for repeated-measured variables, as appropriate for the particular variable. The risk of inflating Type 1 errors for the ANOVA was reduced using Bonferroni adjustments.These MANOVA and ANOVA tests were conducted with R version 3.6.049. The residuals of each variable in each MANOVA and ANOVA model were verified using the Shapiro–Wilk normality test prior to analyses, and normality was improved using arcsine(x) or log (x + 1) transformation when appropriate.Tukey’s multiple comparison test in a one-way ANOVA model (Season/Segment/Reach) was used for comparisons between segments. Any significant changes in values for variables from UD to DD were interpreted as the effects of the dam based on the assumption that conditions in UD and DD were similar prior to dam construction; this was because replenished sediment had not been supplied in March (see before). However, UD may be unsuitable as a reference site compared with TR as the former may be at least partly affected by the dam. This may particularly be the case for benthic invertebrates, such as the interruption of the upstream flight of adult females50. Therefore, UD and TR were treated as reference sites for reservoir and tributary effects, respectively. This was because both were unaffected by the dam, and sediment replenishment as tributaries may function as sites for resource recovery for the dam-affected mainstem of the river37,51,52, despite differing watershed areas. Therefore, the similarity of variables between the TR and UD sites was statistically confirmed such that they could be treated as reference sites. As such, the recovery from March to August could reliably demonstrate the effect of sediment replenishment. For example, although the value at DD differed from that at TR and/or UD in March, it was similar to that at UD and/or TR in August.Multivariate analyses were conducted using the R “vegan” package version 2.5.6 to compare invertebrate assemblage structures between segments. Bray–Curtis coefficients based on species abundance were used to calculate a dissimilarity matrix, and dissimilarities between UD and DD, and between TR and DD in each season were tested using two-way ANOVA and Tukey post-hoc tests.Macroinvertebrate assemblage organization in relation to environmental gradients was analyzed using redundancy analysis (RDA) with the “rda” function of “vegan” package. This was because the preliminary analysis using detrended correspondence analysis (DCA) showed that the gradient lengths of DCA were More

Experimental study and protein domain analysisInsectsHosts and parasitoids were maintained as previously described25. Five Leptopilina heterotoma (Hymenoptera: Figitidae) populations were used for experiments: a population from Japan (Sapporo), two populations from the United Kingdom (1: Whittlesford; 2: Great Shelford) and two populations from Belgium (1: Wilsele; 2: Eupen). Information on collection sites, including GPS coordinates, can be found in25.Determination of host fat contentD. simulans and D. melanogaster hosts were allowed to lay eggs during 24 h in glass flasks containing ~ 50 mL standard medium25. After two days, developing larvae were sieved and ~ 200 were larvae placed in a Drosophila tube containing ~ 10 mL medium. Seven days after egg laying, newly formed pupae were frozen at – 18 °C, after which fat content was determined as described in25, where dry weight before and after neutral fat extraction was used to calculate absolute fat amount (in μg) for each host. The host pupal stage was chosen for estimating fat content, because at this point the host ceases to feed, while the parasitoid starts consuming the entire host36. All data were analysed using R Project version 3.4.360. Fat content of hosts was compared using a one-way ANOVA with host species as fixed factor.Manipulation of host fat contentTo generate leaner D. melanogaster hosts, we adapted our standard food medium25 to contain 100 times less (0.5 g) sugar per litre water. Manipulating sugar content did not alter the structure of the food medium, thus maintaining similar rearing conditions, with the exception of sugar content. Fat content of leaner and fatter D. melanogaster hosts was determined and analysed as described above.Fat synthesis quantification of wasp populationsMated female L. heterotoma were allowed to lay eggs on host fly larvae collected as described above with ad libitum access to honey as a food source until death. Honey consists of sugars and other carbohydrates that readily induce fat synthesis. After three weeks, adult offspring emergence was monitored daily and females were haphazardly placed in experimental treatments. Females were either killed at emergence (to measure teneral lipid reserves) or after feeding for 7 days on honey. Wasps were frozen at − 18 °C after completion of experiments. Fat content was determined as described above for hosts. The ability for fat synthesis was then determined by comparing mean fat levels of recently emerged compared to fed individuals, similar to procedures described in15,25,28. An increase in fat levels after feeding is indicative of active fat synthesis; equal or lower fat levels suggest fat synthesis did not take place. Each population tested on D. melanogaster or D. simulans represented an independent dataset that was analysed separately, as in Visser et al. 201825, because we are interested in the response of each population on each host species. We used T-tests when data was normally distributed and variances equal, log-transformed data for non-normal data, and a Welch’s t-test when variances were unequal. We corrected for multiple testing using Benjamini and Hochberg’s False Discovery Rate61.Fat synthesis quantification using a familial design and GC–MS analysesTo tease apart the effect of wasp genotype and host environment, we used a split-brood design where the offspring of each mother developed on lean D. simulans or fat D. melanogaster hosts in two replicated experiments (experiment 1 and 2). In both experiments, mothers were allowed to lay eggs in ~ 200 2nd to 3rd instar host larvae of one species for four days, after which ~ 200 host larvae of the other species were offered during four days. The order in which host larvae were presented was randomized across families. Following offspring emergence, daughters were allocated into two treatment groups: a control where females were fed a mixture of honey and water (1:2 w/w) or a treatment group fed a mixture of honey and deuterated water (Sigma Aldrich) (1:2 w/w; stable isotope treatment) for 7 days. Samples were prepared for GC–MS as described in 28. Incorporation of up to three deuterium atoms can be detected, but percent incorporation is highest when only 1 deuterium atom is incorporated. As incorporation of a single atom unequivocally demonstrates active fat synthesis, we only analysed percent incorporation (in relation to the parent ion) for the abundance of the m + 1 ion. Percent incorporation was determined for five fatty acids, C16:1 (palmitoleic acid), C16:0 (palmitate), C18:2 (linoleic acid), C18:1 (oleic acid), and C18:0 (stearic acid), and the internal standard C17:0 (margaric acid). Average percent incorporation for C17:0 was 19.4 (i.e. baseline incorporation of naturally occurring deuterium) and all values of the internal standard remained within 3 standard deviations of the mean (i.e. 1.6). Percent incorporation of control samples was subtracted from treatment sample values to correct for background levels of deuterium (i.e. only when more deuterium is incorporated in treatment compared to controls fatty acids are actively being synthesized). For statistical analyses, percent incorporation was first summed for C16:1, C16:0, C18:2, C18:1 and C18:0 to obtain overall incorporation levels, as saturated C16 and C18 fatty acids are direct products of the fatty acid synthesis pathway (that can subsequently be desaturated).Data (presented in Fig. 1) was analysed by means of a linear mixed effects model (GLMM, lme4 package) with host (lean D. simulans and fat D. melanogaster) and experiment (conducted twice) as fixed effect, family nested within population (Japan, United Kingdom 1 and 2, Belgium 1 and 2) as random factor, and percentage of incorporation of stable isotopes as dependent variable (log transformed; n = 138). Non-significant terms (i.e., experiment) were sequentially removed from the model to obtain the minimal adequate model as reported in Table 2. When referring to “families,” we are referring to the comparison of daughters of singly inseminated females, which (in these haplodiploid insects) share 75% of their genome.Identification of functional acc and fas genes in distinct parasitoid speciesTo obtain acc and fas nucleotide sequences for L. clavipes, G. legneri, P. maculata and A. bilineata, we used D. melanogaster mRNA ACC transcript variant A (NM_136498.3 in Genbank) and FASN1-RA (FBtr0077659 in FlyBase) and blasted both sequences against transcripts of each parasitoid (using the blast function available at http://www.parasitoids.labs.vu.nl62,63). Each nucleotide sequence was then entered in the NCBI Conserved Domain database64 to determine the presence of all functional protein domains. All sequences were then translated using the Expasy translate tool (https://web.expasy.org/translate/), where the largest open reading frame was selected for further use and confirming no stop codons were present. Protein sequences were then aligned using MAFFT v. 7 to compare functional amino acid sequences between all species (Supplementary files 1 and 2)65.Simulation studyWe consider the general situation where phenotypic plasticity is only sporadically adaptive and ask the question whether and under what circumstances plasticity can remain functional over long evolutionary time periods when the regulatory processes underlying plasticity are gradually broken down by mutations. We consider a regulatory mechanism that switches on or off a pathway (like fat synthesis) in response to environmental conditions (e.g., host fat content).Fitness considerationsWe assume that the local environment of an individual is characterized by two factors: fat content F and nutrient content N, where nutrients represent sugars and other carbohydrates that can be used to synthesize fat. Nutrients are measured in units corresponding to the amount of fat that can be synthesized from them. We assume that fitness (viability and/or fecundity) is directly proportional to the amount of fat stored by the individual. When fat synthesis is switched off, this amount is equal to F, the amount of fat in the environment. When fat synthesis is switched on, the amount of fat stored is assumed to be (N – c + (1 – k)F). This expression reflects the following assumptions: (i) fat is synthesized from the available nutrients, but this comes at a fitness cost c; (ii) fat can still be absorbed from the environment, but at a reduced rate ((1 – k)). It is adaptive to switch on fat synthesis if (N – c + (1 – k)F) is larger than F, or equivalently if (F < tfrac{1}{k}(N - c)).The right-hand side of this inequality is a straight line, which is illustrated by the blue line in Fig. 4. The three boxes in Fig. 4 illustrate three types of environmental conditions. Red box low-fat environments. Here, (F < tfrac{1}{k}(N - c)) is always satisfied, implying that fat synthesis should be switched on constitutively. Yellow box high-fat environments. Here, (F > tfrac{1}{k}(N – c)), implying that fat synthesis should be switched off constitutively.

Orange box intermediate-fat environments. Here, fat synthesis should be plastic and switched on if for the given environment (N, F) the fat content is below the blue line and switched off otherwise.

Figure 4Environmental conditions encountered by the model organisms. For a given combination of environmental nutrient content N and environmental fat content F, it is adaptive to switch on fat synthesis if (N, F) is below the blue line (corresponding to (F < tfrac{1}{k}(N - c))) and to switch it off otherwise. The three boxes illustrate three types of environment: a low-fat environment (red) where fat synthesis should be switched on constitutively; a high-fat environment (yellow) where fat synthesis should be switched off constitutively; and an intermediate-fat environment (orange) where a plastic switch is selectively favoured.Full size imageThe simulations reported here were all run for the parameters (k = tfrac{1}{2}{text{ and }}c = tfrac{1}{4}). We also investigated many other combinations of these parameters; in all cases, the results were very similar to those reported in Fig. 3.Gene regulatory networks (GRN)In our model, the switching device was implemented by an evolving gene regulatory network (as in van Gestel and Weissing66). The simulations shown in Fig. 3 of the main text are based on the simplest possible network that consists of two receptor nodes (sensing the fat and the nutrient content in the local environment, respectively) and an effector node that switches on fat synthesis if the combined weighted input of the two receptor nodes exceeds a threshold value T and switches it off otherwise. Hence, fat synthesis is switched on if (w_{F} F + w_{N} N > T) (and off otherwise). The GRN is characterized by the weighing factors (w_{F} {text{ and }}w_{N}) and the threshold T. These parameters are transmitted from parents to offspring, and they evolve subject to mutation and selection. We also considered alternative network structures (all with two receptor nodes and one effector node, but with a larger number of evolvable weighing factors67, and obtained very similar results, see below).For the simple GRN described above, the switching device is 100% adaptive when the switch is on (i.e., (w_{F} F + w_{N} N > T)) if (F < tfrac{1}{k}(N - c)) and off otherwise. A simple calculation yields that this is the case if: (w_{N} > 0{, }w_{F} = – k{kern 1pt} w_{N} {text{ and }}T = c{kern 1pt} w_{N}).Evolution of the GRNFor simplicity, we consider an asexual haploid population with discrete, non-overlapping generations and fixed population size (N = 10,000). Each individual has several gene loci, each locus encoding one parameter of the GRN. In case of the simple network described above, there are three gene loci, each with infinitely many alleles. Each individual harbours three alleles, which correspond to the GRN parameters (w_{F} {, }w_{N} {text{ and }}T), and hence determine the functioning of the genetic switch. In the simulations, each individual encounters a randomly chosen environment ((N{, }F)). Based on its (genetically encoded) GRN, the individual decides on whether to switch on or off fat synthesis. If synthesis is switched on, the individual’s fitness is given by (N – c + (1 – k)F); otherwise its fitness is given by F. Subsequently, the individuals produce offspring, where the number of offspring produced is proportional to the amount of fat stored by an individual. Each offspring inherits the genetic parameters of its parent, subject to mutation. With probability μ (per locus) a mutation occurs. In such a case the parental value (in case of a simple network: the parent’s allelic value (w_{F} {, }w_{N} {text{ or }}T)) is changed to a mutated value ((w_{F} { + }delta {, }w_{N} { + }delta {text{ or }}T + delta)), where the mutational step size δ is drawn from a normal distribution with mean zero and standard deviation σ. In the reported simulations, we chose (mu = 0.001) and (sigma = 0.1). The speed of evolution is proportional to (mu cdot sigma^{2}), implying that the rate of change in Fig. 3 (both the decay of plasticity and the rate of regaining adaptive plasticity) are positively related to μ and σ.Preadaptation of the GRNsStarting with a population with randomly initialized alleles for the GRN parameters, we first let the population evolve for 10,000 generations in the intermediate-fat environment (the orange box in Fig. 4). In all replicate simulations, a “perfectly adapted switch” (corresponding to (w_{N} > 0{, }w_{F} = – k{kern 1pt} w_{N} {text{ and }}T = c{kern 1pt} w_{N})) evolved, typically within 1,000 generations. Still, the evolved GRNs differed across replicates, as they evolved different values of (w_{N} > 0). These evolved networks were used to seed the populations in the subsequent “decay” simulations.Evolutionary decay of the GRNsFor the decay experiments reported in Fig. 3 of the main text, we initiated a large number of monomorphic replicate populations with one of the perfectly adapted GRNs from the preadaptation phase. These populations were exposed for an extended period of time (1,000,000 generations) to a high-fat environment (the yellow box in Fig. 4), where all preadapted GRNs switched off fat synthesis. However, in some scenarios, the environmental conditions changed back sporadically (with probability q) to the intermediate-fat environment (the orange box in Fig. 4), where it is adaptive to switch on fat metabolism in 50% of the environmental conditions (when (N, F) is below the blue line in Fig. 4). In Fig. 3, we report on the changing rates (q = 0.0) (no changing back; red), (q = 0.001) (changing back once every 1,000 generations; purple), and (q = 0.01) (changing back once every 100 generations; pink). When such a change occurred, the population was exposed to the intermediate-fat environment for t generations (Fig. 3 is based on t = 3).Throughout the simulation, the performance of the network was monitored every 100 generations as follows: 100 GRNs were chosen at random from the population, and each of these GRNs was exposed to 100 randomly chosen environmental conditions from the intermediate-fat environment (orange box in Fig. 4). From this, we could determine the average percentage of “correct” decisions (where the network should be switched on if and only if (F < tfrac{1}{k}(N - c)). 1.0 means that the GRN is still making 100% adaptive decisions; 0.5 means that the GRN only makes 50% adaptive decision, as would be expected by a random GRN or a GRN that switches the pathway constitutively on or off. This measure for performance in the “old” intermediate-fat environment was determined for 100 replicate simulations per scenario and plotted in Fig. 3 (mean ± standard deviation).Evolving robustness of the GRNsThe simulations in Fig. 3 are representative for all networks and parameters considered. Whenever (q = 0.0), the performance of the regulatory switch eroded in evolutionary time, but typically at a much lower rate in case of the more complex GRNs. Whenever (q = 0.01), the performance of the switch went back to levels above 90% and even above 95% for the more complex GRNs. Even for (q = 0.001), a sustained performance level above 75% was obtained in all cases.Intriguingly, in the last two scenarios the performance level first drops rapidly (from 1.0 to a much lower level, although this drop is less pronounced in the more complex GRNs) and subsequently recovers to reach high levels again. Apparently, the GRNs have evolved a higher level of robustness, a property that seems to be typical for evolving networks8. For the simple GRN studied in Fig. 3, this outcome can be explained as follows. The initial network was characterized by the genetic parameters (w_{N} > 0{, }w_{F} = – k{kern 1pt} w_{N} {text{ and }}T = c{kern 1pt} w_{N}) (see above), where (w_{N}) was typically a small positive number. In the course of evolutionary time, the relation between the three evolving parameters remained approximately the same, but (w_{N}) (and with it the other parameters) evolved to much larger values. This automatically resulted in an increasingly robust network, since mutations with a given step size distribution affect the performance of a network much less when the corresponding parameter is large in absolute value.Costs of plasticityPhenotypically plastic organisms can incur different types of costs68. In our simple model, we only consider the cost of phenotype-environment mismatching, that is, the costs of expressing the ‘wrong’ phenotype in a given environment. When placed in a high-fat environment, the preadapted GRNs in our simulations take the ‘right’ decision to switch off fat metabolism. Accordingly, they do not face any costs of mismatching. Yet, the genetic switch rapidly decays (as indicated in Fig. 3 by the rapid drop in performance when tested in an intermediate-fat environment), due to the accumulation of mutations.It is not unlikely that there are additional fitness costs of plasticity, such as the costs for the production and maintenance of the machinery underlying plasticity68. In the presence of such constitutive costs, plasticity will be selected against when organisms are living in an environment where only one phenotype is optimal (as in the high- and low-fat environments in Fig. 4). This would obviously affect the evolutionary dynamics in Fig. 3, but the size of the effect is difficult to judge, as the constitutive costs of plasticity are notoriously difficult to quantify. In case of the simple switching device considered in our model, we consider the constitutive costs of plasticity as marginal, but these costs might be substantial in other scenarios. More

A reef manta ray visits a cleaning station at Lady Elliot Island, Australia. Credit: A. O. Armstrong et al./Ecol. Evol. (CC BY 4.0)

Ecology
08 April 2021
What manta rays remember: the best spots to get spruced up

Giant fish preserve a mental map of where cleaning fish provide the highest-quality pest removal.

Share on Twitter
Share on Twitter

Share on Facebook
Share on Facebook

Share via E-Mail
Share via E-Mail

Even sea creatures need pampering. Manta rays make regular visits to ‘cleaning stations’, where small fish rid the rays of skin parasites at the coral-reef equivalent of a day spa. Now it seems that rays can identify and remember spots where they have received quality cleaning.Cleaning stations are often centred on corals inhabited by cleaner shrimp or fish. To understand how these stations influence rays’ movements, Asia Armstrong at the University of Queensland in St Lucia, Australia, and her colleagues tracked 34 reef manta rays (Mobula alfredi) off the coast of eastern Australia for about 1.5 years.The highest density of rays was found at places where cleaning fish called blue-streak cleaner wrasses (Labroides dimidiatus) were most abundant. Rays typically visited cleaning stations during the day, when cleaner wrasses are most active, and favoured stations close to foraging regions.Rays are thought to prefer stations that provide superior cleaning — where the cleaners don’t bite them, for example. The rays’ behaviour suggests that they have a mental map of spots that offer both high-quality cleaning and proximity to foraging grounds.

Ecol. Evol. (2021)

Ecology More

1.Gross, V. et al. Miniaturization of tardigrades (water bears): Morphological and genomic perspectives. Arthr. Struct. Dev. 48, 12–19 (2019).Article 

Google Scholar 
2.Møbjerg, N. et al. Survival in extreme environments – on the current knowledge of adaptations in tardigrades. Acta Physiol. 202, 409–420 (2011).Article 
CAS 

Google Scholar 
3.Giribet, G. & Edgecombe, G. D. Current understanding of Ecdysozoa and its internal phylogenetic relationships. Integr. Comp. Biol. 57, 455–466 (2017).PubMed 
Article 
PubMed Central 

Google Scholar 
4.Campbell, L. I. et al. MicroRNAs and phylogenomics resolve the relationships of Tardigrada and suggest that velvet worms are the sister group of Arthropoda. Proc. Natl Acad. Sci. USA 108, 15920–15924 (2011).CAS 
PubMed 
Article 
ADS 
PubMed Central 

Google Scholar 
5.Jørgensen, A., Møbjerg, N. & Kristensen, R. M. Phylogeny and evolution of the Echiniscidae (Echiniscoidea, Tardigrada) – an investigation of the congruence between molecules and morphology. J. Zool. Syst. Evol. Res. 49(Suppl. 1), 6–16 (2011).Article 

Google Scholar 
6.Bertolani, R. et al. Phylogeny of Eutardigrada: new molecular data and their morphological support lead to the identification of new evolutionary lineages. Mol. Phyl. Evol. 76, 110–126 (2014).Article 

Google Scholar 
7.Fujimoto, S., Jørgensen, A. & Hansen, J. G. A. molecular approach to arthrotardigrade phylogeny (Heterotardigrada, Tardigrada). Zool. Scr. 46, 496–505 (2017).Article 

Google Scholar 
8.Gąsiorek, P., Stec, D., Morek, W. & Michalczyk, Ł. Deceptive conservatism of claws: distinct phyletic lineages concealed within Isohypsibioidea (Eutardigrada) revealed by molecular and morphological evidence. Contrib. Zool. 88, 78–132 (2019).Article 

Google Scholar 
9.Hortal, J. et al. Seven shortfalls that beset large-scale knowledge of biodiversity. Ann. Rev. Ecol. Evol. Syst. 46, 523–549 (2015).Article 

Google Scholar 
10.Bartels, P. J., Apodaca, J. J., Mora, C. & Nelson, D. R. A global biodiversity estimate of a poorly known taxon: phylum Tardigrada. Zool. J. Linn. Soc. 178, 730–736 (2016).Article 

Google Scholar 
11.McInnes, S. J. Zoogeographic distribution of terrestrial/freshwater tardigrades from current literature. J. Nat. Hist. 28, 257–352 (1994).Article 

Google Scholar 
12.Morek, W., Stec, D., Gąsiorek, P., Surmacz, B. & Michalczyk, Ł. Milnesium tardigradum Doyère, 1840: The first integrative study of interpopulation variability in a tardigrade species. J. Zool. Syst. Evol. Res. 57, 1–23 (2019).Article 

Google Scholar 
13.Gąsiorek, P., Blagden, B. & Michalczyk, Ł. Towards a better understanding of echiniscid intraspecific variability: A redescription of Nebularmis reticulatus (Murray, 1905) (Heterotardigrada: Echiniscoidea). Zool. Anz. 283, 242–255 (2019).Article 

Google Scholar 
14.Gąsiorek, P. et al. Echiniscus virginicus complex: the first case of pseudocryptic allopatry and pantropical distribution in tardigrades. Biol. J. Linn. Soc. 128, 789–805 (2019).
Google Scholar 
15.Cesari, M., McInnes, S. J., Bertolani, R., Rebecchi, L. & Guidetti, R. Genetic diversity and biogeography of the south polar water bear Acutuncus antarcticus (Eutardigrada : Hypsibiidae) – evidence that it is a truly pan-Antarctic species. Invertebr. Syst. 30, 635–649 (2016).Article 

Google Scholar 
16.Guidetti, R., McInnes, S. J., Cesari, M., Rebecchi, L. & Rota-Stabelli, O. Evolutionary scenarios for the origin of an Antarctic tardigrade species based on molecular clock analyses and biogeographic data. Contrib. Zool. 86, 97–110 (2017).Article 

Google Scholar 
17.Stec, D., Krzywański, Ł, Zawierucha, K. & Michalczyk, Ł. Untangling systematics of the Paramacrobiotus areolatus species complex by an integrative redescription of the nominal species for the group, with multilocus phylogeny and species delineation in the genus Paramacrobiotus. Zool. J. Linn. Soc. 188, 694–716 (2020).Article 

Google Scholar 
18.Thulin, G. Beiträge zur Kenntnis der Tardigradenfauna Schwedens. Ark. Zool. 7, 1–60 (1911).
Google Scholar 
19.Kristensen, R. M. Generic revision of the Echiniscidae (Heterotardigrada), with a discussion of the origin of the family. In Biology of Tardigrada (ed. Bertolani, R.) 261–335 (U.Z.I. Modena, 1987).
Google Scholar 
20.Vecchi, M. et al. Integrative systematic studies on tardigrades from Antarctica identify new genera and new species within Macrobiotoidea and Echiniscoidea. Invertebr. Syst. 30, 303–322 (2016).Article 

Google Scholar 
21.Cesari, M. et al. An integrated study of the biodiversity within the Pseudechiniscus suillus–facettalis group (Heterotardigrada: Echiniscidae). Zool. J. Linn. Soc. 188, 717–732 (2020).
Google Scholar 
22.Tumanov, D. V. Analysis of non-morphometric morphological characters used in the taxonomy of the genus Pseudechiniscus (Tardigrada: Echiniscidae). Zool. J. Linn. Soc. 188, 753–775 (2020).
Google Scholar 
23.Grobys, D. et al. High diversity in the Pseudechiniscus suillus–facettalis complex (Heterotardigrada: Echiniscidae) with remarks on the morphology of the genus Pseudechiniscus. Zool. J. Linn. Soc. 188, 733–752 (2020).Article 

Google Scholar 
24.Roszkowska, M. et al. Integrative description of five Pseudechiniscus species (Heterotardigrada: Echiniscidae: the suillus-facettalis complex). Zootaxa 4763, 451–484 (2020).Article 

Google Scholar 
25.Gąsiorek, P. et al. New Asian and Nearctic Hypechiniscus species (Heterotardigrada: Echiniscidae) signalise a pseudocryptic horn of plenty. Zool. J. Linn. Soc. (in press).26.Fontoura, P. & Morais, P. Assessment of traditional and geometric morphometrics for discriminating cryptic species of the Pseudechiniscus suillus complex (Tardigrada, Echiniscidae). J. Zool. Syst. Evol. Res. 49(Suppl. 1), 26–33 (2011).Article 

Google Scholar 
27.Yu, Y., Harris, A. J. & He, X. S-DIVA (Statistical Dispersal-Vicariance Analysis): a tool for inferring biogeographic histories. Mol. Phyl. Evol. 56, 848–850 (2010).Article 

Google Scholar 
28.Matzke, N. J. Probabilistic historical biogeography: new models for founder- event speciation, imperfect detection, and fossils allow improved accuracy and model-testing. Front. Biogeogr. 5, 242–248 (2013).Article 

Google Scholar 
29.Pagel, M., Meade, A. & Barker, D. Bayesian estimation of ancestral character states on phylogenies. Syst. Biol. 53, 673–684 (2004).PubMed 
Article 
PubMed Central 

Google Scholar 
30.Phillips, S. J., Anderson, R. P. & Schapire, R. E. Maximum entropy modelling of species geographic distributions. Ecol. Model. 190, 231–259 (2006).Article 

Google Scholar 
31.Rocha, A., Doma, I., Gonzalez-Reyes, A. & Lisi, O. Two new tardigrade species (Echiniscidae, Doryphoribiidae) from Salta province (Argentina). Zootaxa 4878, 267–286 (2020).Article 

Google Scholar 
32.Ehrenberg, C. G. Diagnoses novarum formarum. Verhandl. König. Preuss. Akad. Wiss Berlin 8, 526–533 (1853).
Google Scholar 
33.Mihelčič, F. Zwei neue Tardigradenarten aus Spanien. Zool. Anz. 155, 309–311 (1955).
Google Scholar 
34.Mihelčič, F. Beitrag zur Systematik de Tardigraden . Arch. Zool. Ital. 36, 57–103 (1951).
Google Scholar 
35.Iharos, A. Zwei neue Tardigraden-Arten. Zool. Anz. 115, 218–220 (1936).
Google Scholar 
36.Murray, J. Some South African Tardigrada. J. R. Microsc. Soc. 12, 515–524 (1907).Article 

Google Scholar 
37.Yang, T. Three new species and one new record of the Tardigrada from China. Acta Hydrobiol. Sin. 26, 504–507 (2002).
Google Scholar 
38.Mihelčič, F. Beiträge zur Kenntnis der Tardigrada Jugoslawiens. Zool. Anz. 121, 95–96 (1938).
Google Scholar 
39.Bartoš, E. Eine neue Tardigradenart aus der Tschechoslowakei. Zool. Anz. 106, 175–176 (1934).
Google Scholar 
40.Richters, F. Beitrag zur Kenntnis der Moosfauna Australiens und der Inseln des Pazifischen Ozeans. Zool. Jahrb. Abt. Syst. Ökol. Geogr. Tiere 26, 196–213 (1908).
Google Scholar 
41.Vončina, K., Kristensen, R. M. & Gąsiorek, P. Pseudechiniscus in Japan: re-description of Pseudechiniscus asper Abe et al., 1998 and description of Pseudechiniscus shintai sp. nov. Zoosyst. Evol. 96, 527–536 (2020).Article 

Google Scholar 
42.Wang, L. Tardigrades from the Yunnan-Guizhou Plateau (China) with description of two new species in the genera Mixibius (Eutardigrada: Hypsibiidae) and Pseudechiniscus (Heterotardigrada: Echiniscidae). J. Nat. Hist. 43, 2553–2570 (2009).Article 

Google Scholar 
43.Hulings, N. C. & Gray, J. S. A manual for the study of meiofauna. Smithson. Contrib. Zool. 78, 1–84 (1971).
Google Scholar 
44.Gąsiorek, P. & Michalczyk, Ł. Revised Cornechiniscus (Heterotardigrada) and new phylogenetic analyses negate echiniscid subfamilies and tribes. R. Soc. Open Sci. 7, 200581 (2020).PubMed 
PubMed Central 
Article 
ADS 
CAS 

Google Scholar 
45.Dastych, H. Notes on the revision of the genus Mopsechiniscus (Tardigrada). Zool. Anz. 240, 299–308 (2001).Article 

Google Scholar 
46.Rebecchi, L., Altiero, T., Eibye-Jacobsen, J., Bertolani, R. & Kristensen, R. M. A new discovery of Novechiniscus armadilloides (Schuster, 1975) (Tardigrada, Echiniscidae) from Utah, USA with considerations on non-marine Heterotardigrada phylogeny and biogeography. Org. Divers. Evol. 8, 58–65 (2008).Article 

Google Scholar 
47.Binda, M. G. & Kristensen, R. M. Notes on the genus Oreella (Oreellidae) and the systematic position of Carphania fluviatilis Binda, 1978 (Carphaniidae fam. nov., Heterotardigrada). Animalia 13, 9–20 (1986).
Google Scholar 
48.Binda, M. G. Risistemazione di alcuni Tardigradi con l’instituzione di un nuovo genere di Oreellidae e della nuova famiglia Archechiniscidae. Animalia 5, 307–314 (1978).
Google Scholar 
49.Kristensen, R. M. & Hallas, T. E. The tidal genus Echiniscoides and its variability, with erection of Echiniscoididae fam. n. (Tardigrada). Zool. Scr. 9, 113–127 (1980).Article 

Google Scholar 
50.Møbjerg, N., Kristensen, R. M. & Jørgensen, A. Data from new taxa infer Isoechiniscoides gen. nov. and increase the phylogenetic and evolutionary understanding of echiniscoidid tardigrades (Echiniscoidea: Tardigrada). Zool. J. Linn. Soc. 178, 804–818 (2016).Article 

Google Scholar 
51.Møbjerg, N., Jørgensen, A. & Kristensen, R. M. Ongoing revision of Echiniscoididae (Heterotardigrada: Echiniscoidea), with the description of a new interstitial species and genus with unique anal structures. Zool. J. Linn. Soc. 188, 663–680 (2020).Article 

Google Scholar 
52.Gąsiorek, P., Suzuki, A. C., Kristensen, R. M., Lachowska-Cierlik, D. & Michalczyk, Ł. Untangling the Echiniscus Gordian knot: Stellariscus gen. nov. (Heterotardigrada: Echiniscidae) from Far East Asia. Invertebr. Syst. 32, 1234–1247 (2018).Article 

Google Scholar 
53.Dastych, H. Echiniscus rackae sp. n., a new species of Tardigrada from the Himalayas. Entomol. Mitt. Zool. Mus. Hamburg 8, 246–250 (1986).
Google Scholar 
54.McInnes, S. J. Tardigrades from Signy Island, South Orkney Islands, with particular reference to freshwater species. J. Nat. Hist. 29, 1419–1445 (1995).Article 

Google Scholar 
55.Dastych, H. Two new species of Tardigrada from the Canadian Subarctic with some notes on sexual dimorphism in the family Echiniscidae. Entomol. Mitt. Zool. Mus. Hamburg 8, 319–334 (1987).
Google Scholar 
56.Pilato, G., Binda, M. G. & Lisi, O. Remarks on some Echiniscidae (Heterotardigrada) from New Zealand with the description of two new species. Zootaxa 1027, 27–45 (2005).Article 

Google Scholar 
57.Pilato, G., Binda, M. G., Napolitano, A. & Moncada, E. Notes on South American tardigrades with the description of two new species: Pseudechiniscus spinerectus and Macrobiotus danielae. Trop. Zool. 14, 223–231 (2001).Article 

Google Scholar 
58.Fontoura, P., Pilato, G. & Lisi, O. First record of Tardigrada from São Tomé (Gulf of Guinea, Western Equatorial Africa) and description of Pseudechiniscus santomensis sp. nov. (Heterotardigrada: Echiniscidae). Zootaxa 2564, 31–42 (2010).Article 

Google Scholar 
59.Bartoš, E. Die Tardigraden der Chinesischen und Javanischen Moosproben. Acta Soc. Zool. Bohem. 27, 108–114 (1963).
Google Scholar 
60.Beijerinck, M.W. De infusies en de ontdekking der backteriën. Jaarboek van de Koninklijke Akademie v. Wetenschappen. Amsterdam: Müller (1913).61.Baas-Becking, L. G. M. Geobiologie of inleiding tot de milieukunde (W.P. Van Stockum & Zoon, 1934).
Google Scholar 
62.Wallace, A. R. The geographical distribution of animals: with a study of the relations of living and extinct faunas as elucidating the past changes of the Earth’s surface (Macmillan and Company, 1876).
Google Scholar 
63.Niedbała, W. The ptyctimous mites fauna of the Oriental and Australian regions and their centres of origin (Acari: Oribatida). Genus Suppl. 10, 1–493 (2000).
Google Scholar 
64.Niedbała, W. Ptyctimous mites (Acari: Oribatida) of South Africa. Ann. Zool. 56(Suppl. 1), 1–97 (2006).
Google Scholar 
65.Janion-Scheepers, C., Deharveng, L., Bedos, A. & Chown, S. Updated list of Collembola species currently recorded from South Africa. ZooKeys 503, 55–88 (2015).Article 

Google Scholar 
66.Kisielewski, J. Inland-water Gastrotricha from Brazil. Ann. Zool. 43(Suppl. 2), 1–168 (1991).
Google Scholar 
67.Tuxen, S. L. Ecology and zoogeography of the Brazilian Protura (Insecta). Stud. Neotrop. Fauna Environ. 12, 225–247 (1977).Article 

Google Scholar 
68.Greenslade, P. Why are there so many exotic springtails in Australia? A review. Soil Org. 90, 141–156 (2018).
Google Scholar 
69.Smit, H. Australian water mites of the subfamily Notoaturinae Besch (Acari: Hydrachnidia: Aturidae), with the description of 24 new species. Int. J. Acarol. 36, 101–146 (2010).Article 

Google Scholar 
70.Moir, M. L., Brennan, K. E. C. & Harvey, M. S. Diversity, endemism and species turnover of millipedes within the south-western Australian global biodiversity hotspot. J. Biogeogr. 36, 1958–1971 (2009).Article 

Google Scholar 
71.Harvey, M. S., Abrams, K. M., Beavis, A. S., Hillyer, M. J. & Huey, J. A. Pseudoscorpions of the family Feaellidae (Pseudoscorpiones: Feaelloidea) from the Pilbara region of Western Australia show extreme short-range endemism. Invertebr. Syst. 30, 491–508 (2016).Article 

Google Scholar 
72.Claxton, S.K. The taxonomy and distribution of Australian terrestrial tardigrades. PhD thesis, Macquarie University: Sydney (2004).73.Simpson, G. G. Tempo and mode in evolution (Columbia University Press, 1944).
Google Scholar 
74.Dastych, H. The Tardigrada of Poland. Monogr. Faun. Pol. 16, 1–255 (1988).
Google Scholar 
75.Peters, M. K. et al. Predictors of elevational biodiversity gradients change from single taxa to the multi-taxa community level. Nat. Commun. 7, 13736 (2016).CAS 
PubMed 
PubMed Central 
Article 
ADS 

Google Scholar 
76.Petersen, B. The tardigrade fauna of Greenland. Medd. Grønl. 150, 1–94 (1951).
Google Scholar 
77.de Bruyn, M. et al. Borneo and Indochina are major evolutionary hotspots for Southeast Asian biodiversity. Syst. Biol. 63, 879–901 (2014).PubMed 
Article 
PubMed Central 

Google Scholar 
78.Ugland, K. I., Gray, J. S. & Ellingsen, K. E. The species–accumulation curve and estimation of species richness. J. Anim. Ecol. 72, 888–897 (2003).Article 

Google Scholar 
79.Casquet, J. T., Thebaud, C. & Gillespie, R. G. Chelex without boiling, a rapid and easy technique to obtain stable amplifiable DNA from small amounts of ethanol-stored spiders. Mol. Ecol. Resour. 12, 136–141 (2012).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
80.Stec, D., Kristensen, R.M. & Michalczyk, Ł. An integrative description of Minibiotus ioculator sp. nov. from the Republic of South Africa with notes on Minibiotus pentannulatus Londoño et al., 2017 (Tardigrada: Macrobiotidae). Zool. Anz. 286, 117–134 (2020).81.Hall, T. A. BioEdit: a user-friendly biological sequence alignment editor and analysis program for Windows 95/98/NT. Nucleic Acids Symp. Ser. 41, 95–98 (1999).CAS 

Google Scholar 
82.Nguyen, L.-T., Schmidt, H. A., von Haeseler, A. & Minh, B. Q. IQ-TREE: A fast and effective stochastic algorithm for estimating maximum likelihood phylogenies. Mol. Biol. Evol. 32, 268–274 (2015).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
83.Hoang, D. T., Chernomor, O., von Haeseler, A., Minh, B. Q. & Vinh, L. S. UFBoot2: Improving the ultrafast bootstrap approximation. Mol. Biol. Evol. 35, 518–522 (2018).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
84.Kalyaanamoorthy, S., Minh, B. Q., Wong, T. K. F., von Haeseler, A. & Jermiin, L. S. ModelFinder: Fast model selection for accurate phylogenetic estimates. Nat. Meth. 14, 587–589 (2017).CAS 
Article 

Google Scholar 
85.Zhang, J., Kapli, P., Pavlidis, P. & Stamatakis, A. A general species delimitation method with applications to phylogenetic placements. Bioinformatics 29, 2869–2876 (2013).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
86.Katoh, K., Misawa, K., Kuma, K. & Miyata, T. MAFFT: a novel method for rapid multiple sequence alignment based on fast Fourier transform. Nucleic Acids Res. 30, 3059–3066 (2002).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
87.Katoh, K. & Toh, H. Recent developments in the MAFFT multiple sequence alignment program. Brief. Bioinform. 9, 286–298 (2008).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
88.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 
89.Vaidya, G., Lohman, D. J. & Meier, R. SequenceMatrix: concatenation software for the fast assembly of multi-gene datasets with character set and codon information. Cladistics 27, 171–180 (2011).Article 

Google Scholar 
90.Chernomor, O., von Haeseler, A. & Minh, B. Q. Terrace aware data structure for phylogenomic inference from supermatrices. Syst. Biol. 65, 997–1008 (2016).PubMed 
PubMed Central 
Article 

Google Scholar 
91.Ronquist, F. & Huelsenbeck, J. P. MrBayes 3: Bayesian phylogenetic inference under mixed models. Bioinformatics 19, 1572–1574 (2003).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
92.Suchard, M.A. et al. Bayesian phylogenetic and phylodynamic data integration using BEAST 1.10 Virus Evol. 4, vey016 (2018).93.Rambaut, A., Suchard, M.A., Xie, D. & Drummond, A.J. Tracer v1.6 (2014). Available from http://beast.bio.ed.ac.uk/Tracer.94.Drummond, A. J. & Suchard, M. A. Bayesian random local clocks, or one rate to rule them all. BMC Biol. 8, 114 (2010).PubMed 
PubMed Central 
Article 

Google Scholar 
95.Ferreira, M. A. R. & Suchard, M. A. (2008) Bayesian analysis of elapsed times in continuous-time Markov chains. Can. J. Stat. 36, 355–368 (2008).MATH 
Article 

Google Scholar 
96.Münkemüller, T. et al. How to measure and test phylogenetic signal. Meth. Ecol. Evol. 3, 743–756 (2012).Article 

Google Scholar 
97.Yu, Y., Harris, A. J., Blair, C. & He, X. J. RASP (Reconstruct Ancestral State in Phylogenies): a tool for historical biogeography. Mol. Phyl. Evol. 87, 46–49 (2015).Article 

Google Scholar 
98.Yu, Y., Blair, C. & He, X. J. RASP 4: ancestral state reconstruction tool for multiple genes and characters. Mol. Biol. Evol. 37, 604–606 (2020).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
99.Jombart, T., Balloux, F. & Dray, S. adephylo: new tools for investigating the phylogenetic signal in biological traits. Bioinformatics 26, 1907–1909 (2010).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
100.Pennell, M. W. et al. geiger v2.0: an expanded suite of methods for fitting macroevolutionary models to phylogenetic trees. Bioinformatics 30, 2216–2218 (2014).CAS 
PubMed 
Article 
PubMed Central 

Google Scholar 
101.Matzke, N. J. Model selection in historical biogeography reveals that founder-event speciation is a crucial process in island clades. Syst. Biol. 63, 951–970 (2014).PubMed 
Article 
PubMed Central 

Google Scholar 
102.Ree, R. H. & Sanmartín, I. Conceptual and statistical problems with the DEC+J model of founder-event speciation and its comparison with DEC via model selection. J. Biogeogr. 45, 741–749 (2018).Article 

Google Scholar 
103.Phillips, S. J., Anderson, R. P., Dudík, M., Schapire, R. E. & Blair, M. E. Opening the black box: an open-source release of Maxent. Ecography 40, 887–893 (2017).Article 

Google Scholar 
104.Cobos, M. E., Peterson, A. T., Barve, N. & Osorio-Olvera, L. kuenm: an R package for detailed development of ecological niche models using Maxent. PeerJ 7, e6281 (2019).PubMed 
PubMed Central 
Article 

Google Scholar 
105.R Core Team. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org/ (2020).106.Fick, S. E. & Hijmans, R. J. WorldClim 2: new 1km spatial resolution climate surfaces for global land areas. Int. J. Climatol. 37, 4302–4315 (2017).Article 

Google Scholar 
107.Escobar, L. E., Lira-Noriega, A., Medina-Vogel, G. & Peterson, A. T. Potential for spread of the white-nose fungus (Pseudogymnoascus destructans) in the Americas: use of Maxent and NicheA to assure strict model transference. Geospat. Health 9, 221–229 (2014).PubMed 
Article 
PubMed Central 

Google Scholar 
108.Peterson, A. T., Papeş, M. & Soberón, J. Rethinking receiver operating characteristic analysis applications in ecological niche modeling. Ecol. Model. 213, 63–72 (2008).Article 

Google Scholar 
109.Anderson, R. P., Lew, D. & Peterson, A. T. Evaluating predictive models of species’ distributions: criteria for selecting optimal models. Ecol. Model. 162, 211–232 (2003).Article 

Google Scholar 
110.QGIS Development Team. QGIS Geographic Information System. Open Source Geospatial Foundation Project (2020). More