Ecology

Study lakesThe study was conducted in two water bodies created after aquatic restorations of mining pits, Lakes Milada (2.5 km2, high structural complexity; 50° 39′ N, 13° 58′ E) and Most (3.1 km2, low structural complexity lake; 50° 32′ N, 13° 38′ E), in the Czech Republic (Fig. 1). Aquatic restoration took place from 2001 to 2010 in Milada and from 2008 to 2014 in Most. Both lakes are medium-sized (surface area = 252 and 311 hectares, respectively), relatively deep (Zmean = 16 and 22 m, Zmax = 25 and 75 m), oligotrophic (mean summer total phosphorus  More

Because biological materials are often viscoelastic composites, with properties dependent on orientation as well as spatial and temporal scales, all tests were designed to mimic natural conditions of tool use. For example, because of possible anisotropies, indentations were made on contact surfaces instead of cross sections of the “tools”; abrasion resistance was measured in the contact direction. Fracture tests were designed to mimic the tension on the side of a tooth or sting subjected to lateral forces; the velocities in our impact tests were between 0.1 and 1 m/s, typical for tool interactions, though 1/10 as fast as for some organisms48,49.OrganismsAll organisms, except the salmon, were housed alive in our laboratory until just before testing. Samples were measured within 24 h of removal from the organism, and were maintained in high-humidity environments during the preparation period. This is particularly important because mechanical properties can depend strongly on water content and this dependence can differ between materials. For example, non-HEBs can harden more than HEBs with drying30,50,51.

1.

Leafcutter ants, Atta cephalotes, were obtained from colonies we collected in Flores, Guatemala and Arena Forest Reserve, Trinidad. They were maintained on Himalayan blackberry leaves, Rubus armeniacus, Portugese laurel, Prunus lusitanica, and Japanese spurge, Pachysandra terminalis.

2.

Nereid worms, Neanthes brandti (synonymous with Alitta brandti), were collected from trenches dug in bars near the mouth of the Coos River, Charleston, Oregon. They were kept in sea-water in a laboratory refrigerator at about 5 °C.

3.

Scorpions, Hadrurus arizonensis, from Arizona, were obtained from commercial suppliers (e.g. Bugs of America, http://bugsofamerica.com) and fed House crickets, Acheta domesticus, and mealworms, Tenebrio molitor larvae.

4.

Spiders, Araneus diadematus, were collected seasonally around the University of Oregon campus, and kept alive in a laboratory refrigerator at about 5 °C. Tarantulas, Aphonopelma hentzi were obtained from Carolina (https://www.carolina.com).

5.

Chitons, Katharina tunicata, and Cryptochiton stelleri were collected in rocky intertidal regions along the coast near Charleston, Oregon and kept alive in a laboratory refrigerator.

6.

Salmon heads, Oncorhynchus sp., were obtained fresh from seafood stores.

7.

Leaf cutter bees, Megachile rotundata, were obtained from Mason Bees for Sale (www.masonbeesforsale.com).

Hardness, modulus of elasticity and damping measurementsHardness, modulus of elasticity and dynamic mechanical property measurements were made by pressing a sharp diamond probe into specimens and measuring the resulting indentation as it changed in time. A higher modulus of elasticity indicates that a structure is stiffer and suffers less elastic (quickly recovered) deformation. A higher hardness value indicates that the material will undergo less plastic (non-recovered) deformation and thus will have a smaller pit left behind after the indentation. A material with a higher loss tangent will absorb more energy of vibration (higher damping), and is characterized by lagging surface deformation and recovery (viscoelasticity) as the indention force changes. Damping can reduce damage because the energy absorbed and converted to heat is not available for breaking bonds in the material.We used an Atomic Force Microsocpe (AFM; NanoScope IIIa, Digital Instruments, Santa Barbara, CA) with an add-on force/displacement transducer (TriboScope, Hysitron Inc., Minneapolis, Minnesota). The Hysitron transducer held a polished diamond probe in place with capacitors that were used to sense the position of the probe and to impart vertical forces for indenting and imaging the specimen. Measurement regions were selected for minimal slope and surface topography as evidenced by the depth variation in AFM scans and the symmetry of the residual indents.In order to make the most biologically relevant measurements, indentations were made on regions of the external surface of structures that directly contact the environment. However, we also made measurements on cross-sections of the arthropod “tools” as part of preliminary SEM and indentation investigations to ensure that the thickness of epicuticle or other surface material would not distort the results for HEBs. Figure 2A shows an indentation on an un-polished surface—the original surface topography is visible as linear scratches that are small compared to the pyramidal indentation. When we could not avoid indentation-scale topography, we hand-polished the surface with 2000–12,000 grit sandpaper (Micro-mesh sheets, http://micro-surface.com), to smooth the surface on the scale of the test indents.Figure 2Images of testing samples for each of the measured properties. (A) AFM image of a residual indentation on the natural surface of the Mn-HEB region of the sting of a scorpion, Hadrurus arizonensis. Indentations were used in measuring hardness, modulus of elasticity and damping properties. Natural “scratches” are visible on the original surface around the triangular indentation made by the cube-cornered indenter. (B) Mandible of an ant, Atta cephalotes, before and after an abrasion testing session. The tip of the zinc-rich distal tooth in the “before” image has been flattened during a preliminary abrasion session. (C) Images from before and after energy of fracture testing of a 12 μm thick Zn-HEB test piece made from the nearly-flat side-surface of the fang of a spider, Araneus diadematus. The original fang shape is evident with the tip of the fang towards the top and the proximal side to the left. (D) Images from before and after impact resistance testing of a 12 μm thick Zn-HEB test piece made from the fang of a spider, Araneus diadematus. The test piece covers a 50 μm diameter backing-pit milled (using FIB-SEM) into a (reflective) silicon chip. The piece has been shattered from the impact in the image on the right.Full size imageSpecimens were mounted on atomic force microscopy (AFM) specimen disks (TedPella Inc., Redding, CA) in a mound of epoxy composite. The composite was prepared by mixing approximately 0.45 g of 400-grit aluminum oxide powder (Buehler, Evanston, Ill., Ted Pella Inc., Redding CA, or Kramer Industries, Piscataway, NJ—the later preferred because it was less reflective) with 0.075 mL each of resin and hardener (Quick Set Epoxy; Loctite, Rocky Hill, CT, and 5 m Quik-Cure Epoxy, Bob Smith Industries, Atascadero CA). The composite was stiff enough that the unpolished specimen had to be pressed in and could be oriented before curing so that the desired indentation region would retro-reflect a light beam sent through the eyepiece of a dissecting microscope back into the microscope, ensuring that the desired region would be flat for AFM scanning and indentation. The mounted specimens were placed in an oven at 39 °C for at least 1 h to cure the epoxy composite. Samples that were not tested immediately were kept on moist paper in a container in a refrigerator and were tested within 24 h to avoid dehydration and other changes. Additionally, the epoxy composite served as a barrier to reduce loss of water through the cut surfaces.The epoxy composite mounting technique was tested to ensure that small (about the same size as the biological specimens) “floating” glass cover-slip pieces would yield the same hardness and modulus of elasticity values as large, flat-mounted cover-slip pieces, and re-checked when relevant products were changed.In order to test whether there were any rapid changes in hardness or modulus of elasticity of HEBs, we tested a H. arizonensis sting, that was mounted for AFM measurements while still attached to an anaesthetized scorpion. We did not find a significant trend with time for 14 measurements made between 10 or 20 min (Zn-HEB and Mn-HEB respectively) after separating the live scorpion from the sting, and 5 h (largest R2: 0.006). These fast results were also comparable to results made using the standard technique, indicating that the technique described above was sufficient to prevent significant changes from dehydration. We also made preliminary measurements on scorpion joint cuticle, armour teeth, and other non-HEB regions of the cuticle, and found none that were harder than the region at the base of the sting, used here to represent non-HEB cuticle.We used a pyramid-shaped diamond probe with cubic corner facets (90° between the three faces)14,31. The steeper angle of the cubic tip, relative to a more commonly used Berkovich tip, made it easier to avoid surface features such as hairs. The diamond probe was positioned on the specimen using a 30× extra short focus monocular (M1030, Specwell Corporation, Tokyo, Japan). The indentation sequence began with the force being ramped linearly from 0 to 2 milliNewtons (mN) in 0.1 s, maintained at 2 mN for 10 s. The force was then ramped down to 1.5 mN over 0.1 s and then the force was varied sinusoidally at 10 Hz (for 25 cycles) with a peak-to-peak amplitude of 1.0 mN, in order to measure dynamic properties. The force was then ramped to 0 in 0.1 s.Probe-extension (Oliver–Pharr) and image-based measurementsTwo methods of obtaining the modulus of elasticity and hardness were employed14. In the first method, values were obtained only from force–displacement curves using the Oliver–Pharr technique52,53. The modulus of elasticity was obtained from the slope of the force–displacement curve at the beginning of withdrawal of the indenter. Oliver–Pharr hardness values were calculated from the intercept of this sloped line with the line of zero force.In addition to the hardness value obtained from the force–displacement curves, we also obtained a hardness value based on measurements of the size of the residual indentation. These image-based hardness values were calculated as H = F/A, where “F” is the maximum force applied to the probe, and “A” is the projected area of the residual indentation, obtained from the perimeter of the indentation measured on an AFM image (e.g. Fig. 2A) made by scanning the indenting probe itself minutes after indenting the specimen.The Oliver–Pharr method is inaccurate if the indentation force causes the surface of the specimen to move, such as for improperly backed specimens, because it assumes probe extension is a measure of indentation depth. In contrast, the image-based method is nearly insensitive to global displacements of the specimen because it is based only on the applied force and measurements of the residual indentation. We found it useful to obtain both values to check each other: on several occasions differences between the two measured values indicated support problems. This is important for biological specimens with multiple layers, voids, etc. For example, if there is a lumen under the shaft but not the tip, the tip may appear harder because it displaces less. In addition, calibration problems, such as from fractured silica, were quickly identified by differences in the image and Oliver–Pharr values. Finally, the image method is not sensitive to other artifacts, such as “pile-up”, that are associated with estimating contact region from probe extension53,54,55,56.There is a potential difference between the Oliver–Pharr and image hardness values associated with the different time scales. The Oliver–Pharr hardness is measured during the probe withdrawal while the image of the residual indent is obtained a couple of minutes after indentation. If the indent partially recovers in the interim, the image technique would be based on a smaller residual indentation resulting in a higher hardness value. We prefer the image method not only because it is robust to imperfectly supported specimens, but also because we would like our hardness measurement to reflect the long-term indentation damage done to the tips and blades. To test that the indent size had stabilized by the time we measured it, we re-measured an indentation in the zinc-region of a scorpion sting after more than 6 months and found that the indentation diameter had decreased little, by about 15%.We also calculated an image-based modulus of elasticity value that has been suggested for materials that produce “pile-up” artifacts53. The area measured from the image of the residual indentation (used for the image-based hardness), was substituted for the contact area calculated from probe extension in Eq. (6) of Oliver & Farr, 2004 53.Notwithstanding the differences in image-based and probe-extension based (Oliver–Pharr) measurements, there was little practical difference in results, as shown in Fig. 3A,B. The metals and plastics that we measured for Fig. 3 tended to have slightly higher Oliver–Pharr hardness values than image-based values, possibly because of “pile-up” artifacts.Figure 3Comparison of probe-extension (Oliver–Pharr) and image-based techniques for hardness (A) and reduced modulus of elasticity (B) for our indentation data. Pile-up may account for higher (above the line) Oliver–Pharr values of some of the metals and plastics.Full size imageIn the “Results” section, we plot the values obtained using the images, but the Oliver–Pharr values are included along with image-based values in the results table, Table 1.Loss tangentDynamic mechanical properties were measured from the sinusoidal segments of the indentation sequence by comparing the amplitude and phase of the displacement to the applied sinusoidal force57. Phase lags associated with the transducer and electronics were determined assuming zero true-lag from a fused silica standard obtained using the same indentation sequence and force. The loss tangents obtained in this way were for the high-stress regimes associated with indentation, as compared to low-stress tests involving bending without plastic deformation.Calibration for Oliver–Pharr measurementsIn order to obtain the contact area from the indent depth, the shape of the indenting tip must be known. We characterized the indenting tip shape directly using scanning electron microscope images, so that we could make deep, micron-scale indents that would be evident on un-polished biological surfaces. We could not make indents in silica as deep as the indents desired for our biological specimens without fracturing the silica (fracture for the cubic cornered tip began at about 3 mN) so we could not use the usual technique of estimating the tip shape at depth by calibrating with fused silica. The tip shape was characterized by three measurements: first, the angle of the three-sided pyramidal tip (α), second, a measure of the bluntness of the tip (B), the distance between the apex of the tip if it were an ideal pyramid and the actual blunt tip, and third, a measure of the distance from the blunt tip beyond which the shape of the tip was not distinguishable from an ideal pyramid (I). The value “I” was used as a limit: only indents with greater depth than “I” were used to calculate mechanical properties. For these deeper indents, the following description of the projected area (A) of the contact region between the tip and the specimen was used:$$A= frac{(0.433)(4){(D+B)}^{2}}{frac{1}{{mathrm{tan}}^{2}left(frac{alpha }{2}right)}-0.3333},$$where D is the depth of the indent, determined by the extension of the indenting probe, 0.433 is the ratio of the area of an equilateral triangle to the square of the length of a side, and 0.3333 is tan2 (30°). As an example, the tip used for the majority of measurements was characterized by α = 89.9°, B = 115 nm, I = 100 nm. Thus, for indents with a depth greater than 100 nm, for our tip, A = 2.58 (D + 115 nm)2.Measurement of residual indentation areaThe area of the residual indent was measured using an AFM image, obtained minutes after the indentation, using the indenting tip as the imaging tip. To minimize inaccuracies in indent perimeter determination, caused by finite size of the imaging probe or other systematic errors, we calibrated our area measurements so that we obtained a median value of 70 GPa for measurements of the modulus of elasticity on fused silica. Because there is some subjectivity in measuring the size of the indentation, operator-specific calibrations, based on each operator’s measurements of fused silica, were used for most of the measurements.Test piece preparation for impact resistance and energy of fracture measurementsWe measured resistance to impact and fracture using custom miniature versions of testing devices that fracture or damage standardized “test pieces” of materials. We prepared test pieces as follows: the fresh (usually immediately after removal from the organism) specimens were adhered to one end of a glass slide using a marine epoxy (Loctite, Rocky Hill, CT) which required a curing time of 2 h at 39 °C, or with cyanoacrylate adhesive (Krazy Glue, all purpose, Elmer’s Products)14, which required no extra curing time. A flat region ( > 100 μm diameter, but not wide enough to reduce the thickness of the HEB region in the center to less than 12 μm) was polished on a specimen by grinding the slide with the specimen against a sequence of flat 2000, 6000 and 12,000 grit sandpaper (Micro-mesh sheets, http://micro-surface.com). The specimen was removed from the adhesive using a scalpel, inverted, and the polished-flat region was adhered to the glass with a thin film of water and surrounded by a small bead of marine epoxy. The water film kept the epoxy from getting pulled under the specimen by capillary action. The epoxy was cured and the specimen polished to a thickness of 12 ± 2 μm as determined with a digital micrometer, using the same sandpaper sequence. The resulting test piece was then freed by scraping the epoxy from around the edges using a scalpel blade. The area of the pieces varied according to the size of flat regions and was, typically, hundreds of microns on a side.Maintaining hydration was especially important for these test pieces because they were only 12 μm thick and so they could dry quickly. To reduce artifacts from drying or other changes in the tested material, all preparation and testing took place in a ~ 15 m3 enclosure maintained at greater than 90% relative humidity.Although we used FIB-SEM (Focused Ion Beam-Scanning Electron Microscope) to shape specimens of the materials for molecular fragment analysis, we did not use this technique for preparing our micron-scale test pieces for several reasons: potential material property changes caused by beam damage, subjection to vacuum, and because some of the test pieces needed to be large enough that they would be difficult to make with FIB-SEM.Impact resistance measurementsA custom testing device was built to compare the energy required for a swinging pendulum to shatter test pieces of the different materials (Fig. 2D). A 12 ± 2 μm thick test piece was adhered by the moisture in the high-humidity enclosure and held in place with an adhesive (spots of cyanoacrylate, Krazy Glue, all purpose, Elmer’s Products, www.elmers.com, or 5-min epoxy gel) over a 50 μm-diameter circular pit milled in a silicon wafer using a FIB-SEM apparatus. The pendulum, made of carbon fiber and aluminum (length: 0.2 m, moment of inertia: 4.25 x 10−6 kg m2) with a diamond impactor tip polished to a diameter of 20 μm, was held by miniature bearings and electronically released from increasing heights until the test piece fractured. The energy required to fracture the specimen was calculated from the release height from which the pendulum fractured the test piece. This energy was normalized by the measured thickness of the specimens to give joules required per meter of thickness. Nevertheless, we consider this test to be a relative test that is not expected to be generalizable to all impacts, as the energy to fracture is likely to depend not only on thickness but also on variables such as the diameter of the impactor tip and the diameter of the backing hole.This impact test differs from Charpy and Izod tests in that the energy required to fracture the specimen was measured by releasing the pendulum from increasing heights until the specimen fractured, rather than by releasing it from a height sufficient to fracture all specimens, and measuring the residual energy of the pendulum after impact. The advantage of our threshold technique is that the threshold of fracture is likely the biologically important quantity, and direct determination of the fracture threshold avoids the possibility that the energy deposited in a single highly-energetic impact might be partially expended in plastic deformation, leading to an overestimate of the threshold energy. A drawback of our technique is that impacts that do not break the specimen may produce damage that weakens the specimen for subsequent impacts. Nevertheless, all specimens were subjected to the same series of increasingly energetic impacts, until fracture, and were thus comparable.Energy of fracture measurementsWe measured the energy or work required to slowly break a test piece in two (Fig. 2c) using a custom fracture toughness measuring device14. The device drove apart two microscope cover slips bridged by the test piece until it split in two, while recording the required force and the displacement (work is the product of force and incremental distance). This work, divided by the area of the new post-fracture surfaces, is the energy of fracture, reported in Joules per meter squared. It is a measure of the energy required, per unit area, to break the bonds that originally held the two pieces together (as long as the kinetic energy is relatively small—the pieces do not fly away), and is one indication of the resistance of a material to fracture.The length of the fracture was measured using a microscope and multiplied by the thickness of the test piece to obtain the fracture area. This area was used to normalize force–displacement curves. The work of fracture per unit area of the fracture was obtained by numerically integrating these normalized force–displacement curves. The load cell was designed to be stiff in order to minimize storage of energy within the apparatus as the specimen underwent tension58. Fracture planes were perpendicular to the original surface and approximately perpendicular to the long axis of the “tool” (Fig. 1C).The test protocol was altered from that used previously14 because the test pieces used here were smaller. Test pieces were not notched in order to avoid fractures from the notching process, and the specimens were adhered to the test apparatus in place (Krazy Glue, all purpose, Elmer’s Products, www.elmers.com) to avoid premature fracture. To improve the bonding of the cyanoacrylate adhesive to the glass cover slips in the high humidity atmosphere, we treating the cover slips with a 10-s dip in a 2% (by volume) 3-aminopropyltriethoxysilane (Sigma Chemical Co., www.sigmaaldrich.com), 98% acetone solution, followed by rinses in deionized water and air drying. The test protocol for the ant mandibular teeth varied from the others in that whole teeth were fractured instead of 12μ polished test pieces.A consistency test with AFM data was developed to identify cases of imperfect bonding to the cover slips, when part of the measured energy was expended in partially pulling the test piece out of the cyanoacrylate adhesive. For samples subject to this problem (usually specimens with small adhesive contact areas, such as the fang specimen in Fig. 1C), we required that the force–displacement curve be consistent with the stiffness of the test piece, expected from a model based on the shape of the individual test piece and the slope of the force–displacement curve for the insertion portion of the nano-indentation sequence for that material. When a piece partially pulled out and failed this test, the apparent stiffness was much lower than the expected stiffness (from nano-indentation) and slight stretch marks were often visible in the adhesive on close inspection.This stiffness consistency test was also found to be useful in identifying cases where part of the fracture was pre-existing but had not been visible in the test piece. The pre-existing fracture would tend to reduce the effective width of the specimen and thus could be identified by a lower than expected stiffness under tension.Abrasion resistance measurementsWe measured the energy required to abrade away a volume of material from our specimens by holding them against a rotating abrasive disk. The energy used in eroding the material is given by the force of friction multiplied by the distance traveled over the abrasive paper (work is the product of the force and incremental displacement), with units of Joules per meter cubed of volume worn away.The “pin on disk” type testing device, developed for testing pieces of crab cuticle14, was used with modified procedures for the smaller specimens here. Instead of cylindrical core samples, whole, approximately conical tips of teeth, fangs or stings were used (Fig. 2B). The samples were affixed with cyanoacrylate gel adhesive (Maxi-Cure, Bob Smith Industries, Atascadero CA) to a steel pin held in the head of the wear tester. This head was mounted on a custom-made load cell that measured the horizontal force produced by friction between the specimen pin and the abrasive turntable. During the wear test, the specimen pin was held against the turntable with adjustable weights that, for the standard test, produced a downward force of 0.019 Newtons. The surface of the turntable was covered with 600 grit abrasive paper (#413Q, 3 M Corporation, www.3M.com). The turntable rotation period was usually set to about 4 s, resulting in an interaction velocity of 0.027 m/s.The volume worn away was calculated from “before” and “after” measurements of microscope images (image J software) taken from the side (e.g. Fig. 2B) to measure the height of the approximate cone of worn material, and from face-on in order to measure the area of the base and top of the frustum of worn-away material. The horizontal force was recorded continuously during the wear period. The wear rate (w), defined as the volume worn away per unit energy expended, was approximated as follows:$$w= frac{V}{Fd} ,$$where F is the average force of friction measured during the wear period by the load cell, d is the distance traveled by the pin over the abrasive paper and “V” is the worn volume, approximated as a frustum:$$V = 1{/}3,Delta L , (A1 + left( {A1*A2} right)^{1/2} + A2)$$where “A1” and “A2” are the areas of the worn surface before and after the wear sessions (the area of any voids or internal lumens was subtracted from the area of the cross sections) and ΔL, the change in the length of the specimen due to wear. We defined wear resistance as the inverse of the wear rate, 1/w. While we expect this test to be most useful for relative comparisons, and the value is expected to vary somewhat with abrasive properties and normal forces, we found no statistically distinguishable difference in values from 4 samples that were re-run using a ten-times greater force to press them against the abrasive paper14.Molecular composition and nanometer-scale structureIn order to better understand the composition and structure of the HEBs—down to an atomic scale—we examined a representative HEB using Atom Probe Tomography (APT). We checked the APT results and studied their generality using Time-of-Flight–Secondary Ion Mass Spectrometry (ToF–SIMS). Both of these techniques use a pulsed beam (laser and ion respectively) to break the specimen into molecular fragments that are accelerated to a detector; for a particular charge, heavier fragments travel more slowly and arrive later at the detector. The arrival time differences are used to identify the fragments by their mass, giving information about, for example, the atoms attached to zinc atoms in the specimen and, from APT, the spatial distribution of zinc atoms on a nanometer scale.Atom probe tomography (APT)APT is a 3D nanoscale characterization method in which field evaporated ions from a sharpened needle specimen are analyzed by a position-sensitive single-particle detector, in order to provide an isotopically resolved three-dimensional representation of the real-space specimen elemental distribution59. The field evaporation of non-conductive samples is achieved using a pulsed laser focused on the needle specimen apex.A FIB-SEM based lift-out procedure was used to prepare needle-shaped APT specimens using FEI Helios 600i at the University of Oregon CAMCOR facility, and a Helios Dual Beam Nanolab 600 FIB-SEM housed at Environmental Molecular Sciences Laboratory, PNNL.The APT analysis was carried out using a CAMECA LEAP (local electrode atom probe) 4000X HR system equipped with a 355 nm wavelength picosecond pulsed UV laser. A 30 K sample base temperature and a 100 or 200 kHz laser pulse repetition rate was used. Atom probe data reconstruction and analysis was performed using Cameca IVAS software.
Development of APT techniques for these organic materials.APT has not typically been used to examine organic materials, so we began by examining standards (such as zinc picolinate) and adjusting beam current densities and other parameters in both the FIB-SEM preparation of APT samples and in APT itself in order to minimize physical damage detected with SEM and to minimize differences between the chemical formulae and APT results for standards60. We found that current densities often employed in FIB-SEM milling were much too high for our organic specimens, resulting in beam damage visible in SEM.Based on the standards and SEM evidence of damage, we used FIB-SEM currents for producing the sample needles, and APT laser pulses for promoting evaporation, that were similar to or smaller than those used in other investigations of organic materials61,62,63,64,65,66,67,68. We used ≤ 21 pA for the electron beam, ≤ 80 pA for ion beam “trenching”, 7.7 pA for ion beam imaging and for cutting the cantilever “liftout” piece, and ≤ 24 pA for sharpening needles. The results reported here are based on samples analyzed using 10, 20 or 100 pJ laser pulses.Identification of molecular fragments from mass-to-charge ratios is particularly difficult for organic materials because of the many possible fragments of similar mass, and several techniques have been developed to aid in this analysis61,62,63,64,68,69,70,71,72,73.Our identification of zinc-containing fragments was simplified by the pattern of the 3 or 4 main zinc isotopes (Fig. 6A). In addition, we used resources with lists of fragments as a function of mass (e.g. https://webbook.nist.gov/chemistry/mw-ser/). We also cross-checked fragment identification using a Time-of-Flight–Secondary Ion Mass Spectrometry (ToF–SIMS) system.Time of flight–secondary ion mass spectrometry (ToF–SIMS)We used a ToF–SIMS system that had a higher mass-to-charge resolution than the APT system (although it had micron- instead of nanometer-scale spatial resolution) in order to check APT fragment identification. The higher mass sensitivity of the ToF–SIMS system provided additional evidence that, for example, the fragments identified as ZnCN were not actually ZnC2H2, which is only about 0.02% lighter. We also used ToF–SIMS to study the larger-scale spatial distribution of fragments, and similarities with HEBs from other species.We used an ION-TOF ToF–SIMS IV, manufactured by ION-TOF GmbH, Muenster, Germany. The primary ion beam was Bi3+ (25 kV, 10 kHz, 0.4 pA); the static limit (2 × 1012 ions/cm2) was not exceeded. The dimensions of the analysis area varied, but were between 100 × 100 μm and 300 × 300 μm. A low energy electron beam was used for charge neutralization. The spectra were analyzed using the vendor’s software. Chemical maps of peaks of interest were created from the total spectra and used as a basis for retrospective analysis—i.e., pixel-specific extraction of spectra in order to determine the chemical makeup of features of interest. More

1.Sivapalan, M., Savenije, H. H. & Blöschl, G. Socio-hydrology: A new science of people and water. Hydrol. Process. 26(8), 1270–1276 (2012).ADS 
Article 

Google Scholar 
2.Babel, M. S., Das Gupta, A. & Nayak, D. K. A model for optimal allocation of water to competing demands. Water Resour. Manag. 19(6), 693–712 (2005).Article 

Google Scholar 
3.Banihabib, M. E., Zahraei, A. & Eslamian, S. An integrated optimisation model of reservoir and irrigation system applying uniform deficit irrigation. Int. J. Hydrol. Sci. Technol. 5(4), 372–385 (2015).Article 

Google Scholar 
4.Ghahreman, B. & Sepaskhah, A. R. Optimal water allocation of water from a single reservoir to an irrigation project with pre-determined multiple cropping patterns. Irrig. Sci. 21(3), 127–137 (2002).Article 

Google Scholar 
5.Xevi, E. & Khan, S. A multi-objective optimisation approach to water management. J. Environ. Manag. 77(4), 269–277 (2005).CAS 
Article 

Google Scholar 
6.Divakar, L., Babel, M. S., Perret, S. R. & Das Gupta, A. Optimal allocation of bulk water supplies to competing use sectors based on economic criterion—An application to the Chao Phraya River Basin, Thailand. J. Hydrol. 401(1–2), 22–35 (2011).ADS 
Article 

Google Scholar 
7.Roozbahani, R., Abbasi, B., Schreider, S. & Ardakani, A. A multi-objective approach for transboundary river water allocation. Water Resour. Manag. 28(15), 5447–5463 (2014).Article 

Google Scholar 
8.Schnegg, M. & Kiaka, R. D. The economic value of water: The contradictions and consequences of a prominent development model in Namibia. Econ. Anthropol. 6(2), 264–276 (2019).
Google Scholar 
9.Langarudi, S. P., Maxwell, C. M., Bai, Y., Hanson, A. & Fernald, A. Does socioeconomic feedback matter for water models?. Ecol. Econ. 159, 35–45 (2019).Article 

Google Scholar 
10.Keshavarz, M., Karami, E. & Vanclay, F. The social experience of drought in rural Iran. Land Use Policy 30(1), 120–129 (2013).Article 

Google Scholar 
11.Dean, A. J., Fielding, K. S., Lindsay, J., Newton, F. J. & Ross, H. How social capital influences community support for alternative water sources. Sustain. Cities Soc. 27, 457–466 (2016).Article 

Google Scholar 
12.Scanlon, T. et al. The role of social actors in water access in Sub-Saharan Africa: Evidence from Malawi and Zambia. Water Resour. Rural Dev. 8, 25–36 (2016).Article 

Google Scholar 
13.Popovic, T., Kraslawski, A., Heiduschke, R. & Repke. J. Indicators of social sustainability for wastewater treatment processes. In Computer Aided Chemical Engineering, Vol. 723(28) (2014).14.El-Gafy, I. K. E. D. The water poverty index as an assistant tool for drawing strategies of the Egyptian water sector. Ain Shams Eng. J. 9(2), 173–186 (2018).Article 

Google Scholar 
15.Bui, N. T. et al. Social sustainability assessment of groundwater resources: A case study of Hanoi, Vietnam. Ecol. Indic. 93, 1034–1042 (2018).Article 

Google Scholar 
16.Li, C. et al. Three decades of changes in water environment of a large freshwater lake and its relationship with socio-economic indicators. J. Environ. Sci. 77, 156–166 (2019).Article 

Google Scholar 
17.Ahmadi, A., Karamouz, M., Moridi, A. & Han, D. Integrated planning of land use and water allocation on a watershed scale considering social and water quality issues. J. Water Resour. Plan. Manag. 138(6), 671–681 (2012).Article 

Google Scholar 
18.Tu, Y. et al. Administrative and market-based allocation mechanism for regional water resources planning. Resour. Conserv. Recycl. 95, 156–173 (2015).MathSciNet 
Article 

Google Scholar 
19.Kelly, R. A. et al. Selecting among five common modelling approaches for integrated environmental assessment and management. Environ. Model. Softw. 47, 159–181 (2013).Article 

Google Scholar 
20.Wang, K., Davies, E. G. & Liu, J. Integrated water resources management and modeling: A case study of Bow river basin, Canada. J. Clean. Prod. 240, 118242 (2019).Article 

Google Scholar 
21.Zhao, D. et al. Quantifying economic-social-environmental trade-offs and synergies of water-supply constraints: An application to the capital region of China. Water Res. 195, 116986 (2021).CAS 
Article 

Google Scholar 
22.Navarro-Ramírez, V., Ramírez-Hernandez, J., Gil-Samaniego, M. & Rodríguez-Burgueño, J. E. Methodological frameworks to assess sustainable water resources management in industry: A review. Ecol. Indic. 119, 106819 (2020).Article 

Google Scholar 
23.Iran Environment Organization. https://www.doe.ir/portal/home (2013).24.Madani, K. & Mariño, M. A. System dynamics analysis for managing Iran’s Zayandeh-Rud river basin. Water Resour. Manag. 23(11), 2163–2187 (2009).Article 

Google Scholar 
25.Ahmad, S. & Simonovic, S. P. System dynamics modeling of reservoir operations for flood management. J. Comput. Civ. Eng. 14(3), 190–198 (2000).Article 

Google Scholar 
26.Ahmad, S. & Simonovic, S. P. Spatial system dynamics: New approach for simulation of water resources systems. J. Comput. Civ. Eng. 18(4), 331–340 (2004).Article 

Google Scholar 
27.Ahmad, S. & Simonovic, S. P. An intelligent decision support system for management of floods. Water Resour. Manag. 20(3), 391–410 (2006).Article 

Google Scholar 
28.Dong, Q., Zhang, X., Chen, Y. & Fang, D. Dynamic management of a water resources-socioeconomic-environmental system based on feedbacks using system dynamics. Water Resour. Manag. 33(6), 2093–2108 (2019).Article 

Google Scholar 
29.National Statistical Center of Iran. Summary of Industrial Workshop Statistics by Activity. http://www.amar.org.ir (2013).30.Rezaee, A., Bozorg-Haddad, O. & Singh, V. P. Water and society. in Economical, Political, and Social Issues in Water Resources, 257–271 (Elsevier, 2021).Chapter 

Google Scholar 
31.Ministry of Energy, Water and Wastewater Macro Planning Office. Studies on Updating the Country’s Comprehensive Water Plan in Aras, Urmia, Talesh-Anzali Wetland, Sefidrood-Haraz, Haraz-Qarasu, Gorganrood and Atrak Watersheds (Consumption and Needs of Industrial and Mining Water and Production Wastewater in the Base Year (2006)) in the Catchment Area of Urmia), Vol. 8 (2013).32.Hwang, C. L. & Yoon, K. Methods for multiple attribute decision making. in Multiple Attribute Decision Making: Lecture Notes in Economics and Mathematical Systems (Springer, 1981).Chapter 

Google Scholar 
33.Zeyaeyan, S., Fattahi, E., Ranjbar, A. & Vazifedoust, M. Classification of rainfall warnings based on the TOPSIS method. Climate 5(2), 33 (2017).Article 

Google Scholar 
34.Zolghadr-Asli, B., Bozorg-Haddad, O., Enayati, M. & Goharian, E. Developing a robust multi-attribute decision-making framework to evaluate performance of water system design and planning under climate change. Water Resour. Manag. 35(1), 279–298 (2021).Article 

Google Scholar 
35.Opricovic, S. & Tzeng, G. H. Compromise solution by MCDM methods: A comparative analysis of VIKOR and TOPSIS. Eur. J. Oper. Res. 156(2), 445–455 (2004).Article 

Google Scholar 
36.Shannon, C. E. A mathematical theory of communication. Bell Syst. Tech. J. 27(3), 379–423 (1948).MathSciNet 
Article 

Google Scholar  More

1.Brevik, E. C. et al. The interdisciplinary nature of SOIL. Soil 1(1), 117–129. https://doi.org/10.5194/soil-1-117-2015 (2015).Article 

Google Scholar 
2.Liu, X. et al. Heavy metal concentrations of soils near the large opencast coal mine pits in China. Chemosphere 244, 125360. https://doi.org/10.1016/j.chemosphere.2019.125360 (2020).ADS 
CAS 
Article 
PubMed 

Google Scholar 
3.Imin, B., Abliz, A., Shi, Q., Liu, S. & Hao, L. Quantitatively assessing the risks and possible sources of toxic metals in soil from an arid, coal-dependent industrial region in NW China. J. Geochem. Explor. https://doi.org/10.1016/j.gexplo.2020.106505 (2020).Article 

Google Scholar 
4.Doran, J. W. & Parkin, T. B. Defining and assessing soil quality. Defin. Soil Qual. Sustain. Environ. 35, 1–21. https://doi.org/10.2136/sssaspecpub35.c1 (1994).Article 

Google Scholar 
5.Sun, H. et al. Effects of soil quality on effective ingredients of Astragalus mongholicus from the main cultivation regions in China. Ecol. Indic. 114, 106296. https://doi.org/10.1016/j.ecolind.2020.106296 (2020).CAS 
Article 

Google Scholar 
6.Alloway, B. J. Sources of Heavy Metals and Metalloids in Soils. Heavy Metals in Soils 11–50 (Springer, 2013).Book 

Google Scholar 
7.Yang, Q. Q. et al. A review of soil heavy metal pollution from industrial and agricultural regions in China: Pollution and risk assessment. Sci. Total Environ. 642, 690–700. https://doi.org/10.1016/j.scitotenv.2018.06.068 (2018).ADS 
CAS 
Article 
PubMed 

Google Scholar 
8.Huang, Y., Kuang, X., Cao, Y. & Bai, Z. The soil chemical properties of reclaimed land in an arid grassland dump in an opencast mining area in China. RSC Adv. 8(72), 41499–41508. https://doi.org/10.1039/c8ra08002j (2018).ADS 
CAS 
Article 

Google Scholar 
9.Liu, Z. J. et al. Soil quality assessment of Albic soils with different productivities for eastern China. Soil Till. Res. 140, 74–81. https://doi.org/10.1016/j.still.2014.02.010 (2014).Article 

Google Scholar 
10.Bhardwaj, A. K., Jasrotia, P., Hamilton, S. K. & Robertson, G. P. Ecological management of intensively cropped agro-ecosystems improves soil quality with sustained productivity. Agr. Ecosyst. Environ. 140(3–4), 419–429. https://doi.org/10.1016/j.agee.2011.01.005 (2011).Article 

Google Scholar 
11.Mendham, D. S. et al. Soil analyses as indicators of phosphorus response in young eucalypt plantations. Soil Sci. Soc. Am. J. 66(3), 959–968. https://doi.org/10.2136/sssaj2002.9590 (2002).ADS 
CAS 
Article 

Google Scholar 
12.Shukla, M. K., Lal, R. & Ebinger, M. Determining soil quality indicators by factor analysis. Soil Till. Res. 87(2), 194–204. https://doi.org/10.1016/j.still.2005.03.011 (2006).Article 

Google Scholar 
13.Vasu, D. et al. Soil quality index (SQI) as a tool to evaluate crop productivity in semi-arid Deccan plateau. India. Geoderma. 282, 70–79. https://doi.org/10.1016/j.geoderma.2016.07.010 (2016).ADS 
CAS 
Article 

Google Scholar 
14.Mishra, G. et al. Soil quality assessment under shifting cultivation and forests in Northeastern Himalaya of India. Arch. Agron. Soil Sci. 63(10), 1355–1368. https://doi.org/10.1080/03650340.2017.1281390 (2017).CAS 
Article 

Google Scholar 
15.Li, X. Y., Wang, D. Y., Ren, Y. X., Wang, Z. M. & Zhou, Y. H. Soil quality assessment of croplands in the black soil zone of Jilin Province, China: Establishing a minimum data set model. Ecol. Indic. 107, 105251. https://doi.org/10.1016/j.ecolind.2019.03.028 (2019).CAS 
Article 

Google Scholar 
16.Zhao, Q. Q. et al. Effects of freshwater inputs on soil quality in the Yellow River Delta. China. Ecol. Indic. 98, 619–626. https://doi.org/10.1016/j.ecolind.2018.11.041 (2019).CAS 
Article 

Google Scholar 
17.Li, F. P., Liu, W., Lu, Z. B., Mao, L. C. & Xiao, Y. H. A multi-criteria evaluation system for arable land resource assessment. Environ. Monit. Assess. https://doi.org/10.1007/s10661-019-8023-x (2020).Article 
PubMed 
PubMed Central 

Google Scholar 
18.Raiesi, F. A minimum data set and soil quality index to quantify the effect of land use conversion on soil quality and degradation in native rangelands of upland arid and semiarid regions. Ecol. Indic. 75, 307–320. https://doi.org/10.1016/j.ecolind.2016.12.049 (2017).Article 

Google Scholar 
19.Zhou, Y. et al. Assessment of soil quality indexes for different land use types in typical steppe in the loess hilly area, China. Ecol. Indic. 118, 106743. https://doi.org/10.1016/j.ecolind.2020.106743 (2020).CAS 
Article 

Google Scholar 
20.Cheng, W. et al. Geographic distribution of heavy metals and identification of their sources in soils near large, open-pit coal mines using positive matrix factorization. J. Hazard. Mater. 387, 121666. https://doi.org/10.1016/j.jhazmat.2019.121666 (2020).CAS 
Article 
PubMed 

Google Scholar 
21.Zhao, X., Tong, M., He, Y., Han, X. & Wang, L. A comprehensive, locally adapted soil quality indexing under different land uses in a typical watershed of the eastern Qinghai-Tibet Plateau. Ecol. Ind. 125, 107445. https://doi.org/10.1016/j.ecolind.2021.107445 (2021).CAS 
Article 

Google Scholar 
22.Zhang, W. S. et al. Comprehensive assessment methodology of soil quality under different land use conditions. Trans. Chin. Soc. Agric. Eng. 26(12), 311–318. https://doi.org/10.3969/j.issn.1002-6819.2010.12.053 (2010).Article 

Google Scholar 
23.Batjargal, T., Otgonjargal, E., Baek, K. & Yang, J. S. Assessment of metals contamination of soils in Ulaanbaatar, Mongolia. J. Hazard. Mater. 184(1–3), 872–876. https://doi.org/10.1016/j.jhazmat.2010.08.106 (2010).CAS 
Article 
PubMed 

Google Scholar 
24.Ngole-Jeme, V. M. Heavy metals in soils along unpaved roads in south west Cameroon: Contamination levels and health risks. Ambio 45(3), 374–386. https://doi.org/10.1007/s13280-015-0726-9 (2016).CAS 
Article 
PubMed 

Google Scholar 
25.China Soil Census Office. China Soil Census Data[M] (China National Agricultural Press, Beijing, 1997).26.Chen, H., Teng, Y., Lu, S., Wang, Y. & Wang, J. Contamination features and health risk of soil heavy metals in China. Sci. Total Environ. 512, 143–153. https://doi.org/10.1016/j.scitotenv.2015.01.025 (2015).ADS 
CAS 
Article 
PubMed 

Google Scholar 
27.Wang, Y., Duan, X. & Wang, L. Spatial distribution and source analysis of heavy metals in soils influenced by industrial enterprise distribution: Case study in Jiangsu Province. Sci. Total Environ. 710, 134953. https://doi.org/10.1016/j.scitotenv.2019.134953 (2020).ADS 
CAS 
Article 
PubMed 

Google Scholar 
28.Bao, S. D. Soil Agrochemical Analysis 25–114 (China Agricultural Press, 2000).
Google Scholar 
29.Wang, M. E., Peng, C., & Chen, W. P. Impacts of industrial zone in arid area in Ningxia province on the accumulation of heavy metals in agricultural soils. Chin. J. Envir. Sci., 37(9), 3532–3539. https://doi.org/10.13227/j.hjkx.2016.09.035 (2016). Article 

Google Scholar 
30.Xu, Z. et al. Characteristics and sources of heavy metal pollution in desert steppe soil related to transportation and industrial activities. Environ. Sci. Pollut. Res. 27, 38835–38848. https://doi.org/10.1007/s11356-020-09877-9 (2020).CAS 
Article 

Google Scholar 
31.Qi, Y. B. et al. Evaluating soil quality indices in an agricultural region of Jiangsu Province. China. Geoderma. 149(3–4), 325–334. https://doi.org/10.1016/j.geoderma.2008.12.015 (2009).ADS 
Article 

Google Scholar 
32.Hu, Q., Chen, W. F., Song, X. L., Dong, Y. J. & Liu, Z. Q. Effects of reclamation/cultivation on soil quality of Saline-alkali Soils in the yellow river delta. Acta Pedol. Sin. 57(4), 824–833. https://doi.org/10.11766/trxb201905050105 (2020).Article 

Google Scholar 
33.Qu, X. G., Sun, Y. X. & Fu, X. Y. Soil quality and stripping depth evaluation of tillage layer for construction of Qingdao new airport. Bull. Soil Water Conserv. 38(4), 202–206. https://doi.org/10.13961/j.cnki.stbctb.2018.04.033 (2018).Article 

Google Scholar 
34.Abd-Elwahed, M. S. Influence of long-term wastewater irrigation on soil quality and its spatial distribution. Ann. Agric. Sci. 63(2), 191–199. https://doi.org/10.1016/j.aoas.2018.11.004 (2018).Article 

Google Scholar 
35.CNEMC (China National Environmental Monitoring Center). The Background Values of Elements in Chinese Soils. 330–493 (Environmental Science Press of China, 1990).36.Cheng, J. L., Shi, Z., Zhu, Y. W., Liu, C. & Li, H. Y. Differential characteristics and appraisal of heavy metals in agricultural soils of Zhejiang Province. J. Soil Water Conserv. 20(1), 103–107. https://doi.org/10.1016/S1872-2032(06)60052-8 (2006).Article 

Google Scholar 
37.Jin, G. Q. et al. Source apportionment of heavy metals in farmland soil with application of APCS-MLR model: A pilot study for restoration of farmland in Shaoxing City Zhejiang. China. Ecotox. Environ. Safe. 184, 109495. https://doi.org/10.1016/j.ecoenv.2019.109495 (2019).CAS 
Article 

Google Scholar 
38.Marzaioli, R., D’Ascoli, R., De Pascale, R. A. & Rutigliano, F. A. Soil quality in a Mediterranean area of Southern Italy as related to different land use types. Appl. Soil Ecol. 44(3), 205–212. https://doi.org/10.1016/j.apsoil.2009.12.007 (2010).Article 

Google Scholar 
39.Zhao, N., Meng, P., Zhang, J. S., Lu, S. & Cheng, Z. Q. Soil quality assessment of Robinia psedudoacia plantations with various ages in the Grain-for-Green Program in hilly area of North China. Yingyong Shengtai Xuebao https://doi.org/10.13287/j.1001-9332.2014.0038 (2014).Article 
PubMed 

Google Scholar 
40.Zheng, Q. et al. Comprehensive method for evaluating soil quality in cotton fields in Xinjiang. China. Chin. J. Appl. Ecol. 29(4), 1291–1301. https://doi.org/10.13287/j.1001-9332.201804.029 (2018).Article 

Google Scholar 
41.Turrión, M. B. et al. Soil phosphorus forms as quality indicators of soils under different vegetation covers. Sci. Total Environ. 378(1–2), 195–198. https://doi.org/10.1016/j.scitotenv.2007.01.037 (2007).ADS 
CAS 
Article 
PubMed 

Google Scholar 
42.Barbosa, E. R. M. et al. Short-term effect of nutrient availability and rainfall distribution on biomass production and leaf nutrient content of Savanna tree species. PLoS ONE 9(3), e92619. https://doi.org/10.1371/journal.pone.0092619 (2014).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
43.Marty, C., Houle, D., Gagnon, C. & Courchesne, F. The relationships of soil total nitrogen concentrations, pools and C: N ratios with climate, vegetation types and nitrate deposition in temperate and boreal forests of eastern Canada. CATENA 152, 163–172. https://doi.org/10.1016/j.catena.2017.01.014 (2017).CAS 
Article 

Google Scholar 
44.Chen, Z. F. et al. Evaluation on cultivated-layer soil quality of sloping farmland in Yunnan based on soil management assessment framework (SMAF). Trans. Chin. Soc. Agric. Eng. 35(03), 256–267. https://doi.org/10.11975/j.issn.1002-6819.2019.03.032 (2019).Article 

Google Scholar 
45.Ding, J. X. et al. Spatial distribution of the herbaceous layer and its relationship to soil physical–chemical properties in the southern margin of the Gurbantonggut Desert, northwestern China. Acta Ecol. Sin. 36(5), 327–332. https://doi.org/10.1016/j.chnaes.2016.06.006 (2016).Article 

Google Scholar 
46.Güntner, A., Seibert, J. & Uhlenbrook, S. Modeling spatial patterns of saturated areas: An evaluation of different terrain indices. Water Resour. Res. https://doi.org/10.1029/2003wr002864 (2004).Article 

Google Scholar 
47.Yenilmez, F., Kuter, N., Emil, M. K. & Aksoy, A. Evaluation of pollution levels at an abandoned coal mine site in Turkey with the aid of GIS. Int. J. Coal Geol. 86(1), 12–19. https://doi.org/10.1016/j.coal.2010.11.012 (2011).CAS 
Article 

Google Scholar 
48.Kronbauer, M. A. et al. Geochemistry of ultra-fine and nano-compounds in coal gasification ashes: A synoptic view. Sci. Total Environ. https://doi.org/10.1016/j.scitotenv.2013.02.066 (2013).Article 
PubMed 

Google Scholar 
49.Masto, R. E. et al. Assessment of environmental soil quality around Sonepur Bazari mine of Raniganj coalfield, India. Solid. Earth. 6(3), 811. https://doi.org/10.5194/se-6-811-2015 (2015).ADS 
Article 

Google Scholar 
50.Han, Y. et al. Effects of opencast coal mining on soil properties and plant communities of grassland. Chin. J. Ecol. 38(11), 3425–3422. https://doi.org/10.13292/j.1000-4890.201911.011 (2019).Article 

Google Scholar 
51.Liu, J., Wu, L. C., Chen, D., Li, M. & Wei, C. J. Soil quality assessment of different Camellia oleifera stands in mid-subtropical China. Appl. Soil Ecol. 113, 29–35. https://doi.org/10.1016/j.apsoil.2017.01.010 (2017).ADS 
Article 

Google Scholar 
52.Yu, P. J., Liu, S. W., Zhang, L., Li, Q. & Zhou, D. W. Selecting the minimum data set and quantitative soil quality indexing of alkaline soils under different land uses in northeastern China. Sci. Total Environ. 616–617, 564–571. https://doi.org/10.1016/j.scitotenv.2017.10.301 (2018).ADS 
CAS 
Article 
PubMed 

Google Scholar 
53.Liu, Q. Q., Zhang, T., Wang, C. & Liu, J. H. Comparison of vegetation composition and soil fertility quality inside and outside the wind farm. J. Inner Mongolia Agric. Univ. (nat. Sci. Edn.) 41(02), 30–36. https://doi.org/10.16853/j.cnki.1009-3575.2020.02.006 (2020).CAS 
Article 

Google Scholar 
54.Sheldrick, W., Syers, J. K. & Lingard, J. Contribution of livestock excreta to nutrient balances. Nutr. Cycl. Agroecosys. 66(2), 119–131. https://doi.org/10.1023/a:1023944131188 (2003).Article 

Google Scholar 
55.Kasahara, M., Fujii, S., Tanikawa, T. & Mori, A. S. Ungulates decelerate litter decomposition by altering litter quality above and below ground. Eur. J. Forest Res. 135(5), 849–856. https://doi.org/10.1007/s10342-016-0978-3 (2016).Article 

Google Scholar 
56.Zhan, T. Y. et al. Meta-analysis demonstrating that moderate grazing can improve the soil quality across China’s grassland ecosystems. Appl. Soil Ecol. 147, 103438. https://doi.org/10.1016/j.apsoil.2019.103438 (2020).Article 

Google Scholar 
57.Liu, X. Y., Bai, Z. K., Zhou, W., Cao, Y. G. & Zhang, G. J. Changes in soil properties in the soil profile after mining and reclamation in an opencast coal mine on the Loess Plateau. China. Ecol. Eng. 98, 228–239. https://doi.org/10.1016/j.ecoleng.2016.10.078 (2017).Article 

Google Scholar 
58.Sun, L. et al. Levels, sources, and spatial distribution of heavy metals in soils from a typical coal industrial city of Tangshan, China. CATENA 175, 101–109. https://doi.org/10.1016/j.catena.2018.12.014 (2019).CAS 
Article 

Google Scholar 
59.Yang, S. L., Zhou, D. Q., Yu, H. Y., Wei, R. & Pan, B. Distribution and speciation of metals (Cu, Zn, Cd, and Pb) in agricultural and non-agricultural soils near a stream upriver from the Pearl River. China. Environ. Pollut. 177, 64–70. https://doi.org/10.1016/j.envpol.2013.01.044 (2013).CAS 
Article 
PubMed 

Google Scholar 
60.Zhao, F. J., Ma, Y., Zhu, Y. G., Tang, Z. & McGrath, S. P. Soil Contamination in China: Current Status and Mitigation Strategies. Environ. Sci. Technol. 49(2), 750–759. https://doi.org/10.1021/es5047099 (2014).ADS 
CAS 
Article 

Google Scholar 
61.Wang, Y. Z., Duan, X. J. & Wang, L. Spatial distribution and source analysis of heavy metals in soils influenced by industrial enterprise distribution: Case study in Jiangsu Province. Sci. Total Environ. https://doi.org/10.1016/j.scitotenv.2019.134953 (2019).Article 
PubMed 
PubMed Central 

Google Scholar 
62.Nehrani, S. H. et al. Quantification of soil quality under semi-arid agriculture in the northwest of Iran. Ecol. Indic. 108, 105770. https://doi.org/10.1016/j.ecolind.2019.105770 (2020).CAS 
Article 

Google Scholar 
63.Huang, Y. et al. Heavy metal pollution and health risk assessment of agricultural soils in a typical peri-urban area in southeast China. J. Environ. Manage. 207, 159–168. https://doi.org/10.1016/j.jenvman.2017.10.072 (2018).CAS 
Article 
PubMed 

Google Scholar 
64.Qu, C. S. et al. Spatial distribution, risk and potential sources of lead in soils in the vicinity of a historic industrial site. Chemosphere 205, 244–252. https://doi.org/10.1016/j.chemosphere.2018.04.119 (2018).ADS 
CAS 
Article 
PubMed 

Google Scholar 
65.Charlesworth, S., Everett, M., McCarthy, R., Ordóñez, A. & de Miguel, E. A comparative study of heavy metal concentration and distribution in deposited street dusts in a large and a small urban area: Birmingham and Coventry, West Midlands, UK. Environ. Int. 29(5), 563–573. https://doi.org/10.1016/s0160-4120(03)00015-1 (2003).CAS 
Article 
PubMed 

Google Scholar 
66.Liang, J. et al. Facile synthesis of alumina-decorated multi-walled carbon nanotubes for simultaneous adsorption of cadmium ion and trichloroethylene. Chem. Eng. J. 273, 101–110. https://doi.org/10.1016/j.cej.2015.03.069 (2015).CAS 
Article 

Google Scholar 
67.Liang, J. et al. Spatial distribution and source identification of heavy metals in surface soils in a typical coal mine city, Lianyuan. China. Environ. Pollut. 225, 681–690. https://doi.org/10.1016/j.envpol.2017.03.057 (2017).CAS 
Article 
PubMed 

Google Scholar 
68.Chen, H., Lu, X. W., Li, L. Y., Gao, T. N. & Chang, Y. Y. Metal contamination in campus dust of Xi’an, China: A study based on multivariate statistics and spatial distribution. Sci. Total. Environ. 484, 27–35. https://doi.org/10.1016/j.scitotenv.2014.03.026 (2014).ADS 
CAS 
Article 
PubMed 

Google Scholar 
69.Adachi, K. & Tainosho, Y. Characterization of heavy metal particles embedded in tire dust. Environ. Int. 30(8), 1009–1017. https://doi.org/10.1016/j.envint.2004.04.004 (2004).CAS 
Article 
PubMed 

Google Scholar 
70.Garcia-Guinea, J. et al. Influence of accumulation of heaps of steel slag on the environment: Determination of heavy metals content in the soils. An. Acad. Bras. Cienc. 82(2), 267–277. https://doi.org/10.1590/S0001-37652010000200003 (2010).CAS 
Article 
PubMed 

Google Scholar 
71.Fan, X. G., Mi, W. B., Ma, Z. N. & Wang, T. Y. Spatial and temporal characteristics of heavy metal concentration of surface soil in Hebin industrial park in Shizuishan northwest China. Chin. J. Envir. Sci. 34(5), 1887–1894. https://doi.org/10.13227/j.hjkx.2013.05.033 (2013).Article 

Google Scholar 
72.Huang, T., Yue, X. J., Ge, X. Z. & Wang, X. D. Evaluation of soil quality on gully region of loess plateau based on principal component analysis. Agri. Res. Arid Areas. 28(03), 141–147. https://doi.org/10.1016/S1002-0160(10)60014-8 (2010).Article 

Google Scholar 
73.Jiang, L. B. et al. Co-pelletization of sewage sludge and biomass: The density and hardness of pellet. Bioresour. Technol. 166, 435–443. https://doi.org/10.1016/j.biortech.2014.05.077 (2014).CAS 
Article 
PubMed 

Google Scholar 
74.Oumenskou, H. et al. Multivariate statistical analysis for spatial evaluation of physicochemical properties of agricultural soils from Beni-Amir irrigated perimeter, Tadla plain, Morocco. Geol. Ecol. Landsc. 3(2), 83–94 (2019).Article 

Google Scholar 
75.Liu, Y., Wang, L., Liu, B. H. & Henderson, M. Observed changes in shallow soil temperatures in Northeast China, 1960–2007. Clim. Res. 67(1), 31–42. https://doi.org/10.3354/cr01351 (2016).Article 

Google Scholar 
76.Jiang, Y. F. et al. Distribution, compositional pattern and sources of polycyclic aromatic hydrocarbons in urban soils of an industrial city, Lanzhou. China. Ecotox. Environ. Safe. 126, 154–162. https://doi.org/10.1016/j.ecoenv.2015.12.037 (2016).CAS 
Article 

Google Scholar 
77.Frohne, T. & Rinklebe, J. Biogeochemical fractions of mercury in soil profiles of two different floodplain ecosystems in Germany. Water Air Soil Poll. 224(6), 1591. https://doi.org/10.1007/s11270-013-1591-4 (2013).ADS 
CAS 
Article 

Google Scholar 
78.Stefanowicz, A. M., Kapusta, P., Zubek, S., Stanek, M. & Woch, M. W. Soil organic matter prevails over heavy metal pollution and vegetation as a factor shaping soil microbial communities at historical Zn–Pb mining sites. Chemosphere 240, 124922. https://doi.org/10.1016/j.chemosphere.2019.124922 (2020).ADS 
CAS 
Article 
PubMed 

Google Scholar  More

OverviewOur study investigates the patterns and drivers of global marine plankton diversity by simultaneously modeling the spatial distribution of 860 phyto- and zooplankton species, based on the widest and most recent compilations of in situ observations available. These observations were associated with various sets of relevant predictors to train a range of statistical species distribution models (SDMs) on a monthly resolution. The SDMs were used to estimate contemporary and future levels of global surface species richness (SR) for total plankton, phytoplankton and zooplankton. We explore how, and why, global phyto- and zooplankton SR and community composition are affected by future climate change under the RCP8.5 scenario of greenhouse gas (GHG) emissions. We also summarize regional patterns of climate change impacts on plankton diversity by clustering the global ocean and examine how hotspots of climate change impacts might overlap with the current provision of marine ecosystem services. All data manipulation and analyses were performed under the R programming language39. The R packages used are mentioned below in their corresponding section.Plankton species observationsFirst, to model global, open ocean plankton diversity from species-level field observations, comparable datasets of phytoplankton and zooplankton occurrences (i.e., presences) had to be compiled. We refer to as open ocean all those regions where the seafloor depth exceeds 200 m. We made use of the large dataset of phytoplankton occurrences recently compiled by Righetti et al.11. For zooplankton, a new dataset was compiled following the same methodology. Both occurrence datasets were based on publicly available data from online biodiversity repositories, as well as some additional published datasets. The R packages mainly used for implementing the datasets are those constituting the tidyverse package.Phytoplankton occurrences
For the phytoplankton occurrences used here, Righetti et al.40 compiled data from various sources: the Global Biodiversity Information Facility (GBIF; https://www.gbif.org), the Ocean Biogeographic Information System (OBIS; https://www.obis.org), the data from Villar et al.41, and the MAREDAT initiative42. Righetti et al.40 gathered >106 presences from nearly 1300 species sampled through various methodologies within the monthly climatological mixed-layer depth, at an average depth of 5.41 ± 6.95 m (mean ± sd), between 1800 and 2015. The species names were corrected and harmonized following the reference list of Algaebase (http://www.algaebase.org/) and were further validated by expert opinion. The final species list spanned most of the extant phytoplankton taxa composing the biodiversity of the euphotic zone of the global ocean. Fossil records, sedimentary records, and occurrences associated with senseless metadata were removed. This dataset has been mined to effectively obtain phytoplankton SR estimates for the global open ocean that were: (i) robust to sampling spatial-temporal biases11, and (ii) validated against independent data11.

Zooplankton occurrences
A new dataset of global zooplankton species occurrences was compiled in a comparable fashion to that put together for phytoplankton. Prior to retrieving the occurrence data online, we first identified the phyla (Order/Class/Family) that comprise the bulk of extant oceanic zooplankton communities: Copelata (i.e., appendicularians), Ctenophora, Cubozoa (i.e., box jellyfish), Euphausiidae (i.e., krill), Foraminifera, Gymnosomata (i.e., sea angels, pteropods), Hydrozoa (i.e. jellyfish), Hyperiidea (i.e., amphipods), Myodocopina (i.e., ostracods), Mysidae (i.e., small pelagic shrimps resembling krill), Neocopepoda, Podonidae and Penilia avirostris (i.e., cladocerans), Sagittoidea (i.e., chaetognaths), Scyphozoa (i.e., jellyfish), Thaliacea (i.e., salps, doliolids and pyrosomes), Thecosomata (i.e., sea snails, pteropods), and four families of pelagic Polychaeta (i.e., worms) that are often found in the zooplankton and whose species are known to display holoplanktonic lifecycles (Tomopteridae, Alciopidae, Lopadorrhynchidae, Typhloscolecidae). The presence data associated with species belonging to these groups were retrieved from OBIS and GBIF between the 12/04/2018 and the 18/04/2018 using online queries via the R packages RPostgreSQL, robis and rgbif. Since the Neocopepoda infra-class comprise several thousands of benthic and parasitic taxa43, a preliminary selection of the non-parasitic planktonic species had to be carried out prior to the downloading using the species list of Razouls et al.43 as a reference. The spatial distributions of the groups cited above were first inspected using GBIF’s and OBIS’s online mapping tools to evaluate the potential number of overlapping observations between the two databases. As a result of their relatively low contributions to total observations/diversity, and very high overlap between databases, the occurrences of Cladocera and Polychaeta were retrieved from OBIS only (which usually harbors more occurrences). On top of the data collected from OBIS and GBIF, the copepod occurrences from Cornils et al.44 and the pteropod occurrences from the MAREDAT initiative45 were added to the dataset. This initial collection of zooplankton observations gathered 4,899,151 occurrences worldwide.
Then, similar criteria as Righetti et al.11,40 were applied to progressively remove those presences that would be discordant with estimates of contemporary open ocean zooplankton diversity. The number of observations and species discarded after each main step and for each initial dataset are reported in Supplementary Data 1. We discarded records that: (i) presented at least one missing spatial coordinate, (ii) were associated with an incomplete sampling date (d/m/y), (iii) were associated with a year of collection older than 1800, (iv) were not associated with any sampling depth, (v) were not identified down to the species level, and (vi) were issued from drilling holes or sediment core data. For step (vi), a list of keywords (Supplementary Note 3) was used to identify the names of those original datasets that contained either fossil or sedimentary records. These first steps resulted in the removal of 1,766,783 occurrences (~36%). Like for phytoplankton, the remaining occurrences were associated with surface salinity values from the World Ocean Atlas (WOA) 201346 and bathymetry levels from the National Oceanic and Atmospheric Administration (NOAA) using the marmap R package. Occurrences associated with salinity levels 500 m. The average depth was used when maximal depth was not provided in the metadata. Therefore, the maximal depth of a zooplankton species occurrence allowed in our dataset is 500 m. This way, we tried to account for the zooplankton community that frequently performs diel vertical migration across the euphotic zone or the mixed layer, and that often co-occurs with species inhabiting surface layers. This removed 109,582 (~6%) occurrences. Next, for each phylum, OBIS and GBIF datasets were merged and the list of species names were extracted. Every species name was then carefully examined and compared to the taxonomic reference list of the World Register of Marine Species (WoRMS; http://www.marinespecies.org) for all taxa but copepods, for which we used the taxonomic reference of Razouls et al.43. This way, we rigorously harmonized and corrected the species names across all datasets. In addition, we used the notes and attributes of WoRMS to identify whether species were holoplanktonic or meroplanktonic (i.e., those species that have at least one benthic phase in their life cycle). Jellyfish species usually display a fixed polyp phase during their life cycle, therefore we used the dataset of Gibbons et al.47 to remove the species that were not holoplanktonic. Overall, these steps discarded 37,234 occurrences (~2%). One of two duplicate occurrences were removed from the dataset if they displayed the same species name, sampling depth, sampling date, and if they occurred within the same 0.25° x 0.25° cell grid. This last step removed 900,446 occurrences (~54%), highlighting the high overlap between the two main data sources. The remaining 766,033 presences were binned into the monthly 1° x 1° grid cell of the WOA to match the spatial resolution of the environmental predictors. The average maximum (±std) sampling depth was 73 ± 109 m and the average sampling year was 1985 ± 21. Observation densities were spatially biased towards the North Atlantic Ocean and the Southern Ocean (Supplementary Figure 1). The data reflected the historical seasonal sampling bias towards spring and summer. In the northern hemisphere, observations were equally distributed from March to October but constituted 78% of the data. In the southern hemisphere, 75% of occurrences were sampled between November and March. The final dataset gathered occurrences for 2034 different species (576 genera, 161 families) spanning all the major zooplankton phyla and several size classes. The only notable missing taxa are those belonging to the Cercozoa and Radiozoa because they present little to no species-level observations in online biodiversity repositories as they have been historically overlooked by traditional sampling techniques48.
Contemporary environmental conditionsA comprehensive set of environmental variables that are known to affect the physiology and constrain the distribution of plankton was prepared to define the candidate predictors for the SDMs. The R packages mainly used were raster and ncdf4. First, twelve primary variables that are relevant for modeling the distribution of both phytoplankton and zooplankton taxa were identified49,50,51,52. These were then aggregated into twelve monthly climatologies at a 1° x 1° resolution (i.e., the spatial cell grid of the WOA). The first six primary variables were sea surface temperature (SST, °C), sea surface salinity (SSS), nitrates (NO3-), phosphates (PO43-) and silicic acid (Si(OH)4) surface concentrations (µM), as well as dissolved oxygen concentration (dO2, ml l−1). Oxygen limitations and oxygen minimum zones (OMZ) are key factors controlling the horizontal and vertical distribution of zooplankton53,54. However, the effects of oxygen are often confounded with those of temperature because surface oxygen scales linearly with SST on a global scale. Therefore, dO2 at 175 m depth was used instead of surface dO2. For all six variables, the twelve monthly climatologies of the WOA13v2 (https://www.nodc.noaa.gov/OC5/woa13/woa13data.html) were used. In addition, satellite observations stemming from the Sea-viewing Wide Field-of-view Sensor (SeaWiFS; https://oceancolor.gsfc.nasa.gov/) over the 1997 to 2007 time period were used to derive monthly climatologies of Photosynthetically Active Radiation (PAR; µmol m−2 s−1) and chlorophyll (Chl; mg m−3), the latter serving as a proxy for surface phytoplankton biomass. Monthly climatologies of mixed-layer depth (MLD, m) based on the temperature criterion of55 from the Argo floats data (http://mixedlayer.ucsd.edu/) were also considered. Climatologies of surface wind stress (m s−1) were obtained from the Cross-Calibrated Multi-Platform56, using data from 1987 to 2011 (https://podaac.jpl.nasa.gov/). Climatologies of surface carbon dioxide partial pressure (pCO2; atm) were obtained from the Surface Ocean CO2 Atlas (SOCATv2; https://www.socat.info/) and made available by Landschützer et al.57. Lastly, a variable depicting sub-mesoscale dynamics and the strength of sea currents was derived from the daily satellite altimetry observations over the 1993–2012 period (https://cds.climate.copernicus.eu/#!/home): mean Eddy Kinetic Energy (EKE, m2 s−2). EKE was computed from the northward and eastward components of surface geostrophic seawater velocity (assuming the sea level as geoid), following the method of Qiu & Chen58. Such variable enabled us to account for the potentially important role of sub-mesoscale activity in structuring plankton biodiversity59.Then, nine secondary predictors were derived from some of the predictors described above. PAR over the MLD (MLPAR, µmol m−2 s−1) was calculated following Brun et al.49 An estimate of annual range of SST (dSST) was added by computing the difference between the warmest and the coldest temperature across the 12 months. The excess of nitrate to phosphate (N*, µM) relative to the Redfield ratio was computed as [NO3-] −16[PO43-]. Changes in N* represent varying conditions of denitrification and remineralization from N2-fixing organisms7. The excess of silicates to nitrates (Si*=[Si(OH)4]-[NO3-], µM) was also computed to represent regions where silicates are in excess compared to what diatoms would need to use up the nitrates7. Si*  > 0 are indicative of conditions where diatoms can grow healthy. Since the distribution of macronutrients concentrations, chlorophyll concentration, and EKE values were all skewed towards lower values, we considered their logarithmic values (logNO3, logPO4, logSiOH4, logEKE, and logChl), based on either natural log or base 10, as additional predictors because they were much closer to a normal distribution.Species distribution modelingSDMs refer to a wide range of statistical algorithms that link an observed biological response variable (i.e., presence-only, presence/absence, abundance) to contextual environmental variables in the form of a response curve60. The latter is used to explore how a species’ environmental niche is realized in space and time24. In short, SDMs mainly rely on the following assumptions: (i) species distributions are not strongly limited by dispersal at a macroecological scale, an assumption valid for plankton considering the very strong connectivity of ocean basins through surface current on decadal scales61,62, which enables plankton species to display very large spatial ranges;11,47,60 (ii) species distributions are primarily shaped by the combinations of environmental factors that define the conditions allowing a species to develop. The latter assumption has been supported on macroecological scales, where the imprint of biological interactions (and dispersal) has been found to be relatively small63,64. Neither comparable abundance data nor presence/absence data were available from our datasets. In addition, presence-only data are less sensitive to discrepancies in species detection across various plankton sampling techniques. Therefore, based on species presences, we developed an exhaustive SDM framework to estimate plankton diversity patterns from an ensemble modeling approach65 that addresses the underlying main sources of uncertainties66,67.We follow the methodology developed in ref. 11, but simplify this approach to accommodate the limited predictor availability in the future model projections (see section “Choice of environmental predictors”), and the large number of diverse species we model in this work (sections “Background data (pseudo-absences)”, “Choice of environmental predictors”, and “SDMs evaluation and projections of monthly plankton species community composition”). We further derive SR based on habitat suitability rather than observed presence–absence data (section “Ensemble projection of global plankton species richness”). All methodological choices led to a minimization of computational cost and model complexity, while preserving all crucial patterns reported in ref. 11. Each methodological choice was carefully evaluated against other options, see sections below.
Background data (pseudo-absences)
Since we aimed at training correlative SDMs to model species distributions from presence data and environmental predictors, background data had to be simulated to indicate those conditions where a species is likely not to occur (i.e., pseudo-absences68). The generation of background data is a critical step in niche modeling experiments, and though no single optimal method has been identified by the niche modeling community, this step must address the important spatial and temporal sampling biases inherent to field-based observations. To do so, we made use of the target-group approach of Philipps et al.69, which has been shown to efficiently model phytoplankton distributions11. This method was found appropriate for our study because it generates background data according to the density distribution of the presence data, and therefore it: (i) does not induce additional bias to the initial biases in the presence data; and (ii) does not misclassify regions lacking observations (e.g., South Atlantic and Subtropical Pacific, Supplementary Fig. 1) as regions of absences.
For phytoplankton, we followed the background selection procedure described in Righetti et al.11. The authors used either the total pool of occurrences as a target-group, or defined three target groups based on the taxonomic groups contributing most to species diversity and observations (Bacillariophyta, Dinoflagellata, and Haptophyta). Background data of each species were randomly drawn based on the monthly resolved 1° x 1° occurrences of both their corresponding target groups, after applying an environmental stratification based on the SST and MLD gradients. This way, a species’ background is located at the sites where its lack of presence is most likely to reflect an actual absence. For each species, ten times more background data than presences were generated following the guidelines of Barbet-Massin et al.68. The amount of background data sampled from a specific SSTxMLD stratum was proportional to the number of monthly 1° x 1° cells provided by the target-group in this very stratum, thereby reflecting original sampling efforts. Both the total target-group background data (drawn from all sampling sites together) and the group-specific target background data were considered for our study, but both led to comparable estimates of phytoplankton species diversity (but see Fig. S3B of ref. 11).
The same method was applied for zooplankton species. First, we defined target groups based on their sampling distribution and broad taxonomic classification: Arthropoda (mainly copepods, but also krill and amphipods, that are sampled through similar techniques), Pteropoda, Chaetognatha, Cnidaria, Ctenophora, Chordata, Foraminifera, and Annelida. Unfortunately, the last three target groups displayed too few occurrences for drawing ten times more background data than presences. Consequently, their background data were drawn from the total pool of occurrences. Ctenophora also showed very few observations so they were merged with the Cnidaria as they are often considered together as jellyfish and collected in similar ways.
Total and target-group background data were drawn for all zooplankton species presenting more than 100 presences to run preliminary SDMs based on a preliminary set of predictors (Supplementary Fig. 2). These SDMs were then used to project preliminary diversity patterns for the four months that represent each season in the northern hemisphere (April, July, October, and January). These projections were examined for every group, and the predictive skills of the SDMs were evaluated using a repeated ten times split-sample test (see below). These tests showed that the total target-group background and the group target-group background converged towards SDMs of comparable skills and similar diversity patterns, except for the Chaetognatha, for which the target-group background leads to models of much poorer predictive skills (Supplementary Fig. 2). Furthermore, both background choices led to very similar latitudinal diversity gradients for phytoplankton (Fig. S3B in ref. 11). Therefore, to generate diversity patterns that are robust and consistent across the two trophic levels, the total target background data were used as standard background. The phytoplankton diversity pattern obtained with the total target background approach was only slightly lower in the Indo-Pacific and at very high latitudes11. Once they were generated, all background data were matched with monthly values for the 21 environmental variables described above.

Choice of the SDMs algorithms
The choice of the statistical method is a main source of uncertainty when projecting biodiversity scenarios through niche modeling66,67. Therefore, an ensemble forecasting strategy was adopted based on four types of SDMs that cover the range of algorithms types and model complexity that are commonly used:67,69 Generalized Linear Models (GLM), Generalized Additive Models (GAM), Random Forest (RF), and Artificial Neural Networks (ANN). The level of complexity of those models was constrained to avoid model overfitting70, a common pitfall when dealing with noisy and spatially biased data. SDMs including numerous predictors and parameterization features are more likely to fit spurious relationships and to be less transferable70,71. Consequently, the number of predictors was limited relative to the number of presences (see below) and the SDMs were tuned to fit relatively simple response curves. The GLM followed a binomial logit link, including linear and quadratic terms, and a stepwise bi-directional predictor selection procedure. The GAM also followed a binomial logit link. Smoothing terms with five dimensions, estimated by penalized regression splines without penalization to zero for single variables, were applied. Interaction levels between environmental predictors were set to zero for both GLM and GAM. The RF included 750 trees, and terminal node size was fixed at 10 to avoid having single occurrences as end members of some trees. The number of variables randomly sampled as candidates at each tree split (mtry parameter) was equal to the number of predictors used divided by three. The numbers of units in the hidden layers of the ANN, as well as the decay parameter, were optimized through five different cross-validations and a maximum of 200 iterations. Background data were weighted inverse-proportional to that of presence data (total weight = 1).

Choice of environmental predictors
To select for parsimonious and ecologically relevant sets of environmental predictors, a three-stage hierarchical selection framework was developed: (i) the distribution of the predictors’ values fitted to the presences were compared to their realized distribution between the main ocean basins to check whether one predictor could bias SDMs outputs towards a particular basin; (ii) pair-wise rank correlations between variables were examined, and one of two collinear variables was discarded where necessary; and (iii) models were trained to evaluate the explanatory power of several predictors sets of increasing parsimony, and rank the predictors within those sets at species-level. This selection procedure was carried out by separating phytoplankton from zooplankton since: (i) the two groups show different sampling distributions, and (ii) their niche dimensions might differ because of differences in their lifecycles (few days for phytoplankton, months to years for zooplankton) and biological requirements (photoautotrophy vs. heterotrophy and respiration). For those tests, only well-observed species with >100 occurrences were selected (nphytoplankton = 328; nzooplankton = 372). Ultimately, to account for the uncertainty in predictors choice, several final sets of predictors were defined based on the steps of the selection framework, and ensemble forecasting was adopted again (i.e., diversity estimates will be averaged across the sets of predictors).
Removal of variables impacted by sampling imbalances across ocean basins: Imbalance of sampling effort in geographical space can lead to sampling imbalance in environmental space if portions of an environmental gradient are strongly connected to an ocean basin that has been surveyed more extensively than others. To avoid such issues, the distributions of the annual values of the predictors were examined between the main basins (Arctic, Southern, Pacific, Indian and Atlantic Oceans). The most spatially imbalanced predictors were SSS and pCO2: the former is on average higher in the Atlantic Ocean, while the latter exhibits many of its most extreme values in the Peruvian upwelling system (Supplementary Note 3). The Peruvian upwelling is a hotspot of phytoplankton observations with clearly skewed observations (and the number of species sampled) towards pCO2 values  > 400 atm (Supplementary Note 3). Plus, the pCO2 data do not cover the Arctic Ocean, the Mediterranean Sea, and the Red Sea. Consequently, pCO2 was discarded from the list of predictors to avoid strong sampling bias effects on SDM projections. A majority of the zooplankton data are concentrated in the Atlantic Ocean (Supplementary Fig. 1). As a result, the distribution of SSS values fitted to zooplankton occurrences is skewed towards SSS values >35 (Supplementary Note 3). As SSS is commonly used as a predictor for modeling the distribution of zooplankton47,48,49, we wanted to further examine its potential to act as a basin indicator rather than a predictor meant to represent an actual environmental control on species distribution. To do so, we performed ensemble SDMs projections for the zooplankton species, based on three variables sets: (i) without SSS, (ii) with SSS, and (iii) with Longitude (0°–360°) instead of SSS. Variables sets (ii) and (iii) led to very similar global zooplankton SR patterns with hotspots in the Atlantic Ocean. On the contrary, (i) led to more balanced zooplankton SR between basins without significantly lowering SDMs skills (Supplementary Note 3). We interpreted this as a bias in environmental space towards the conditions prevailing in the Atlantic Ocean, therefore we chose to discard SSS from the list of predictors.
Removal of collinear variables: Strong correlations among predictors can mislead the ranking of variable importance in SDMs72, so it has become common practice to exclude one of two variables that are highly collinear. Pair-wise Spearman’s rank correlation coefficients (⍴) were computed based on the predictors’ values fitted to the presences. When two variables exhibited a |⍴|  > 0.70, the one displaying the distribution closest to a normal distribution was kept. From phytoplankton occurrences, we identified two clusters of strongly correlated variables: one comprising MLD, PAR, MLPAR (by construction), and wind stress (but PAR and MLD were only correlated at ⍴ = −0.66); and the other one comprising [NO3-], [PO43-] and their logged values. Similar clusters were found from zooplankton data, except that PAR was slightly less correlated to Wind stress (⍴ = −0.66) and MLD (⍴ = −0.58), and that [NO3-], [PO43-], plus their logged versions, showed stronger correlations with SST (⍴ = −0.80). As [NO3-] is a key factor for structuring planktonic systems7, and because we aimed to keep the variables sets as consistent as possible between species, logNO3 was kept as a candidate predictor. The variables retained for phytoplankton were: SST, dSST, logEKE, Si*, N*, logSiOH4, logNO3, logChl, wind stress, PAR, MLPAR, and MLD. The last four variables were kept to explore the outcome from alternative choices in the variables sets (but see below). The variables retained for zooplankton were the same but with the addition of dO2.
Examination of the explanatory power of predictors sets and ranking of predictors: To further evaluate which subset of these variable subsets are key to model species distributions, GLM and RF were performed for each species for several sets of decreasing complexity (from ten to five predictors), and the adjusted R2 of the models, as well as the ranking of predictors within each set, was extracted (Supplementary Note 1). GLM and RF were used here because they are part of the SDMs that will be used for projections afterwards and because they represent maximally different inherent model complexities among the SDM types used73. For GLM, predictor importance was determined according to their absolute t-statistic using the caret R package. For RF, predictor importance was based on the Gini index, which measures the mean decrease in node impurity by summing over the number of splits (across all trees) that includes a variable, proportionally to the number of samples it splits. The ranger R package was used for assessing variable importance with RF models. To keep the variables ranks comparable across predictors sets, rank values were normalized to their maximum. For each model type, the distributions of the models’ R2 and the distribution of the predictors’ ranks were examined for phytoplankton and zooplankton separately. The same was done between the main groups constituting the phytoplankton (Bacillariophyta, Dinoflagellata and Haptophyta) and the zooplankton (Copepoda, Chaetognatha, Pteropoda, Malacostraca, Jellyfish, Chordata and Foraminifera). This allowed us to identify the most important predictors for modeling the species distributions and to evaluate if a decrease in the models’ skill was linked to the removal of certain variables. Group patterns allowed us to test whether different groups differed in their main environmental drivers.
For phytoplankton species, 14 sets of variables were examined (Supplementary Note 1). The first nine aimed to test: (i) the impact of alternative choices between variables that were identified as collinear (wind vs. MLPAR vs. MLD + PAR); (ii) the impact of progressively discarding variables that initially presented lower ranks (logEKE, Si*), and (iii) the impact of choosing logNO3 over logSiOH4, two variables representing global macronutrients availability and that present relatively high correlation coefficient (⍴ = 0.59). The last five sets of predictors (10–14) aimed to test the impact of alternatively removing those variables that presented relatively high ranks in the previous sets: SST, dSST, N*, logSiOH4, logChl, PAR. In a similar fashion, 15 sets of variables were tested for zooplankton (Supplementary Note 1). The first ten aimed to test: (i) the impact of choosing wind stress over MLPAR or over MLD + PAR; (ii) the impact of selecting PAR over MLD; (iii) the impact of discarding Si*, N*, logEKE; and (iv) the impact of choosing logNO3 over logSiOH4 (⍴ = 0.64). The last five sets of predictors (11–15) aimed to test the impact of alternatively discarding the top five predictors: SST, dSST, dO2, logSiOH4, and logChl.
GLM and RF converged towards similar median variable rankings and evidenced high inter-species variability (Supplementary Note 1). For total phytoplankton, GLM identified the following median ranking across all species: SST  > N*  >logChl  > logSiOH4 and dSST  > logNO3  > logEKE  > Si* > the PAR/MLD/MLPAR/wind stress cluster. RF ranked predictors in the following median order: SST and N*  >logChl  > dSST  > logSiOH4  > logEKE and logNO3  > Si* > the PAR/MLD/MLPAR/wind stress cluster. Yet, both GLM and RF also identified PAR as a major predictor for Haptophyta, which does not appear in the rankings for total phytoplankton because Haptophyta represented ~9% of species composition only. Since adding PAR does not alter the models’ R2 for the Bacillariophyta and Dinoflagellata, it was retained for the final predictors sets. For total zooplankton, GLM ranked predictors in the following median order across all species: SST  > dSST and logSiOH4  > logChl and logEKE  > dO2 and logNO3  > N*  >Si* and MLPAR  > wind stress  > PAR and MLD. RF identified the following median ranks: SST  > dSST and dO2  > logNO3  > logSiOH4  > logChl and logEKE  > N* and Si*  > PAR  > wind stress  > MLD and MLPAR. These median rankings reflected those of the Copepoda since they represented >70% of all zooplankton species. Again, rankings displayed high variance, reflecting high inter-species variability. Overall, based on all the results shown above, eight different final predictors sets were kept for modeling the distribution of phytoplankton (n = 4) and zooplankton (n = 4) consistently. In contrast to ref. 11, predictor ensembles were defined across all species rather than for each species. This was due to multiple reasons: (i) predictor availability for future model projections was limited and did not allow for species-specific variable choices, (ii) computational constraints with regard to the total number of ensemble members that could be projected, (iii) the five sets already contain those predictors that explain a majority of the variability in most models, (iv) recent findings from Righetti et al. (in prep.) that the uncertainty due to predictor choice is low for models with optimized background selection.
Phytoplankton:

1.

SST, dSST, logChl, N*, PAR, and logNO3

2.

SST, dSST, logChl, N*, PAR, and logSiOH4

3.

SST, dSST, logChl, N*, PAR, logNO3 and Si*

4.

SST, dSST, logChl, PAR, and logNO3

Zooplankton:

1.

SST, dSST, dO2, logChl, and logNO3

2.

SST, dSST, dO2, logChl, and logSiOH4

3.

SST, dSST, dO2, logChl, logSiOH4, and N*

4.

SST, dSST, dO2, logChl, logNO3, and Si*

SDMs evaluation and projections of monthly plankton species community composition
Only species with more than 75 presences were considered for modeling plankton species distributions (nphytoplankton = 348; nzooplankton = 541) because we aimed to achieve a relatively high presence-to-predictors ratio (~15, which is the ratio achieved for a species with 75 presences and five to six predictors) to be more conservative than Righetti et al.11 (i.e., minimum 24 presences) since we aimed to project the SDMs in future conditions based on a pool of species for which we have high confidence. This is in line with Guisan et al.60 who suggest to maintain at least a ratio of ten. For each species, each SDM, and each set of predictors, presences and background data were randomly split into a training set (80%) and a testing set (20%) and these evaluation tests were repeated ten times. Therefore, 160 (four SDM types x four predictor sets x ten separate evaluation runs) models were trained per species, resulting in a total of 142,240 SDMs. Model skill was evaluated based on two widely used metrics: the True Skills Statistic (TSS74) and the Area Under the Curve (AUC60). TSS values range between −1 and 1, with null values indicating that models perform no better than at random. AUC ranges between 0 and 1, with values 0.30 were retained for the final ensemble projections (Supplementary Fig. 3).
In total, 860 species were considered as successfully modeled (nphytoplankton = 336; nzooplankton = 524; Supplementary Data 2). For those, each of the 160 SDMs was projected onto the twelve monthly climatologies of its corresponding predictors set and the projections were averaged over the ten cross-evaluation runs. This way, we obtained global maps of monthly mean presence probability for each of the 16 SDM x predictor set combinations. These maps are to be interpreted as habitat suitability patterns that highlight the regions of the global ocean where the environmental conditions are most favorable for a species to develop. Habitat suitability maps were not converted to binary presence–absence maps as probabilistic outputs provide more gradual responses that should better reflect the very dynamic occupancy patterns of plankton and that are better suited than threshold approaches for our purposes75,76,77. For each grid cell of the global ocean, the probabilistic estimates of species habitat suitability were stacked to obtain monthly estimates of species composition. All SDMs were trained, evaluated, and projected using the biomod2 R package.

Ensemble projection of global plankton species richness
For every SDM x predictor set combination and every month, we summed the species habitat suitability to estimate monthly SR. Then, annual average SR was estimated for each cell. The annual average was preferred over the annual integral because of the high latitudes that presented a lot of missing values in winter because of the lower coverage of satellite products. Since this diversity estimate is the sum of habitat suitability indices, it is to be interpreted as the amount of SR that the monthly/annual average environmental conditions should be able to sustain (i.e., potential SR). This way we obtained 16 estimates of annual SR for total plankton (all 860 species together), phytoplankton and zooplankton. Ensemble projections of annual SR were then obtained for these three categories by averaging the annual SR estimates.
As the biological data used to train the SDMs span several decades (mostly between the 1970s and 2000s), our diversity estimates are integrative of changes in SR and species composition (i.e., changes in beta diversity) that occurred during these decades. The phytoplankton species modeled are mainly members of the Bacillariophyceae (45.8%), and the Dinoflagellata (45.8%), which usually rank among the large marine microalgae. Therefore, the phytoplankton SR estimates shown here should be mainly representative of the microphytoplankton (20–200 µm) rather than smaller size fractions. Nearly half (51.9%) of the zooplankton species modeled are Copepoda, making it the most represented groups in the zooplankton SR patterns followed by: Malacostraca (crustacean macrozooplankton such as Euphausiids and Amphipods; 13.9%), Jellyfish (13.1%), Foraminifera (5%), Chaetognatha (5%), Pteropods (4%), Chordata (3%). The last 4% of species modeled are a mix of Annelids and Branchiopods. The full list of the species modeled as well as their taxonomic classification is given in the Supplementary Data 2.
We underline that the 860 species modeled are those for which we have enough records for training reliable SDMs. Therefore, these species are likely to be those that are the most frequently detected by conventional sampling and identification techniques, either because: (i) they are the ones dominating total plankton abundance, which makes their collection more likely; or (ii) they are larger species (in terms of cell volume or body length), which would facilitate their sampling and identification under the binocular or recent imaging systems. We acknowledge that our approach does not allow us to account for rare taxa and thus under samples the true diversity of the marine plankton. Nonetheless, we argue that our approach does allow us to estimate global plankton diversity patterns as the species dominating plankton abundances are those carrying biogeographical information30, meaning their distribution and abundance patterns can be correlated to environmental gradients. Meanwhile, the patterns of rare and non-dominant species, which constitute the majority local SR, exhibit no biogeographical signature30. This has been supported in the previous study of Righetti et al.11 where the authors showed that global SR patterns were robust to the progressive exclusion of taxa with relatively few records.
We also acknowledge our estimates of species distribution might be biased by imbalances in species detection and sampling effort between sampling cruises, as those rely on a wide range of collection and identification methodologies. We argue that such biases are particularly significant when relying on abundance data, and that we mitigate them by: (i) converting all observations to presences and aggregating them onto a 1° x 1° grid; (ii) modeling SR as an emergent property overlapping the distribution of single species with equal weighting rather than modeling SR directly, in which case diversity estimates would be highly sensitive to sampling effort imbalances; (iii) by designing the SDMs in a way that accounts for spatial and temporal sampling biases in geographical and environmental space; and (iv) by tuning down the complexity of the SDMs (i.e., reduced number of features and predictors) in order to avoid model overfitting70.
Future environmental conditions in the global surface oceanThe future monthly fields of the selected environmental predictors were obtained from the projections for the 2012–2100 period of five ESM simulations for the IPCC’s RCP8.5 scenario from the MARine Ecosystem Model Intercomparison Project (MAREMIP, http://pft.ees.hokudai.ac.jp/maremip/index.shtml78) and/or the Coupled Model Intercomparison Project 5 (CMIP579). The model ensemble contained the following five ESMs (with their embedded ocean and ecosystem models indicated indicated in brackets after the semicolons): Community Earth System Model version 1 (CESM1, POP-BEC), Geophysical Fluid Dynamics Laboratory Earth System Model with Modular Ocean Model version 4 (GFDL-ESM2M; MOM-TOPAZ), Institut Pierre Simon Laplace Climate Model version 5A-LR (IPSL-CM5A-LR; NEMO-PISCES), Centre National de Recherches Météorologiques Climate Model version 5 (CNRM-CM5; NEMO-PISCES) and the Model for Interdisciplinary Research on Climate version 5 (MIROC5; MRI.COM-MEM). All ESMs were fully-coupled except for MIROC5 for which the ocean model was forced by the atmospheric component. All of the projections used here were benchmarked, quality-controlled and described in the previous multi-model comparison studies of Laufkötter et al.26,80. Considering the scope of the present study, we refer to these authors’ previous extensive descriptions for the full detail of the ESMs used here. Taken together, the present five ESM ensemble gathers models of various sensitivity to future climate forcing, and thus provides a wide range of alternative environmental conditions projected for the future surface ocean. With the present ESM ensemble, we account for the variability in the choice of the climate model, which is known to be a significant source of uncertainty in biodiversity projections; this source being consistently lower than those associated with SDM choice, though28,66,67.The monthly projections of the five selected ESMs were interpolated on the 1° x 1° cell grid of the WOA (i.e., the one used to train our SDMs) over the 2012–2100 period for all the nine chosen environmental predictors. To obtain future monthly climatologies that span a comparable amount of temporal variability as the in situ climatologies used to train the SDMs (~20 years), a baseline and an end-of-century time periods were first defined (2012–2031 and 2081–2100, respectively) for every ESM projection run. The 12 monthly climatologies were derived based on the models’ monthly projections and monthly anomalies were computed by subtracting the baseline values to the end-of-century ones. For dSST (i.e. annual range of SST), the annual maximum of SST was derived from the monthly climatologies and the difference between the baseline and the end-of-century dSST provided the delta value. These anomalies can be either positive or negative and they represent the difference in the predictors’ condition due to future climate change under the RCP8.5 GHG concentration scenario25. To obtain the final conditions prevailing in the surface ocean for the end-of-century period, the delta values were simply added to the in situ climatologies representing the conditions in the contemporary ocean. The SDMs of the 860 plankton species successfully modeled were then projected onto these future monthly climatologies for each of the ESM. This way, we estimate the monthly probability-based species composition in the future global ocean for each of the 80 combinations of SDMs (n = 4), ESMs (n = 5), and predictor set (n = 4). Overall, our ensemble forecast approach65 generates an unprecedented set of 825,600 species-level estimates of global future habitat suitability patterns. Finally, mean annual SR and community composition were calculated for total plankton, phytoplankton and zooplankton for each of the 80 possible combinations of projections, as described in section “Ensemble projection of global plankton species richness”.Analyses
Ensemble projections of changes in species richness, community composition turnover, and changes in species associations between the contemporary and the future ocean
For each of the 80 projection combinations described above, the mean annual SR estimates for the contemporary ocean were subtracted to their corresponding mean annual SR estimates for the future ocean to compute the percentage difference in mean annual SR (%∆SR) for total plankton, phyto- and zooplankton. The %∆SR represents the emergent change in SR caused by future climate change(s) through changes in species-level habitat suitability patterns. While changes in SR indicate climate change impacts on plankton alpha diversity, these do not inform us on the potential impacts on beta diversity (i.e., changes in community composition81). A community that experiences the replacement of all its constituting species by an equivalent number of newcomers will display a 100% rate of community turnover but no changes in SR. To investigate the amplitude of global plankton species turnover triggered by climate change, we examined future total turnover in annual species composition using Jaccard’s dissimilarity index while decomposing its two additive components: nestedness (i.e., changes in SR) and true turnover (ST), which indicates the % of species that will be replaced in a community27 using the betapart R package.
To do so, the mean annual species habitat suitability patterns used to estimate the ensemble changes in SR had to be converted to presence–absence maps as the Jaccard’s dissimilarity index requires binary inputs80. A range of thresholds (0.10 to 0.80, by steps of 0.01) was first explored for each SDM type (GLM, GAM, ANN, and RF) to infer threshold-based annual SR patterns. Then, we quantified the similarity of the threshold-based annual SR vs. the probability-based annual SR using Spearman’s rank correlation coefficient (⍴) and ordinary linear regressions (R2) to identify the range of thresholds that best match the probability-based estimates. The 0.25–0.40 range provided the most similar global SR patterns for GLMs, GAMs and ANNs (all ⍴  > 0.95, and all R2  > 0.90). The 0.10–0.25 range was chosen for RF models. These ranges largely overlap with the species mean probability thresholds that maximize the TSS/AUC evaluation metrics, which are commonly used to convert habitat suitability into presence–absence maps. However, the maximizing-TSS approach tends to underestimate the natural gradual response of organisms to environmental variations, which is particularly problematic when dealing with SR patterns of widely-dispersed and climate-sensitive ectotherms such as the plankton75,76,77. Therefore, we chose to rely on a range of thresholds instead as it enables us to account for a wider range of possible realizations of community composition and better reflect the dynamic occupancy patterns inherent to planktonic taxa. Consequently, we derived ST estimates in annual species composition for each of the SDM-dependent threshold mentioned above and every of the 80 annual projections combinations, for total plankton, phyto- and zooplankton separately. Again, the ensemble projection in annual ST was derived by averaging those projections.
We further examined how climate change could impact not only community composition but also those species associations within the community that represent potential biotic interactions, which support ecosystem functioning3,20,21. Based on the mean annual species composition estimates used to compute ST rates, a text analysis algorithm82,83 was used to identify pairs of species, which co-occur more frequently than expected given their individual occurrence. In short, the text analysis algorithm assigns an association score to each possible pair of two plankton species in all grid cells based on a likelihood ratio (LLR)82. The latter compares the probability of two species co-occurring together to the probability of one species occurring without their partner (i.e., two alternative probabilities), or when both are projected as absent, based on a combination of Shannon’s entropy indices (H’). LLR values are >0 and they scale with the significance level of the projected species association, whether it is a positive (co-occurrence) or a negative (one-sided occurrence or co-absence) association. To disentangle between those two cases, when a species pair displayed an observed co-occurrence frequency lower than the product of the one-sided occurrence frequencies normalized to sample size, its LLR value was multiplied by −1. This way, we can identify those species pairs whose co-occurrence probability is lower than the products of the two one-sided occurrence probabilities (LLR  More

Geochemistry and ecotoxicology of AMDAMD systems are an important source of metal/metalloid pollution to the receiving hydrosphere with devastating consequences on the biological drivers of affected ecosystems. Environmental menaces of AMD have not been exhaustively reported worldwide. Scanty information exists across Africa and many developing economies. The homogenised mixture of detached biofilm and AMD samples from a derelict coal mine at three sampling periods were assayed for geochemical delineation and analysed for pollution intensity against reference background geochemical values. The measured values of the physical properties and contents of selected HMs in drains from a coal mine in Nigeria were as presented in Supplementary Table A.1. Virtually all the measured parameters exceeded the permissible limits of WHO guidelines for potable water. The AMD water was acidic (pH = 3.1 ± 0.265), and contained characteristic anions that are common to AMD including dissolved sulphides (1.37 ± 0.233 mg l−1), sulphates (313.0 ± 15.9 mg l−1), carbonate (253.0 ± 22.4 mg l−1) and nitrate (86.6 ± 41.0 mg l−1) above the allowable limits of WHO. Although the acidic pH of AMD in the present study compares well with those associated with mines in Russia14, more extreme acidic pH values have been reported in other climes. Negative pH values of − 1.56 and − 3.6 were observed in AMD from Iberian Pyrite Belt20 and Richmond Mine at Iron Mountain, USA21, respectively. The values of physicochemical parameters associated with the AMD from Onyeama were similar to data reported for other mine wastewaters in Nigeria22 and elsewhere4. It is known that sulphide minerals, in presence of water and oxygen, oxidise to sulphate as observed in the elevated sulphate concentration (313 ± 15.9 mg l−1) in the present study. The low pH observed in the AMD is due to the formation of sulphuric acid from sulphate in presence of protons (H+). This consequently causes the leaching of metal/metalloid ions into the drains. The concentrations of dissolved organic matter in AMD tends to be relatively low ( Co  > Pb  > As  > Ni  > Cr  > Fe (Table 1). Enrichment of five HMs was exceptionally high (Cd  > Co  > Pb  > As  > Ni), while Cr and Fe were very high and moderately enriched the AMD water, respectively. The astronomically high contamination and enrichment factors of the AMD signified the enrichment potentials the AMD portends on receiving surface waters. The AMD from the Onyeama coal mine has been reportedly impacting the water qualities of rivers within the location25. It is assumed that the extremely high concentrations of toxic metals/metalloids in the AMD dilutes out upon discharges into nearby rivers, contaminating the surface water and raising the bioavailable metals/metalloids beyond safe thresholds. Further reports of toxic metals/metalloids enrichment of surface waters via inflow of AMDs from other mines in Nigeria26 and other climes3,4,27 are worrisome and oblige mitigations.Table 1 Physico-chemistry, pollution and ecological impact determinants of heavy metals and metalloid contained in the AMD from coal mine.Full size tableThe HMs-enriched environments inadvertently exert ecotoxicity unto the drivers of the ecosystems. The level of HMs accumulation to the organic matter in the AMD, through geo-accumulation (Igeo) index of Fe (7.60 ± 0.779) to Cd (20.9 ± 0.075) (Supplementary Table A.2), was very severe and in a similar order to CF. It possibly implies organic matter in the AMD harbours the mobile toxic metal/metalloid concentrations and make them available to the food web28. Thus, biomagnification of the toxic metals/metalloids along the trophic level becomes palpable and a challenge to the biota of any surface water receiving the AMD and to public health21,28. Ecological risk assessments define and categorise the pollution status of ecosystems with the HMs contained in the AMD. Based on the potential ecological risk factor (Er), Cd exerted an extremely high-risk index (36.3 ± 1.96 × 106), and none of the metals/metalloids exercised less than 1000 risk index (Supplementary Table A.2). All the HMs/metalloid contained in the AMD posed very high ecological risks and could be categorised in the order of Cd  > Co  > Pb  > As  > Ni  > Cr  > Fe. The modified potential ecological risk factor (MEr), however, stipulated that five HMs posed a very high risk in the order: Cd  > Co  > Pb  > As  > Ni, whereas Cr and Fe were determined to be of considerate and low risks, respectively. The HMs exerted high risk to the AMD ecosystem as calculated by ecological risk quotient (RQ) in the order: Pb  > Cd  > As  > Ni  > Co  > Fe  > Cr. The ecological risk index of all the HMs as a whole was very high (375,000 ± 22,400) index as stipulated by the modified potential ecological risk index (Table 1). The prodigiously high ecological risks indexes of the HMs/metalloid in the AMD indicated grave danger the AMD would portend on the surface- and ground-waters.Microbial community structure of AMD from Onyeama coal mineA total of 26,160 and 40,403 valid (filtered) sequence reads were obtained for bacteria and eukarya, respectively, after a quality check of biofilm-water amplicon sequence data. The valid sequences were clustered into 2036 and 1002 operational taxonomic units (OTUs) of bacteria and eukarya domains of life, respectively, as presented in Table 2. Microbial community structures are sensitive descriptors of ecological stressors pivotal to understanding ecosystem functions29. The number of clustered high quality, non-chimeric sequences as OTUs based on UCLUST and CD-HIT against the sequence reads was depicted as asymptotic rarefaction curves (Supplementary Fig. A.1). The curves revealed that higher numbers of OTUs were delineated from valid sequence reads of 16S rRNA genes, unlike the lesser number of OTUs obtained from valid sequence reads of ITS2 region located between 5.8S and 28S rRNA genes of eukaryotes. The OTU richness observed in the rarefaction curves established coverage of the majority of species and was further validated with the richness and diversity estimations presented in Table 2. Despite the higher number of valid sequence reads obtained from the amplified ITS2 (40,403) than that of 16S rRNA genes (26,160), the observed OTUs were more in 16S rRNA genes (2036) than those of ITS2 (1002). More than 99.8% and about 98.5% of the microbial community in AMD from the Onyeama coal mine represented eukarya and bacteria OTUs, respectively, based on estimated Good’s library coverage. The coverage degree of the MiSeq sequencing corroborated the rarefaction curves. Furthermore, the estimated OTU richness (based on higher values obtained from ACE, Chao1 and JackKnife indexes) showed that bacterial phylotypes were richer than those of eukarya. Alpha diversity indexes (NPShannon, Shannon, and inverse Simpson) phylogenetic diversity index revealed that bacteria in the AMD were more diverse than eukarya OTUs.Table 2 Alpha diversity of microbiome evenness, richness and varieties of species in the sediments.Full size tableTaxonomy and phylogeny of microbial OTUs in AMD from coal mineThe taxonomic composition and relative abundances of the AMD microbiome, as shown in Fig. 1, revealed that the bacterial community spanned 10 phyla whose sequence reads were at least 1% (Fig. 1a). Whereas the eukarya domain of life (with sequence reads ≥ 1%) found in the AMD include Fungi, Plantae and Animalia kingdoms (Fig. 1b). Ascomycota, unclassified Fungi phylum (Fungi_p), Basidiomycota, and Mucoromycota represented Fungi kingdom, while Ciliophora and Arthropoda phyla were Animalia and Chlorophyta phylum epitomised Plantae kingdom. Association of the domain Eukarya (comprising Alveolates, Chlorophyta and Fungi as observed in this study) with AMD is reported to a lesser extent when compared with Bacteria30. The Fungi, largely represented by Ascomycota and Basidiomycota, are primarily found in sub-surface low-pH biofilms thriving in AMD31. While the Alveolates are suggested to have acted as primary/secondary consumers, the amoebae were secondary grazers in the AMD ecosystem29,32. Fungi taxa must have participated in carbon cycling as the main decomposers in the microbial community of the AMD. The taxonomic composition and relative abundance of phyla regarded as ‘Others’ (sequence reads  50%). Evolutionary analyses were conducted in MEGA6.Full size imageUrease-producing bacteria instigate insoluble metal-carbonate micro-precipitation through urease activity16. The growth-time courses and urease activities of the bacteria consortium in simulated AMD were presented as curves (Fig. 5). It was observed that the impact of high concentrations of HMs cocktails was not pronounced beyond the early 6 h post-inoculation, which was regarded as the lag phase. The bacteria consortium might have activated necessary genes needed to tolerate and sequester the metals/metalloids toxicity during the lag phase without cell multiplications. Afterwards, the bacteria consortium grew steadily with the production of urease, based on increasing measurement of urease activity, as incubation continued. At 30 h post-inoculation, 245.3 (± 23.7) U ml−1 activity of urease was observed in broth without a toxic metal cocktail. However, more urease activity (255 ± 7.6 U ml−1) by the bacteria consortium was observed in medium amended with low concentrations of metal cocktails unlike lesser activities of 235 (± 7.6) U ml−1 and 193.7 (± 10.7) U ml−1 associated with medium and high metal concentrations, respectively. As the growth remains stationary and pH further increased to  > 8.2, urease activities were at least 253 U ml−1 in all the cultures. Although urease activities at acidic pH have been reported in acid-tolerant human pathogens19, the findings in this report were assumedly the first amongst bacterial strains from AMD-polluted environments. The urease activities at acidic pH compared favourably with activities at alkaline pH in previous studies7,16,42,44. Moreover, the pH of the culture system kept increasing, alleviating the acidity condition that initially prevailed in the AMD system.Figure 5Growth kinetics of bacterial consortium via viable counts extrapolated into optical density at 600 nm wavelength (a) and growth-dependent urease activity of bacterial consortium (b) in TGYM broth without heavy metals (HMs) cocktail, and with low, medium, and high concentrations of HMs cocktails. Low HMs concentrations cocktail comprised (per liter) Cd, 27.9 mg; Pb, 118.7 mg; Co, 16.2 mg; Ni, 16.2 mg; and As, 61.5 mg. While medium HMs concentration contained (per liter) Cd, 55.7 mg; Pb, 237.3 mg; Co, 32.4 mg; Ni, 32.3 mg; and As, 123.1 mg. High HMs concentration contained (per liter) Cd, 139.3 mg; Pb, 593.3 mg; Co, 81.1 mg; Ni, 80.7 mg; and As, 307.6 mg. The mean pH at the beginning of experiment was 3.5 and rose to 8.2–8.4 at 48 h post-inoculation. Growth kinetics at exponential growth phase are in the inserts of panel (a), where ‘Td’ represents doubling time and ‘K’ is the growth rate at exponential growth phase. Error bars represent standard error mean (SEM) of triplicate experiments. The culture conditions were as explained in the “Methods” Section (Growth kinetics and urease activity of bacteria consortium; Determination of bacterial growth-dependent HMs/metalloid sequestration in simulated and natural AMD).Full size imageInterestingly, urease activity was observed in low quantity at acidic pH, unlike higher activity when the pH inclined towards alkaline (Fig. 5). It is proposed that urea finds its way into Onyeama coal mine drains through runoff from agricultural soils fortified with urea fertilizers and animal manures, which are common agricultural practices in Nigeria. The products of urea hydrolysis might have equilibrated in water to form bicarbonate, ammonium and hydroxyl ions that serially increased the culture pH. Ultimately, the bicarbonate equilibrium might have shifted to form carbonate ions (HCO3− + H+ + 2NH4+ + 2OH− ↔ CO32− + NH4+  + 2H2O) that enhanced the metal-carbonate micro-precipitation (Me2+  + Cell → Cell-Me2+ + CO32− → Cell-MeCO3). The gradual increase in pH could have further indulged the formation of CO32− from HCO3−, leading to metal-CO3 precipitation around cells and in culture media. Bicarbonates enrichment with inherent ammonia production was thought to have provided additional acid neutralization of the AMD. The growth kinetics after the presumed lag phase in the early 6 h to late exponential phase at 18 h showed that a low concentration of HMs cocktails did not have an impact on the growth of the bacteria consortium. Consequently, the bacteria consortium exhibited excellent sequestration of multi-component toxic HMs in both the simulated toxic metal-rich AMD and the actual AMD obtained from the Onyeama coal mine (Table 3).Table 3 Growth associated sequestration and precipitation of heavy metals/metalloid cocktail and AMD from Onyeama coal mine.Full size tableThe bacteria consortium displayed more than 94% efficiency of Cd and Pb sequestration in natural AMD, while 100% efficiency was observed in all the simulated AMD treatments (Table 3). Low performance was found with Ni and As, but not less than 70% sequestration efficiency was observed in all treatments. Efficient sequestrations of toxic metals, up to 100% removal efficiency of most toxic metals, observed with the bacteria consortium were similar to findings in a previous study13. Mixed-bacterial cultures are known to be able to perform more complex tasks and survive in more unstable environments than a monoculture. Nevertheless, 89.3–98% removal efficiencies of Ni, Pb, Co, and Cd from solution have been reportedly achievable with urease-producing Sporosarcina koreensis45. Similarly, Bacillus sp. KK1 reportedly mitigated lead-contaminated mines tailings containing mobile Pb (1050 mg kg−1) to form insoluble precipitates of PbS and PbSiO334. Growth-dependent sequestration of HMs cocktails by the bacteria consortium was adduced to be via precipitation. The weight of the precipitates was evaluated to be proportional to concentrations of HMs cocktail present. The bacteria consortium was observed to drive the formation of as much as 15.6 (± 0.92) mg ml−1 precipitates (Table 3) that were assumed to be in form of HMs-carbonates in TGYM supplemented with high concentrations of HMs cocktail within 24 h post-inoculation. In natural AMD bio-stimulated with urea and seeded with bacteria consortium for 24 h, 10.5 (± 0.52) mg ml−1 HMs precipitates was observed unlike 8.57 (± 2.52) mg ml−1 precipitates obtained from natural AMD toxic metals sequestration without urea fortification. It appeared that the quantity of toxic metal precipitate was proportional to quantities of available toxic metals, which corresponded to the number of heterogeneous nucleation sites on the surface of the bacterial cells. Omoregie et al.42 reported a relatively similar quantum of precipitation as CaCO3 with species of ureolytic Firmicutes isolated from limestone caves. As such, there was no correlation between urease activity and quantum of toxic metal precipitation since there is a likelihood that other metabolic activities may be linked to urease activities. Nevertheless, the bioremediation strategies demonstrated in the present study exhibited excellent toxic metal sequestrations unlike insignificant (p  > 0.05) natural attenuation process of the autochthonous community without augmentation with bacteria consortium and stimulation with nutrients (as presented in Table 3).In conclusion, AMD from the Onyeama coal mine is a point source of pollution to the surrounding environments because of its richness in anions and toxic metals/metalloids. It has a high potential of enriching the receiving hydrosphere with toxic metals/metalloids and exerts severe ecological risks (Er  > 320) with Cd and Pb wielding a huge critical risk index (38.1 ± 2.18 × 106) on the biological elements of the ecosystems. The dominance of Proteobacteria (50.8%), Bacteroidetes (18.9%), Ascomycota (60.8%), and Ciliophora (12.6%) characterised the microbial community of the AMD, where unclassified OTUs occurred mostly among the species. Enrichment of the AMDs skewed the bacterial community as depicted in the alpha diversity indexes against that of coal AMD leading to the selection of bacteria consortium with an excellent potential of stemming the toxicants in the AMD. The bacteria consortium efficiently removed toxic metals/metalloids ( > 70%) through precipitation and simultaneously neutralised AMD acidity. The bacteria consortium exhibited appreciable urease activity ( > 190 U ml−1), through which the precipitation was assumed possible via the formation of metal/metalloid-carbonates. The bacteria consortium is suggested to be a sustainable biotechnological candidate in designing a bioremediation strategy for decommissioning AMD before discharge into the surrounding environment. More

Convergence is a common feature of evolution and has great effect on the succession of microbial communities. For natural microbial communities such as the microbiome of gut [1], soil [2], sediment [3], rhizosphere [4], and phyllosphere [5], convergence generally means that different communities converge towards a similar species composition, which is accompanied by species loss and acquisition. Such a convergence can be reproduced in simplified synthetic communities [6,7,8], or even in single-species populations, in which convergence can still be achieved at sub-species level [9, 10]. Unlike the convergence of natural microbial community, those experiments carried out in a sterile laboratory environment only involves the loss of species. Specifically, the main manifestation of convergence in the synthetic community containing stably coexisting species lies in that the relative proportion of species tend to become consistent [7, 8]. Nonetheless, synthetic community opens a window for us to investigate the ecological mechanism. Previous studies of synthetic communities have revealed that the convergence of bacterial community can be regulated by pH [11], mortality [12], and particularly nutrient availability [13, 14]. Most existing studies focus on the changes in species proportions, but there is a lack of in-depth understanding of the gene expression changes driven by the community species interaction.In this study, we constructed a synthetic community with two model microorganisms, Escherichia coli K-12 (EC) and Pseudomonas putida KT2440 (PP), and reproduced a convergent community assembly in closed broth-culture system. In monocultures, the growth curves of both E. coli and P. putida fitted well with the bacterial growth model, and fell into a logarithmic phase at the first 4 h of bacterium culture and a stationary phase at subsequent 20 h (after the first 4 h) (Fig. 1a). When same quantities of bacteria were grown in cocultures, their quantities were basically similar to those in monocultures, particularly in the logarithmic phase (Fig. 1b–d). By contrast, the quantities of minority species in cocultures continued to increase, and they were close to the quantities in monocultures at 24 h post co-cultivation (Fig. 1b–d). Besides, statistical analysis showed that the quantities of P. putida in all three cocultures were overall greater than that in monoculture, while E. coli quantities were no more than its monoculture (Fig. 1b–d), suggesting that P. putida has a negative effect on the growth of E. coli, but E. coli promotes that of P. putida.Fig. 1: Convergence of community structure and gene expression.a–d Growth curves of E. coli and P. putida in monoculture (a) and the “1:1000”, “1:1”, “1000:1” cocultures (b–d). In b–d subplots, the growth curves of monocultures were placed on the background layer (dashed lines), and the significant differences in cell quantity between coculture and corresponding monoculture were shown (ns, non-significant; *p  More

Thank you for visiting nature.com. You are using a browser version with limited support for CSS. To obtain
the best experience, we recommend you use a more up to date browser (or turn off compatibility mode in
Internet Explorer). In the meantime, to ensure continued support, we are displaying the site without styles
and JavaScript. More