Philippa Kaur

Sample collectionThe 267 samples were collected between 2009 and 2018 from 41 locations expanding from the subarctic region to the tropics (Fig. 1, Auxillary Table S1)10 and processed using the same analytical pipeline (Fig. 2). The majority of the samples were collected using a depth-discrete Limnos tube-sampler (Limnos, Poland), with the exception of the samples from La Plata reservoir (Puerto Rico), which were collected using horizontal Van Dorn sampler (5 L capacity) and samples from Lake Loclat, which were collected using a deployed PVC-inlet connected to a peristaltic pump via tubing. Of all the lakes, 29 were sampled during the open water season and the majority of the lakes were sampled once. For 12 of the lakes only surface samples taken during the ice-covered period in winter were available, and one of the Swedish lakes (Lake Lomtjärnan) was sampled twice during the ice-covered period. Moreover, a total of 5 samples (one depth profile) from the time series of the Swiss lake (Loclat) were taken from under the ice. Time series samples were taken for Lake Loclat (seven time points, Auxillary Table S1)10 and for Lake Mekkojärvi (22 time points, see Saarenheimo et al.11 for details). For most lakes and ponds, samples were collected from multiple depths, including samples from the oxic surface layer (epilimnion), the layer with steepest change in oxygen concentration and temperature (metalimnion) and from the layer where oxygen levels were below the detection limit (hypolimnion). The exception to this were the 12 Swedish lakes sampled during ice-covered period, and five shallow ponds in Canada, for which only one sample from the oxic surface layer was taken (see Auxillary Table S1)10.Fig. 2Overview of the workflow from sample collection to mOTUs.Full size imageFrom two of the lakes, Lake Lomtjärnan in Sweden and Lake Alinen Mustajärvi in Finland, samples were collected also for single cell sorting. From both locations samples were preserved in glycerol-TE (gly-TE) and from Lomtjärnan samples were preserved also using phosphate buffered saline (PBS). For both preservants, the samples were flash frozen in liquid nitrogen after first incubating for 1 minute at ambient temperature.Simultaneous to collection of the DNA samples, also samples for environmental variables were taken. Variables included temperature, pH, conductivity, oxygen, total and dissolved nutrients (P and N species), gases (CO2 or dissolved inorganic carbon and methane (CH4)), total or dissolved organic carbon, iron, sulfate and chlorophyll a (Auxillary Table S1 and Auxillary Table S210 for the methods). As the samples were collected during multiple years and by different research groups, there was some variation for the procedures between the different sampling occasions, leading to variation in the final set of environmental data across the samples.DNA extraction and metagenome sequencingMost of the DNA samples were collected on 0.2 µm Sterivex filters (Millipore), except for the time-series samples collected from Loclat, which were collected by vacuum filtration onto 47 mm polycarbonate membrane filters with 0.2 μm pore size, and time series samples from Finnish Lake Mekkojärvi, for which the water for DNA extraction was collected from epilimnion (0–0.5 m), metalimnion (0.5–1 m) and hypolimnion (1–3 m) and pooled samples from each stratum were stored in 100 ml plastic containers and frozen at −20 °C and eventually freeze-dried (Alpha 1–4 LD plus, Christ). For all filter samples, water was filtered until the filter clogged. All filters were stored frozen (−20 to −80 °C) until the extraction of DNA. For all samples, DNA was extracted using PowerSoil DNA extraction kit (MoBio, Carlsbad, CA, USA) following the manufacturer’s instructions and the DNA concentrations were measured using Qubit dsDNA HS kit (Thermo Fisher Scientific Inc.).Sequencing libraries were prepared from 10 or 20 ng of DNA using the ThruPLEX DNA-seq Prep Kit according to the manufacturer’s preparation guide. Briefly, the DNA was fragmented using a Covaris E220 system, aiming at 400 bp fragments. The ends of the fragments were end-repaired and stem-loop adapters were ligated to the 5′ ends of the fragments. The 3′ end of the stem loop were subsequently extended to close the nick. Finally, the fragments were amplified and unique index sequences were introduced using 7 cycles of PCR followed by purification using AMPure XP beads (Beckman Coulter).The quality of the libraries was evaluated using the Agilent Fragment Analyzer system and the DNF-910-kit. The adapter-ligated fragments were quantified by qPCR using the Library quantification kit for Illumina (KAPA Biosystems/Roche) on a CFX384Touch instrument (BioRad) prior to cluster generation and sequencing.The sequencing libraries were pooled and subjected to cluster generation and paired-end sequencing with 150 bp read length S2/S4 flow-cells and the NovaSeq 6000 system (Illumina Inc.) using the v1 chemistry according to the manufacturer’s protocols. Negative controls were included to the sequencing as well as 1% of PhiX control library as a positive control.Base calling was done on the instrument by RTA (v3.3.3, 3.3.5, 3.4.4) and the resulting.bcl files were demultiplexed and converted to fastq format with tools provided by Illumina Inc., allowing for one mismatch in the index sequence. Additional statistics on sequence quality were compiled with an in-house script from the fastq-files, RTA and CASAVA output files. Sequencing was performed by the SNP&SEQ Technology Platform in Uppsala, Sweden.Single-cell sorting and DNA amplificationAll Gly-TE cryopreserved samples were thawed and diluted in 1 xPBS if needed while all plates with PBS were UV-treated with a dose of 2 J prior to sorting. Samples collected from both lakes were sorted, and then screened for organisms belonging to candidate phyla radiation. Samples collected from Lake Lomtjärnan were additionally subjected to sorting based on autofluorescence to identify and sequence cells belonging to lineage Chlorobia.For obtaining SAGs from representatives of the candidate phyla radiation (CPR), samples were first stained with 1 x SYBR Green I for approximate 30 minutes. Subsequent single cell sorting was performed with a MoFlo Astrios EQ (Beckman Coulter, USA) cell sorter using a 488 nm laser for excitation, 70 µm nozzle, sheath pressure of 60 psi and 0.1 µm sterile filtered 1x PBS as sheath fluid. Individual cells were deposited into empty 384-well plates (Biorad, CA USA) UVed at 2 Joules using a CyCloneTM robotic arm and the most stringent single cell sort settings (single mode, 0.5 drop envelope). Green fluorescence (488–530/40) was used as trigger and sort decisions were made based on combined gates of 488–530/40 Height log vs 488–530/40 Area log and 488–530/40 Height log vs SSC with increasing side scatter divided up in three different regions. Flow sorting data was interpreted and displayed using the associated software Summit v 6.3.1. Next, individual cells were subject to lysis, neutralization and whole genome amplification using MDA based on the protocol and workflow described by Rinke et al.12 but with several modifications. Reagent mastermixes were added using the MANTIS liquid dispenser (Formulatrix) and the LV or HV silicone chips. The lysozyme, D2 buffer, stop solution and MDA-mastermix were each dispensed with its own chip. Most MDA-reactions were run using the phi29 from ThermoFisher but a few were run with a more heat-stable phi29, EquiPhi also provided by ThermoFisher. The MDA reaction was carried out in a total volume of 5.2 µl. Thawed, sorted cells were first pre-treated with 400 nl/well of 12 U/µl of Ready-Lyse™ Lysozyme Solution (R1804M, Lucigen) at room temperature for 15 minutes before adding 400 nl Qiagen lysis buffer D2 followed by incubation at 95 °C for 10 seconds and 10 minutes on ice. Reactions were neutralized by adding 400 nl Qiagen Stop solution. Four µl of MDA mix containing 1x reaction buffer, 0.4 mM dNTP, 0.05 mM exonuclease-resistant Hexamers, 10 mM DTT, 1.7 U phi29 DNA polymerase (ThermoFisher Scientific) and 0.5 µM Syto13 was added to a final reaction volume of 5.2 µl. All reagents except SYTO13 were UV decontaminated at 2 Joules in a UV crosslinker. The whole genome amplification was run at 30 °C for 7 or 10 h followed by an inactivation step at 65 °C for 5 min. The reaction was monitored in real time by detection of SYTO13 fluorescence every 15 minutes using a FLUOstar® Omega plate reader (BMG Labtech, Germany) or a qPCR instrument. The EquiPhi protocol was run as previously described for ThermoFisher phi29 with the following exceptions; the EquiPhi polymerase was added in 1U/reaction, reaction buffer included with the polymerase was used and the reaction was carried out at 45 °C. The single amplified genome (SAG) DNA was stored at −20 °C until further PCR screening, library preparation and Illumina sequencing.The CPR SAGs were screened using the bacterial PCR primers targeting the 16 S rRNA gene, Bact_341 F and Bact_805 R13. The reactions were run in a LightCycler 480 PCR machine (ROCHE, MA USA) in 10 µl and a final concentration of 1 x LightCycler480 SYBR Green I Master mix, 0.25 µM of each primer and 2 µl of 60 to 80 times diluted SAGs. Following a 3 min denaturation at 95 °C, targets were amplified for 40 cycles of 95 °C for 10 s, 55 °C for 20 s, 72 °C for 30 s and a final 10 min extension at 72 °C followed by melting curve analysis. The products were purified using the NucleoSpin Gel and PCR clean-up purification kit (Macherey-Nagel, Germany), quantified using the Quant-iT TM PicoGreen® dsDNA assay kit (Invitrogen, MA USA) in a FLUOstar® Omega microplate reader (BMG Labtech, Germany) and submitted for identification by Sanger sequencing at Eurofin Genomics. All SAGs were further screened using the newly designed primers targeting the phylum Parcubacteria 684F-OD1 (3′ GTAGKRRTRAAATSCGTT 5′) and 784 R (5′ TAMNVGGGTATCTAATCC -3′). These primers target with good specificity 67% of Parcubacteria in the SILVA database14. Parcu-PCR was run at 3 min at 95 °C, 40 cycles of 95 °C for 10 s, 55 °C for 20 s, 72 °C for 30 s and a final 10 min extension at 72 °C followed by melting curve analysis. The products were purified using the NucleoSpin Gel and PCR clean-up purification kit (Macherey-Nagel, Germany), quantified using the Quant-iT TM PicoGreen® dsDNA assay kit (Invitrogen, MA USA) in a FLUOstar® Omega microplate reader (BMG Labtech, Germany) and submitted for identification by Sanger sequencing at Eurofin Genomics.To recover Chlorobia single amplified genomes, sorting was done in 2016 on a MoFlo™ Astrios EQ sorter (Beckman Coulter, USA) using a 488 and 532 nm laser for excitation, a 70 μm nozzle, a sheath pressure of 60 psi, and 0.1 μm filtered 1x PBS as sheath fluid. An ND filter ND = 1 and the masks M1 and M2 were used. The trigger channel was set to the forward scatter (FSC) at a threshold of 0.025% and sort regions were defined on autofluorescence using laser 532 nm and band pass filters 710/45 and 664/22. Three populations were sorted based on differences in autofluorescence signals. The sort mode was set to single cell with a drop envelope of 0.5. The target populations were sorted at approximately 400 events per second into 96-well plates containing 1 µl 1x PBS per well with either 1 or 10 cells (positive control) deposited. A few wells remained empty (no cell sorted) were kept as negative controls. Sorted plates were stored frozen at −80 °C.The subsequent whole genome amplification was performed in 2018 using the REPLI-g Single Cell kit (QIAGEN) following the instructions provided by the manufacturer but with total reaction volume reduced to 12.5 µl. The denaturation reagent D2, stop solution, water, and reagent tubes and strips were UV-treated at 2.5 J. The lysis was changed slightly to 10 min at 65 °C, followed by 5 min on ice before adding the stop solution. To the master mix containing water, reaction buffer, and the DNA the polymerase we added SYTO 13 (Invitrogen) at a final concentration of 0.5 µM. The amplification was performed at 30 °C for 8 hours in a plate reader with fluorescence readings every 15 min. The reaction was stopped by incubating it for 5 min at 65 °C. The plate was stored for less than a week at −20 °C. Amplified DNA was mixed thoroughly by pipetting up and down 20 times before diluting it 50x and 100x in nuclease-free water. The DNA was screened for bacterial 16 S rRNA applying the primers Bact_341 F (5′- CCTACGGGNGGCWGCAG- 3′) and Bact_805 R (5′- GACTACHVGGGTATCTAATCC-3′)13 using the LightCycler® 480 SYBR Green I Master (Roche) kit. The PCR mix contained 1.5 µl diluted amplified DNA, 1x the LightCycler® 480 SYBR Green I Master mix, 0.25 µM of each primer, and nuclease-free water in a total reaction volume of 10 µl. The PCR cycling (5 min at 95 °C, followed by 40 cycles of 10 sec at 95 °C, 20 sec at 60 °C, 30 sec at 72 °C) was followed by meltcurve analysis on the LightCycler® 480 Instrument (Roche). DNA of confirmed Chlorobia was sent to sequencing as outlined below.Library preparation and Illumina sequencing of the single cellsFor the CPR-targeted analysis, Illumina libraries were prepared from sixty SAGs mainly selected from the screening procedure in a PCR-free workflow using the sparQ DNA Frag & Library Prep Kit (Quantabio) and IDT for Illumina TruSeq UD Indexes (Illumina). Libraries were prepared from 50–250 ng of MDA-products in 25% of the recommended reaction volumes according to manufacturer’s instructions. The MDA-products were fragmented for 7 minutes (5 minutes for 4 samples) without using the DNA Frag Enhancer Solution. Library insert sizes were determined using Bioanalyzer High Sensitivity DNA Kit (Agilent). Each library was quantified using the KAPA Library Quantification kit (Roche) in 5 µl reaction volumes in a 384-well plate run on LightCycler 480 (Roche) to allow equimolar pooling before sequencing on Illumina HiSeqX v2.5 PE 2 × 150 bp including negative and positive (PhiX) controls.For the Chlorobia-targeted sequencing, amplified DNA from 23 SAGs were quantified individually with Qubit dsDNA HS assay kit (ThermoFisher Scientific) and diluted to 0.2 ng/ul in nuclease free water. Sequencing libraries were prepared with Nextera XT DNA Library Preparation Kit and combinatorial combinations of molecular identifiers in the Nextera XT Index Kit (Illumina, CA USA) according to manufacturer’s instructions. Libraries with an average length of 1200 bp were quantified with Qubit dsDNA HS assay kit to allow pooling of equal amounts of the libraries based on mass. The libraries were sequenced on an Illumina MiSeq v3 PE 2 × 300 bp including negative and positive (PhiX) controls.Data processing of the metagenome and single cell sequencesThe metagenome sequencing resulted in a total of ~107 paired-end reads of length 2 × 150 bp, amounting to a total of total 3 Tbp. The raw data was trimmed using Trimmomatic (version 0.36; parameters: ILLUMINACLIP:TruSeq 3-PE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36)15 (Auxillary Table S3)10. The trimmed data was assembled using Megahit (version 1.1.13)16 with default settings. Two types of assemblies were done, single sample assemblies for all the samples individually and a total of 53, mainly lake-wise, co-assemblies (see Auxillary Table S4)10, some samples of the Canadian ponds have also been coassembled with previously sequenced libraries of the same sample (see Auxillary Table S5)10. The relevant quality controlled reads were mapped to all the assemblies using BBmap17 with default settings and the mapping results were used to bin the contigs using Metabat (version 2.12.1, parameters –maxP 93 –minS 50 -m 1500 -s 10000)18. Genes of obtained bins were predicted and annotated using Prokka (version 1.13.3)19 using standard parameters except for the bin containing all the unbinned contigs where the –metagenome flag was used. Single-cell libraries were processed similarly to the metagenomes, but without the binning step, and using the single-cell variant of the SPAdes20 assembler instead of Megahit.Prokaryotic completeness and redundancy of all bins from Metabat and for all assembled single cells were computed using CheckM (version 1.0.13)21 (Auxillary Tables S6 and S7 for MAGs and SAGs, respectively)10. Average Nucleotide Identity (ANI) for all bin-pairs was computed with fastANI (version 1.3)22. The bins were clustered into metagenomic Operational Taxonomic Units (mOTUs) starting with 40% complete genomes with less than 5% contamination. Genome pairs with ANI above 95% were clustered into connected components. Additionally, less complete genomes were recruited to the mOTU if its ANI similarity was above 95%. Bins were taxonomically annotated in a two-step process. GTDB-Tk (version 102 with database release 89)23 was used first with default settings. Using this classification an lca database for SourMASH (version 1.0)24 was made. This database as well as one based on the GTDB release 89 was then used with SourMASH’s lca classifier for a second round of classification of bins that were not annotated with GTDB-tk (Auxillary Table S8)10.The taxonomic diversity of the bacterial (Fig. 3) and archaeal (Fig. 4) mOTUs, respectively, were visualized in a tree format. The trees were computed using GTDB-tk with one representative MAG per mOTU of the stratfreshDB, and one random representative genome per family of the GTDB. Trees were visualized using anvi’o25.Fig. 3Bacterial diversity of the stratfreshDB27. The insert illustrates the quality of the MAGs and SAGs included in the tree. Interactive version of the tree with more information available at https://anvi-server.org/moritzbuck/bacterial_diversity_of_the_stratfreshdb.Full size imageFig. 4Archaeal diversity of the stratfreshDB27. The insert illustrates the quality of the MAGs included in the tree. Interactive version of the tree with more information available at https://anvi-server.org/moritzbuck/archaeal_diversity_of_the_stratfreshdb.Full size image More

1.East, E. M. The role of reproduction in evolution. Am. Nat. https://doi.org/10.1086/279670 (1918).Article 

Google Scholar 
2.Proctor, M., Yeo, P. F. & Lack, A. A Natural History of Pollination. (1996).3.Spigler, R. B. & Ashman, T.-L. Gynodioecy to dioecy: are we there yet?. Ann. Bot. 109, 531–543. https://doi.org/10.1093/aob/mcr170%JAnnalsofBotany (2011).Article 
PubMed 
PubMed Central 

Google Scholar 
4.Barrett, S. Sexual interference of the floral kind. Heredity 88, 154–159 (2002).CAS 
PubMed 
Article 

Google Scholar 
5.Li, Q.-J. et al. Flexible style that encourages outcrossing. Nature 410, 432–432. https://doi.org/10.1038/35068635 (2001).ADS 
CAS 
Article 
PubMed 

Google Scholar 
6.Sun, S., Gao, J. Y., Liao, W. J., Li, Q. J. & Zhang, D. Y. Adaptive significance of flexistyly in Alpinia blepharocalyx (Zingiberaceae): a hand-pollination experiment. Ann. Bot. 99, 661–666 (2007).PubMed 
PubMed Central 
Article 

Google Scholar 
7.Kumar, B. D., Deepika, D. S. & Raju, A. S. Reproductive ecology of the semi-evergreen tree Vitex negundo (Lamiaceae). Phytol. Balcanica 23, 39–53 (2017).
Google Scholar 
8.Faegri, K. & Van Der Pijl, L. Principles of Pollination Ecology (Elsevier, Amsterdam, 2013).
Google Scholar 
9.Darwin, C. The Effects of Cross and Self Fertilisation in the Vegetable Kingdom (D. Appleton, Boston, 1877).
Google Scholar 
10.Baker, H. G. in Cold Spring Harbor Symposia on Quantitative Biology. 177–191 (Cold Spring Harbor Laboratory Press).11.Heithaus, E. R., Opler, P. A. & Baker, H. G. Bat activity and pollination of Bauhinia pauletia: plant-pollinator coevolution. Ecology 55, 412–419 (1974).Article 

Google Scholar 
12.Armbruster, W. S. Can indirect selection and genetic context contribute to trait diversification? A transition-probability study of blossom-colour evolution in two genera. J. Evolut. Biol. 15, 468–486 (2002).Article 

Google Scholar 
13.Bradshaw, H. Jr. & Schemske, D. W. J. N. Allele substitution at a flower colour locus produces a pollinator shift in monkeyflowers. Nature 426, 176 (2003).ADS 
CAS 
PubMed 
Article 

Google Scholar 
14.Gómez, J. M. & Zamora, R. Ecological factors that promote the evolution of generalization in pollination systems, in Plant–pollinator interactions: from specialization to generalization. 145–166 (2006).15.Barrett, S. C. & Harder, L. D. Ecology and evolution of plant mating. Trends Ecol. Evolut. 11, 73–79 (1996).CAS 
Article 

Google Scholar 
16.Elzinga, J. A. et al. Time after time: flowering phenology and biotic interactions. Trends Ecol. Evol. 22, 432–439 (2007).PubMed 
Article 

Google Scholar 
17.Huang, S.-Q., Xiong, Y.-Z. & Barrett, S. C. H. Experimental evidence of insect pollination in Juncaceae, a primarily wind-pollinated family. Int. J. Plant Sci. 174, 1219–1228. https://doi.org/10.1086/673247 (2013).Article 

Google Scholar 
18.Sánchez-Bayo, F. & Wyckhuys, K. A. Worldwide decline of the entomofauna: a review of its drivers. Biol. Conserv. 232, 8–27 (2019).Article 

Google Scholar 
19.Memmott, J., Craze, P. G., Waser, N. M. & Price, M. V. Global warming and the disruption of plant–pollinator interactions. Ecol. Lett. 10, 710–717. https://doi.org/10.1111/j.1461-0248.2007.01061.x (2007).Article 
PubMed 

Google Scholar 
20.Winfree, R., Griswold, T. & Kremen, C. Effect of human disturbance on bee communities in a forested ecosystem. Conserv. Biol. 21, 213–223. https://doi.org/10.1111/j.1523-1739.2006.00574.x (2007).Article 
PubMed 

Google Scholar 
21.Biesmeijer, J. C. et al. Parallel declines in pollinators and insect-pollinated plants in Britain and the Netherlands. Science 313, 351–354 (2006).ADS 
CAS 
PubMed 
Article 

Google Scholar 
22.Ollerton, J., Winfree, R. & Tarrant, S. How many flowering plants are pollinated by animals?. Oikos 120, 321–326. https://doi.org/10.1111/j.1600-0706.2010.18644.x (2011).Article 

Google Scholar 
23.Gallai, N., Salles, J.-M., Settele, J. & Vaissière, B. E. Economic valuation of the vulnerability of world agriculture confronted with pollinator decline. Ecol. Econ. 68, 810–821. https://doi.org/10.1016/j.ecolecon.2008.06.014 (2009).Article 

Google Scholar 
24.Mayer, C. et al. Pollination ecology in the 21st century: key questions for future research. J. Pollinat. Ecol. 3, 8–23 (2011).Article 

Google Scholar 
25.Chavez, D. J. & Lyrene, P. M. Effects of self-pollination and cross-pollination of Vaccinium darrowii (Ericaceae) and other low-chill blueberries. Hortsci. Publ. Am. Soc. Hortic. Sci. 44, 1538–1541 (2009).
Google Scholar 
26.Negussie, A., Achten, W. M. J., Verboven, H. A. F., Hermy, M. & Muys, B. Floral display and effects of natural and artificial pollination on fruiting and seed yield of the tropical biofuel crop Jatropha curcas L. Global Change Biol. Bioenergy 6, 210–218 (2014).Article 

Google Scholar 
27.Okubo, S., Yamada, M., Yamaura, T. & Akita, T. Effects of the pistil size and self-incompatibility on fruit production in Curculigo latifolia (Liliaceae). J. Jpn. Soc. Hortic. Sci. 79, 354–359 (2010).Article 

Google Scholar 
28.Benjamin, F. E. & Winfree, R. Lack of pollinators limits fruit production in commercial blueberry (Vaccinium corymbosum). Environ. Entomol. 43, 1574–1583 (2014).PubMed 
Article 

Google Scholar 
29.Bennett, J. et al. A review of European studies on pollination networks and pollen limitation, and a case study designed to fill in a gap. AoB PLANTS https://doi.org/10.1093/aobpla/ply068 (2018).Article 
PubMed 
PubMed Central 

Google Scholar 
30.Wang, H., Matsushita, M., Tomaru, N., Nakagawa, M. & Arroyo, J. Differences in female reproductive success between female and hermaphrodite individuals in the subdioecious shrub Eurya japonica (Theaceae). Plant Biol. 17, 194–200 (2015).CAS 
PubMed 
Article 

Google Scholar 
31.Wang, H., Matsushita, M., Tomaru, N. & Nakagawa, M. High male fertility in males of a subdioecious shrub in hand-pollinated crosses. AoB PLANTS 8, plw067 (2016).PubMed 
PubMed Central 
Article 

Google Scholar 
32.Shou, C., Wang, J., Zheng, X. & Guo, D. Inhibitory effect of jujuboside A on penicillin sodium induced hyperactivity in rat hippocampal CA1 area in vitro. Acta Pharmacol. Sin. 22, 986–990 (2001).CAS 
PubMed 

Google Scholar 
33.Zhang, M. et al. Inhibitory effect of jujuboside A on glutamate-mediated excitatory signal pathway in hippocampus. Planta Med. 69, 692–695 (2003).CAS 
PubMed 
Article 

Google Scholar 
34.Yue, Y. et al. Wild jujube polysaccharides protect against experimental inflammatory bowel disease by enabling enhanced intestinal barrier function. Food Funct. 6, 2568–2577 (2015).CAS 
PubMed 
Article 

Google Scholar 
35.Han, D. et al. Jujuboside A protects H9C2 cells from isoproterenol-induced injury via activating PI3K/Akt/mTOR signaling pathway. Evidence-Based Complementary Alternative Medicine 2016 (2016).36.Lu, J., Liu, M., Mao, Y. & Shen, L. Effects of vesicular-arbuscular mycorrhizae on the drought resistance of wild jujube (Zizyphs spinosus Hu) seedlings. Front. Agric. China 1, 468–471 (2007).Article 

Google Scholar 
37.Zhang, S. et al. Threshold effects of photosynthetic efficiency parameters of wild jujube in response to soil moisture variation on shell beach ridges, Shandong, China. Plant Biosyst. 148, 140–149 (2014).Article 

Google Scholar 
38.Wang, Q. Y. The developments of embryo and endosperm of Zizyphus jujuba mill. J. Integr. Plant Biol. 25, 32–37 (1983).ADS 

Google Scholar 
39.Cruden, R. W. Pollen-ovule ratios: a conservative indicator of breeding systems in flowering plants. Evolution 31, 32–46. https://doi.org/10.1111/j.1558-5646.1977.tb00979.x (1977).Article 
PubMed 

Google Scholar 
40.Dafni, A. Pollination ecology: A practical approach. (1992).41.Barrett, S. C. H. The evolution of mating strategies in flowering plants. Trends Plant Sci. 3, 335–341. https://doi.org/10.1016/s1360-1385(98)01299-0 (1998).Article 

Google Scholar 
42.Carr, D. E. & Dudash, M. R. Recent approaches into the genetic basis of inbreeding depression in plants. Philos. Trans. R. Soc. Lond. Ser. B: Biol. Sci. 358, 1071–1084 (2003).CAS 
Article 

Google Scholar 
43.Lloyd, D. G. & Webb, C. The avoidance of interference between the presentation of pollen and stigmas in angiosperms I. Dichogamy. NZ J. Bot. 24, 135–162 (1986).Article 

Google Scholar 
44.Ren, M. Stamen movements in hermaphroditic flowers: diversity and adaptive significance. J. Plant Ecol. (Chin. Vers.) 34, 867–875 (2010).
Google Scholar 
45.Xiao, C.-L. et al. Sequential stamen maturation and movement in a protandrous herb: mechanisms increasing pollination efficiency and reducing sexual interference. AoB PLANTS 9, plx019. https://doi.org/10.1093/aobpla/plx019 (2017).Article 
PubMed 
PubMed Central 

Google Scholar 
46.Nagy, E. S., Strong, L. & Galloway, L. F. Contribution of delayed autonomous selfing to reproductive success in mountain laurel, Kalmia latifolia (Ericaceae). Am. Midl. Nat. 142, 39–47 (1999).Article 

Google Scholar 
47.Liu, K.-W. et al. Pollination: self-fertilization strategy in an orchid. Nature 441, 945 (2006).ADS 
CAS 
PubMed 
Article 

Google Scholar 
48.Ye, Z.-M., Jin, X.-F., Yang, J., Wang, Q.-F. & Yang, C.-F. Accurate position exchange of stamen and stigma by movement in opposite direction resolves the herkogamy dilemma in a protandrous plant, Ajuga decumbens (Labiatae). AoB PLANTS https://doi.org/10.1093/aobpla/plz052 (2019).Article 
PubMed 
PubMed Central 

Google Scholar 
49.Brantjes, N. & De Vos, O. The explosive release of pollen in flowers of Hyptis (Lamiaceae). New Phytol. 87, 425–430 (1981).Article 

Google Scholar 
50.Guerrina, M., Casazza, G., Conti, E., Macrì, C. & Minuto, L. Reproductive biology of an Alpic paleo-endemic in a changing climate. J. Plant. Res. 129, 477–485 (2016).PubMed 
Article 

Google Scholar 
51.Bawa, K. S. & Beach, J. H. Evolution of sexual systems in flowering plants. Ann. Mo. Bot. Garden 68, 254–274 (1981).Article 

Google Scholar 
52.Dietzsch, A. C., Stanley, D. A. & Stout, J. C. Relative abundance of an invasive alien plant affects native pollination processes. Oecologia 167, 469–479 (2011).ADS 
PubMed 
Article 

Google Scholar 
53.Ollerton, J. The evolution of pollinator-plant relationships within the arthropods. Evolution and phylogeny of the arthropoda. Entomology Society of Aragon, Zaragoza, 741–758 (1999).54.Blaauw, B. R. & Isaacs, R. Flower plantings increase wild bee abundance and the pollination services provided to a pollination-dependent crop. J. Appl. Ecol. 51, 890–898 (2014).Article 

Google Scholar 
55.Inouye, D. W., Larson, B. M., Ssymank, A. & Kevan, P. G. Flies and flowers III: ecology of foraging and pollination. J. Pollinat. Ecol. 16, 115–133 (2015).
Google Scholar 
56.Rader, R., Bartomeus, I., Garibaldi, L. A., Garratt, M. P. D. & Woyciechowski, M. Non-bee insects are important contributors to global crop pollination. Proc. Natl. Acad. Sci. U.S.A. 113, 146–151 (2016).ADS 
CAS 
PubMed 
Article 

Google Scholar 
57.Corbet, S. A. Pollination and the weather. Isr. J. Plant Sci. 39, 13–30 (1990).
Google Scholar 
58.Tuell, J. K. & Isaacs, R. Weather during bloom affects pollination and yield of highbush blueberry. J. Econ. Entomol. 103, 557–562 (2010).PubMed 
Article 

Google Scholar 
59.Ellis, C. R., Feltham, H., Park, K., Hanley, N. & Goulson, D. Seasonal complementary in pollinators of soft-fruit crops. Basic Appl. Ecol. 19, 45–55 (2017).Article 

Google Scholar 
60.Wang, W., Liu, Y., Chen, F.-D. & Dai, H.-G. Behavior and activity rhythm of flower-visiting insects on Chrysanthemum morifolium in Nanjing suburb. Shengtaixue Zazhi 27, 1167–1172 (2008).
Google Scholar 
61.Waser, N. M., Chittka, L., Price, M. V., Williams, N. M. & Ollerton, J. Generalization in pollination systems, and why it matters. Ecology 77, 1043–1060 (1996).Article 

Google Scholar 
62.Navarro-Pérez, M., López, J., Rodríguez-Riaño, T. & Ortega-Olivencia, A. Reproductive system of two Mediterranean Scrophularia species with large, showy flowers. Bot. Lett. 166, 467–477 (2019).Article 

Google Scholar 
63.Elle, E. Floral adaptations and biotic and abiotic selection pressures. Plant adaptation: Molecular genetics and ecology. National Research Council of Canada, Ottawa, Ontario, 111–118 (2004).64.Redmond, C. M. & Stout, J. C. Breeding system and pollination ecology of a potentially invasive alien Clematis vitalba L. Ireland. J. Plant Ecol. 11, 56–63. https://doi.org/10.1093/jpe/rtw137%JJournalofPlantEcology (2018).Article 

Google Scholar 
65.Byers, D. & Waller, D. Do plant populations purge their genetic load? Effects of population size and mating history on inbreeding depression. Ann. Rev. Ecol. Syst. 30, 479–513 (1999).Article 

Google Scholar 
66.Lande, R. & Schemske, D. W. The evolution of self-fertilization and inbreeding depression in plants I. Genetic models. Evolution 39, 24–40 (1985).PubMed 
Article 

Google Scholar 
67.Zhang, C. et al. The genetic basis of inbreeding depression in potato. Nat. Genet. 51, 374–378 (2019).CAS 
PubMed 
Article 

Google Scholar 
68.Jones, C. E. & Little, R. J. Handbook of Experimental Pollination Biology (Scientific and Academic Editions, New York, 1983).
Google Scholar 
69.Liu, P. et al. Study on the biological basis of pollination in Chinese Jujube (Zizyphus jujuba) and Wild Jujube (Z. spinosa). Journal of Fruit Science 21(3), 224–228 (2004).
Google Scholar 
70.Sun, Y., Wu, C., Wang, D. & Wang, Z. Comparative Study on Floral Organ Structure, Pollen Morphology and Viability of Ziziphus acdiojujuba. Chin. Agric. Sci. Bull. 32(4), 87–91 (2016).
Google Scholar  More

To enhance the legibility of the results, the habituation phase was marked as the H mean score from habituation trials 5 to 7, which served as a reference value for further analyses, while the test trials were marked as T1, T2, and T3, respectively. Novelty, i.e. addition or change of objects in zone C, was introduced in the first test trial T1.The initial four habituation trials have not been presented here, as they served only as a habituation phase and not as an element of the comparative analysis of the animals’ response to novelty.The data was analysed using a General Linear Model procedure GLM, with repeated measurements H, T1, T2, T3 as within-subject factors, followed by an LSD PostHoc test which involved a comparison of the habituation phase H with the three test trials T1, T2 and T3. Bonferroni correction for multiple comparisons was employed. Differences were considered significant for p ≤ 0.05. Data analysis was carried out using JASP v. 0.14.1 software, an open-source project supported by the University of Amsterdam.Time spent in the transporterThe amount of time spent in the transporter, excluding the latency to leave the transporter (that is, the amount of time from the moment the transporter was opened until the rat first entered the experimental apparatus), was measured for each group.In the ADD group, the analysis showed a significant main effect of trial: F(3, 39) = 5.033, p = 0.005, Eta2 = 0.279 (Wilks’ Lambda). A post-hoc analysis showed a significant decrease in the time spent in the transporter in the first and third test trials compared to the habituation phase (T1: p = 0.008, d = 1.090; T3: p = 0.017, d = 0.982).In the CMPLX group, the analysis showed a significant main effect of trial: F(3, 36) = 8.695, p  More

1.Lindenmayer, D., Cunningham, S. & Young, A. Land use intensification: Effects on agriculture, biodiversity and ecological processes (CSIRO Publishing, Collingwood, 2012).Book 

Google Scholar 
2.Gibson, D. J. Grasses and grassland ecology (Oxford University Press, Oxford, 2009).
Google Scholar 
3.White, R., Murray, S., & Rohweder, M. Pilot Analysis of Global Ecosystems: Grassland Ecosystems. (2000). https://doi.org/10.1021/es00328814.Schmidt, A. C., Fraser, L. H., Carlyle, C. N. & Bassett, E. R. L. Does cattle grazing affect ant abundance and diversity in temperate grasslands?. Rangeland Ecol. Manag. 65(3), 292–298. https://doi.org/10.2111/REM-D-11-00100.1 (2012).Article 

Google Scholar 
5.Phifer, C. C., Knowlton, J. L., Webster, C. R., Flaspohler, D. J. & Licata, J. A. Bird community responses to afforested eucalyptus plantations in the Argentine pampas. Biodivers. Conserv. 26(13), 3073–3101. https://doi.org/10.1007/s10531-016-1126-6 (2017).Article 

Google Scholar 
6.Flynn, D. F. B. et al. Loss of functional diversity under land use intensification across multiple taxa. Ecol. Lett. 12(1), 22–33. https://doi.org/10.1111/j.1461-0248.2008.01255.x (2009).Article 
PubMed 

Google Scholar 
7.Sasaki, T. et al. Nestedness and niche-based species loss in moorland plant communities. Oikos 121(11), 1783–1790. https://doi.org/10.1111/j.1600-0706.2012.20152.x (2012).Article 

Google Scholar 
8.Baselga, A. Partitioning the turnover and nestedness components of beta diversity. Glob. Ecol. Biogeogr. 19(1), 134–143. https://doi.org/10.1111/j.1466-8238.2009.00490.x (2010).Article 

Google Scholar 
9.Kraft, N. J. B. et al. Community assembly, coexistence and the environmental filtering metaphor. Funct. Ecol. 29(5), 592–599. https://doi.org/10.1111/1365-2435.12345 (2015).Article 

Google Scholar 
10.Swenson, N. G. & Enquist, J. Opposing assembly mechanisms in a Neotropical dry forest: Implications for phylogenetic and functional community ecology. Ecology 90(8), 2161–2170 (2009).Article 

Google Scholar 
11.Stubbs, W. J. & Wilson, J. B. Evidence for limiting similarity in a sand dune community. J. Ecol. 92, 557–567 (2004).Article 

Google Scholar 
12.Hooper, D. U. et al. Effects of biodiversity on ecosystem functioning: A consensus of current knowledge. America 75(1), 3–35. https://doi.org/10.1890/04-0922 (2005).Article 

Google Scholar 
13.Gámez-Virués, S. et al. Landscape simplification filters species traits and drives biotic homogenization. Nat. Commun. 6, 8568. https://doi.org/10.1038/ncomms9568 (2015).ADS 
CAS 
Article 
PubMed 
PubMed Central 

Google Scholar 
14.Díaz, S. & Cabido, M. Vive la différence: Plant functional diversity matters to ecosystem processes. Trends Ecol. Evol. 16(11), 646–655. https://doi.org/10.1016/S0169-5347(01)02283-2 (2001).Article 

Google Scholar 
15.Bruno, J. F. & Cardinale, B. J. Cascading effects of predator richness. Front. Ecol. Environ. 6(10), 539–546. https://doi.org/10.1890/070136 (2008).Article 

Google Scholar 
16.Avalos, G., Rubio, G. D., Bar, M. E. & González, A. Arañas (Arachnida: Araneae) asociadas a dos bosques degradados del Chaco húmedo en Corrientes, Argentina. Rev. Biol. Trop. 55(3–4), 899–909 (2007).PubMed 

Google Scholar 
17.Downie, I. S. et al. The impact of different agricultural land-uses on epigeal spider diversity in Scotland. J. Insect Conserv. 3(4), 273–286 (1999).Article 

Google Scholar 
18.Salas-Lopez, A., Violle, C., Mallia, L. & Orivel, J. Land-use change effects on the taxonomic and morphological trait composition of ant communities in French Guiana. Insect. Conserv. Divers. 11(2), 162–173. https://doi.org/10.1111/icad.12248 (2018).Article 

Google Scholar 
19.Mousseau, T. A. Ectotherms follow the converse to Bergmann’s rule. Evolution 51(2), 630. https://doi.org/10.2307/2411138 (1997).Article 
PubMed 

Google Scholar 
20.Woolley, C., Thomas, C. F. G., Blackshaw, R. P. & Goodacre, S. L. Aerial dispersal activity of spiders sampled from farmland in southern England. J. Arachnol. 44(3), 347–358. https://doi.org/10.1636/p15-56.1 (2016).Article 

Google Scholar 
21.Rypstra, A. L., Carter, P. E., Balfour, R. A. & Marshall, S. D. Architectural features of agricultural habitats and their impact on the spider inhabitants. J. Arachnol. 27(1), 371–377. https://doi.org/10.2307/3706009 (1999).Article 

Google Scholar 
22.Tuf, I. H., Dedek, P. & Veselý, M. Does the diurnal activity pattern of carabid beetles depend on season, ground temperature and habitat?. Arch. Biol. Sci. 64(2), 721–732. https://doi.org/10.2298/ABS1202721T (2012).Article 

Google Scholar 
23.Entling, W., Schmidt-Entling, M. H., Bacher, S., Brandl, R. & Nentwig, W. Body size-climate relationships of European spiders. J. Biogeogr. 37(3), 477–485. https://doi.org/10.1111/j.1365-2699.2009.02216.x (2010).Article 

Google Scholar 
24.Blandenier, G. Ballooning of spiders (Araneae) in Switzerland: General results from an eleven-year survey. Arachnology 14(7), 308–316. https://doi.org/10.13156/arac.2009.14.7.308 (2014).Article 

Google Scholar 
25.Greenstone, M. H. Determinants of web spider species diversity: Vegetation structural diversity vs. prey availability. Oecologia 62(3), 299–304 (1984).ADS 
Article 

Google Scholar 
26.Morello, J., Matteucci, S. D., & Rodríguez, A. F. Ecorregiones y complejos ecosistémicos de argentina. Orientación Gráfica Editora, Buenos Aires (2012).27.Satorre, E. H. Cambios tecnológicos en la agricultura argentina actua. Ciencia hoy. 15(87), 6 (2005).
Google Scholar 
28.Viglizzo, E., La Pampa, I.C.R., Satorre, E., Solbrig, O.T., Torres, F. & Ingaramo, J. The provision of ecosystem services and human well-being in the Pampas of Argentina. Millennium Ecosystem Assessment: Full Report (2005).29.INTA. Instituto Nacional de Tecnología Agropecuaria (INTA). Plan De Tecnologia Regional 2009–2011, INTA Centro Regional Entre Ríos (2009).30.Santoandré, S., Filloy, J., Zurita, G. A. & Bellocq, M. I. Ant taxonomic and functional diversity show differential response to plantation age in two contrasting biomes. For. Ecol. Manag. 437, 304–313. https://doi.org/10.1016/j.foreco.2019.01.021 (2019).Article 

Google Scholar 
31.Pinto, C. M., Santoandré, S., Zurita, G., Bellocq, M. I. & Filloy, J. Conifer plantations in grassland and subtropical forest: Does spider diversity respond different to edge effect?. J. For. Res. 23(5), 253–259. https://doi.org/10.1080/13416979.2018.1506248 (2018).CAS 
Article 

Google Scholar 
32.Bell, J., Wheater, C. & Cullen, W. The implications of grassland and heathland management for the conservation of spider communities: A review. J. Zool. 255, 377–387. https://doi.org/10.1017/s0952836901001479 (2001).Article 

Google Scholar 
33.Spears, L.R., & MacMahon, J.A. An experimental study of spiders in a shrub-steppe ecosystem: The effects of prey availability and shrub architecture. J. Arachnol. 40(2):218–227 (2012). http://digitalcommons.usu.edu/etd/1207/34.Abràmoff, M. D., Magalhães, P. J. & Ram, S. J. Image processing with ImageJ. Biophotonics Int. 11(7), 36–42 (2004).
Google Scholar 
35.Merrett, P. & Snazell, R. A comparison of pitfall trapping and vacuum sampling for assessing spider faunas on heath-land at Ashdown Forest, south-east England. Bull. Br. Arachnol. Soc. 6(1), 1–13 (1983).
Google Scholar 
36.Lambeets, K., Vandegehuchte, M., Jean-Pierre, M. & Dries, B. Physical defences wear you down: Progressive and. J. Anim. Ecol. 78, 281–291. https://doi.org/10.1111/j.1365-2656.2007.0 (2009).Article 

Google Scholar 
37.Duelli, P., Obrist, M. K. & Schmatz, D. R. Environment Biodiversity evaluation in agricultural landscapes: Above-ground insects (Woodhead Publishing Limited, Cambridge, 1999). https://doi.org/10.1016/B978-0-444-50019-9.50006-6.Book 

Google Scholar 
38.Munévar, A., Rubio, G. D. & Zurita, G. A. Changes in spider diversity through the growth cycle of pine plantations in the semi-deciduous Atlantic forest: The role of prey availability and abiotic conditions. For. Ecol. Manag. https://doi.org/10.1016/j.foreco.2018.03.025 (2017).Article 

Google Scholar 
39.Horváth, R., Lengyel, S., Szinetár, C. & Jakab, L. L. The effect of prey availability on spider assemblages on European black pine (Pinus nigra) bark: Spatial patterns and guild structure. Can. J. Zool. 83(2), 324–335. https://doi.org/10.1139/z05-009 (2005).Article 

Google Scholar 
40.Bonte, D., Borre, J. V., Lens, L. & Maelfait, J.-P. Geographical variation in wolf spider dispersal behaviour is related to landscape structure. Anim. Behav. 72(3), 655–662. https://doi.org/10.1016/j.anbehav.2005.11.026 (2006).Article 

Google Scholar 
41.Legendre, P., Legendre, L. Numerical ecology: Developments in environmental modelling. Developments in Environmental Modelling. 20 (1998)42.R Development Core Team. R: A language and environment for statistical computing. Vienna, Austria. Internet: http://www.R-project.org. 2012.43.Oksanen, J., Blanchet, G., Kindt, R., Legendre, P., Minchin, P.R., O’hara, R.B., Simpson, G.L., Solymos, P., Stevens, M.H.H., Wagner, H. Vegan: community ecology package 2.3–2 (2015).44.Laliberté, E. & Legendre, P. A distance-based framework for measuring functional diversity from multiple traits. Ecology 91(1), 299–305 (2010).Article 

Google Scholar 
45.Lavorel, S. et al. Assessing functional diversity in the field—Methodology matters!. Funct. Ecol. 22(1), 134–147. https://doi.org/10.1111/j.1365-2435.2007.01339.x (2008).Article 

Google Scholar 
46.Leps, J., de Bello, F., Lavorel, S., Berman, S. Quantifying and interpreting functional diversity of natural communities: Practical considerations matter (2006).47.Zuur, A. F., Ieno, E. N., Walker, N. J., Saveliev, A. A. & Smith, G. M. Mixed effects models and extensions in ecology with R (Springer, Berlin, 2009).Book 

Google Scholar 
48.Santoandré, S., Filloy, J., Zurita, G. A. & Bellocq, M. I. Taxonomic and functional β-diversity of ants along tree plantation
chronosequences differ between contrasting biomes. Basic Appl. Ecol. 41, 1–12. https://doi.org/10.1016/j.baae.2019.08.004 (2019).49.Anderson, M. J. & Walsh, D. C. I. PERMANOVA, ANOSIM, and the Mantel test in the face of heterogeneous dispersions- What null hypothesis are you .pdf. Ecol. Monogr. 83(4), 557–574. https://doi.org/10.1890/12-2010.1 (2013).Article 

Google Scholar 
50.Swenson, N. G. Functional and phylogenetic ecology in R (Springer, Berlin, 2014). https://doi.org/10.1007/978-1-4614-9542-0.Article 
MATH 

Google Scholar 
51.Craven, D., Hall, J. S., Berlyn, G. P., Ashton, M. S. & van Breugel, M. Environmental filtering limits functional diversity during succession in a seasonally wet tropical secondary forest. J. Veg. Sci. 29(3), 511–520. https://doi.org/10.1111/jvs.12632 (2018).Article 

Google Scholar 
52.Woodcock, B. A., Pywell, R. F., Roy, D. B., Rose, R. J. & Bell, D. Grazing management of calcareous grasslands and its implications for the conservation of beetle communities. Biol. Cons. 125, 193–202. https://doi.org/10.1016/j.biocon.2005.03.017 (2005).Article 

Google Scholar 
53.Mangels, J., Fiedler, K., Schneider, F. D. & Blüthgen, N. Diversity and trait composition of moths respond to land-use intensification in grasslands: Generalists replace specialists. Biodivers. Conserv. 26(14), 3385–3405. https://doi.org/10.1007/s10531-017-1411-z (2017).Article 

Google Scholar 
54.Martello, F. et al. Homogenization and impoverishment of taxonomic and functional diversity of ants in Eucalyptus plantations. Sci. Rep. 8(1), 1–11. https://doi.org/10.1038/s41598-018-20823-1 (2018).ADS 
CAS 
Article 

Google Scholar 
55.Rubio, G. D., Nadal, M. F., Munévar, A. C., Avalos, G. & Perger, R. Iberá Wetlands: Diversity hotspot, valid ecoregion or transitional area? Perspective from a faunistic jumping spiders revision (Araneae: Salticidae). Species 19, 117–131 (2018).
Google Scholar 
56.Schiapelli, R. E. Arañas argentinas. Museo Argentino de Ciencias Naturales “Bernardino Rivadavia.” (1948).57.Zapata, L. & Grismando, C. Lista sistemática de arañas (Arachnida: Araneae) de la Reserva Ecológica Costanera Sur (Ciudad Autónoma de Buenos Aires, Argentina), con notas sobre su taxonomía y distribución. Rev. Mus. Argentino Cienc. Nat. 17(2), 183–211 (2015).
Google Scholar 
58.Argañaraz, C. I., Rubio, G. D. & Gleiser, R. M. Spider communities in urban green patches and their relation to local and landscape traits. Biodivers. Conserv. 27(4), 981–1009. https://doi.org/10.1007/s10531-017-1476-8 (2018).Article 

Google Scholar 
59.Bao, L., et al. Spider assemblages associated with different crop stages of irrigated rice agroecosystems from eastern Uruguay. Biodivers. Data J. (2018) (6).60.Uetz, G. W. Habitat structure and spider foraging. Habitat Struct. 1948, 325–348. https://doi.org/10.1007/978-94-011-3076-9_16 (1991).Article 

Google Scholar 
61.Balfour, R. A. & Rypstra, A. L. The influence of habitat structure on spider density in a no-till soybean agroecosystem. J. Arachnol. 26, 221–226 (1998).
Google Scholar 
62.Robinson, J. V. The effect of architectural variation in habitat on a spider community: An experimental field study. Ecol. Soc. Am. 62(1), 73–80 (1981).
Google Scholar 
63.Chisté, M. N. et al. Losers, winners, and opportunists: How grassland land-use intensity affects orthopteran communities. Ecosphere 7(11), e01545 (2016).Article 

Google Scholar 
64.Blandenier, G., Bruggisser, O. T., Rohr, R. P. & Bersier, L. F. Are phenological patterns of ballooning spiders linked to habitat characteristics?. J. Arachnol. 41(2), 126–132. https://doi.org/10.1636/P12-48 (2013).Article 

Google Scholar 
65.De Bello, F. et al. Evidence for scale- and disturbance-dependent trait assembly patterns in dry semi-natural grasslands. J. Ecol. 101(5), 1237–1244. https://doi.org/10.1111/1365-2745.12139 (2013).Article 

Google Scholar 
66.Gibb, H. et al. Habitat disturbance selects against both small and large species across varying climates. Ecography 41(7), 1184–1193. https://doi.org/10.1111/ecog.03244 (2018).Article 

Google Scholar 
67.Entling, W., Schmidt, M. H., Bacher, S., Brandl, R. & Nentwig, W. Niche properties of Central European spiders: Shading, moisture and the evolution of the habitat niche. Glob. Ecol. Biogeogr. 16(4), 440–448. https://doi.org/10.1111/j.1466-8238.2006.00305.x (2007).Article 

Google Scholar  More

SampleA total of 620 carnivore tooth pits were included in the present study. These samples included tooth marks produced by;

Brown Bears (Ursus arctos, Ursidae, 69 pits)

Spotted Hyenas (Crocuta crocuta, Hyaenidae, 86 pits)

Wolves (Canis lupus, Canidae, 80 pits)

African Wild Dogs (Lycaon pictus, Canidae, 89 pits)

Foxes (Vulpes vulpes, Canidae, 53 pits)

Jaguars (Panthera onca, Felidae, 77 pits)

Leopards (Panthera pardus, Felidae, 84 pits)

Lions (Panthera leo, Felidae, 82 pits)

Samples originated from a number of different sources, including animals kept in parks as well as wild animals. Samples obtained from wild animals included those produced by foxes as well as wolves. The only sample containing both wild and captive animals was the wolf sample. Preliminary data from these tooth pits revealed animals in captivity to have highly equivalent tooth pit morphologies to wild animals ((vert d vert ) = 0.125, p = 9.0e−14, BFB = 1.4e+11), while tooth scores revealed otherwise ((vert d vert ) = 0.152, p = 0.99, BFB = 3.7e+01 against (H_{a})). Under this premise, and so as to avoid the influence of confounding variables that go beyond the scope of the present study, tooth scores were excluded from the present samples and are under current investigation (data in preperation). Nevertheless, other research have shown tooth pits to be more informative than tooth scores when considering morphology20,23.When working with tooth mark morphologies, preference is usually given to marks found on long bone diaphyses. This is preferred considering how diaphyses are denser than epiphyses, and are thus more likely to survive during carnivore feeding. Nevertheless, when working with captive or semi-captive animals, controlling the bones that carnivores are fed is not always possible. This is due to the rules and regulations established by the institution where these animals are kept64. While this was not an issue for the majority of the animals used within the present study, in the case of P. pardus, animals were only fed ribs in articulation with other axial elements. In light of this, a careful evaluation on the effects this may have on the analogy of our samples was performed (Supplementary Appendix 2). These reflections concluded that in order to maintain a plausible analogy with tooth marks produced by other animals on diaphyses, tooth marks could only be used if found on the shaft of bovine ribs closest to the tuburcle, coinciding with the posterior and posterior-lateral portions of the rib, and farthest away from the costochondral junction65. This area of the rib corresponds to label RI3 described by Lam et al.65. Moreover, with a reported average cortical thickness of 2.3mm (± 0.13 mm) and Bone Mineral Density of (4490 kg/m^{3} [213.5, 334.6])66, bovine ribs are frequently employed in most bone simulation experiments used in agricultural as well as general surgical sciences. Finally, considering the grease, muscle and fat content of typical domestic bovine individuals67, alongside the general size of P. pardus teeth, it was concluded that the use of rib elements for this sample was the closest possible analogy to the tooth marks collected from other animals.Carnivores were fed a number of different sized animals, also dependent in most cases on the regulations established by the institution where these animals are kept64. Nevertheless, recent research has found statistical similarities between tooth marks found on different animals25, with the greatest differences occurring between large and small sized animals. Needless to say, considering the typical size of prey some of these carnivores typically consume, this factor was not considered of notable importance for the present study25 (Supplementary Appendix 1).For the purpose of comparisons, animals were split into 5 groups according to ecosystem as well as taxonomic family. From an ecological perspective, two datasets were defined; (1) the Pleistocene European Taxa dataset containing U. arctos, V. vulpes, C. crocuta, P. pardus, P. leo and C. lupus; and (2) the African Taxa dataset containing C. crocuta, P. pardus, L. pictus and P. leo. When considering taxonomic groupings, animals were separated into 3 groups, including; (1) the Canidae dataset, including V. vulpes, L. pictus and C. lupus; (2) the Felidae dataset, including P. pardus, P. onca and P. leo; and (3) a general Taxonomic Family dataset, including all Canidae in the same group, all Felidae in the same group, followed by Hyaenidae and Ursidae. Some complementary details on each of these carnivores have been included in Supplementary Appendix 1.All experiments involving carnivores were performed in accordance with the relevant ethical guidelines as set forth by park keepers and general park regulations. No animals were sacrificed specifically for the purpose of these experiments. Likewise, carnivores were not manipulated or handled at any point during the collection of samples. Collection of chewed bones were performed directly by park staff and assisted by one of the authors (JY). The present study followed the guidelines set forth by ARRIVE (https://arriveguidelines.org/) wherever necessary. No licenses or permits were required in order to perform these experiments. Finally, in the case of animals in parks, bone samples were provided by the park according to normal feeding protocols. More details can be consulted in the Extended Samples section of the supplementary files.3D modelling and landmark digitisationDigital reconstructions of tooth marks were performed using Structured Light Surface Scanning (SLSS)68. The equipment used in the present study was the DAVID SLS-2 Structured Light Surface Scanner located in the C.A.I. Archaeometry and Archaeological Analysis lab of the Complutense University of Madrid (Spain). This equipment consists of a DAVID USB CMOS Monochrome 2-Megapixel camera and ACER K11 LED projector. Both the camera and the projector were connected to a portable ASUS X550VX personal laptop (8 GB RAM, Intel® CoreTM i5 6300HQ CPU (2.3 GHz), NVIDIA GTX 950 GPU) via USB and HDMI respectively. The DAVID’s Laser Scanner Professional Edition software is stored in a USB Flash Drive. Equipment were calibrated using a 15 mm markerboard, using additional macro lenses attached to both the projector and the camera in order to obtain optimal resolution at this scale. Once calibrated the DAVID SLS-2 produces a point cloud density of up to 1.2 million points which can be exported for further processing via external software.The landmark configuration used for this study consists of a total of 30 landmarks (LMs)21; 5 fixed Type II landmarks18 and a (5 times 5) patch of semilandmarks69 (Fig. S2). Of the 5 fixed landmarks, LM1 and LM2 mark the maximal length (l) of each pit. For the correct orientation of the pit, LM1 can be considered to be the point along the maximum length furthest away from the perpendicular axis marking the maximum width (w). LM2 would therefore be the point closest to said perpendicular axis (see variables (d_{1}) and (d_{2}) in Fig. S2 for clarification). LM3 and LM4 mark the extremities of the perpendicular axis (w) with LM3 being the left-most extremity and LM4 being the right-most extremity. LM5 is the deepest point of the pit. The semilandmark patch is then positioned over the entirety of the pit, so as to capture the internal morphology of the mark.Landmark collection was performed using the free Landmark Editor software (v.3.0.0.6.) by a single experienced analyst. Inter-analyst experiments prior to landmark collection revealed the landmark model to have a robustly defined human-induced margin of error of 0.14 ± 0.09 mm (Median ± Square Root of the Biweight Midvariance). Detailed explanations as well as an instructional video on how to place both landmarks and semilandmarks can be consulted in the Supplementary Appendix and main text of Courtenay et al.21.Geometric morphometricsOnce collected, landmarks were formatted as morphologika files and imported into the R free software environment (v.3.5.3, https://www.r-project.org/). Initial processing of these files consisted in the orthogonal tangent projection into a new normalized feature space. This process, frequently referred to as Generalized Procrustes Analysis (GPA), is a valuable tool that allows for the direct comparison of landmark configurations18,19,70. GPA utilises different superimposition procedures (translation, rotation and scaling) to quantify minute displacements of individual landmarks in space71. This in turn facilitates the comparison of landmark configurations, as well as hypothesis testing, using multivariate statistical analyses. Nevertheless, considering observations made by Courtenay et al.20,21,25 revealed tooth mark size to be an important conditioning factor in their morphology, prior analyses in allometry were also performed72. From this perspective, allometric analyses first considered the calculation of centroid sizes across all individuals; the square root of the sum of squared distances of all landmarks of an object from their centroid18. These calculations were then followed by multiple regressions to assess the significance of shape-size relationships. For regression, the logarithm of centroid sizes were used. In cases where shape-size relationships proved significant, final superimposition procedures were performed excluding the scaling step of GPA (form).In addition to these analyses, preliminary tests were performed to check for the strength of phylogenetic signals73. This was used as a means of testing whether groups of carnivores produced similar tooth pits to other members of the same taxonomic family. For details on the phylogenies used during these tests, consult Fig. S1 and Supplementary Appendix 1.For the visualisation of morphological trends and variations, Thin Plate Splines (TPS) and central morphological tendencies were calculated19,71. From each of these mean landmark configurations, for ease of pattern visualisation across so many landmarks, final calculations were performed using Delaunay 2.5D Triangulation algorithms74 creating visual meshes of these configurations in Python (v.3.7.4, https://www.python.org/).Once normalised, landmark coordinates were processed using dimensionality reduction via Principal Components Analyses (PCA). In order to identify the optimal number of Principal Component Scores (PC Scores) that best represented morphological variance, permutation tests were performed calculating the observed variance explained by each PC with the permuted variance over 50 randomized iterations75. Multivariate Analysis of Variance (MANOVA) tests were then performed on these select PCs to assess the significance of multivariate morphological variance among samples.Geometric Morphometric applications were programmed in the R programming language (Sup. Appendix 8).Robust statisticsWhile GPA is known to normalize data76, this does not always hold true. Under this premise, caution must be taken when performing statistical analyses on these datasets. Taking this into consideration, prior to all hypothesis testing, normality tests were also performed. These included Shapiro tests and the inspection of Quantile–Quantile graphs. In cases where normality was detected, univariate hypothesis tests were performed using traditional parametric Analysis of Variance (ANOVA). For multivariate tests, such as MANOVA, calculations were derived using the Hotelling-Lawley test-statistic. When normality was rejected, robust alternatives to each of these tests were chosen. In the case of univariate testing, the Kruskal–Wallis non-parametric rank test was prefered, while for MANOVA calculations, Wilk’s Lambda was used.Finally, in light of some of the recommendations presented by The American Statistical Association (ASA), as debated in Volume 73, Issue Sup1 of The American Statistician77,78, the present study considers p-values of ( >2sigma ) from the mean to indicate only suggestive support for the alternative hypothesis ((H_{a})). (p ; > ; 0.005), or where possible, (3sigma ) was therefore used as a threshold to conclude that (H_{a}) is “significant”. In addition, Bayes Factor Bound (BFB) values (Eq. 1) have also been included alongside all corresponding p-Values79. Unless stated otherwise, BFBs are reported as the odds in favor of the alternative hypothesis (BFB:1). More details on BFB, Bayes Factors and the (p ; > ; 3sigma ) threshold have been included in Supplementary Appendix 3. General BFB calibrations in accordance with Benjamin and Berger’s Recommendation 0.379, as well as False Positive Risk values according to Colquhoun’s proposals80, have also been included in Table S20 of Supplementary Appendix 3.$$begin{aligned} BFB = frac{1}{-e ; p ; log (p)} end{aligned}$$
(1)
All statistical applications were programmed in the R programming language (Sup. Appendix 8).Computational learningComputational Learning employed in this study consisted of two main types of algorithm; Unsupervised and Supervised algorithms. The concept of “learning” in AI refers primarily to the creation of algorithms that are able to extract patterns from raw data (i.e. “learn”), based on their “experience” through the construction of mathematical functions38,81. The basis of all AI learning activities include the combination of multiple components, including; linear algebra, calculus, probability theory and statistics. From this, algorithms can create complex mathematical functions using many simpler concepts as building blocks38. Here we use the term “Computational Learning” to refer to a very large group of sub-disciplines and sub-sub-disciplines within AI. Deep Learning and Machine Learning are terms frequently used (and often debated), however, many more branches and types of learning exist. Under this premise, and so as to avoid complication, the present study has chosen to summarise these algorithms using the term “Computational”.Similar to the concepts of Deep and Machine Learning, many different types of supervision exist. The terms supervised and unsupervised refer to the way raw data is fed into the algorithm. In most literature, data will be referred to via the algebraic symbol x, whether this be a vector, scalar or matrix. The objective of algorithms are to find patterns among a group of x. In an unsupervised context, x is directly fed into the algorithm without further explanation. Algorithms are then forced to search for patterns that best explain the data. In the case of supervised contexts, x is associated with a label or target usually denominated as y. Here the algorithm will try and find the best means of mapping x to y. From a statistical perspective, this can be explained as (pleft( y vert x right) ). In sum, unsupervised algorithms are typically used for clustering tasks, dimensionality reduction or anomaly detection, while supervised learning is typically associated with classification tasks or regression.The workflow used in the present study begins with dimensionality reduction, as explained earlier with the use of PCA. While preliminary experiments were performed using non-linear dimensionality reduction algorithms, such as t-distributed Stochastic Neighbor Embedding (t-SNE)82 and Uniform Manifold Approximation and Projection (UMAP)83, PCA was found to be the most consistent across all datasets, a point which should be developed in detailed further research. Once dimensionality reduction had been performed, and prior to any advanced computational modelling, datasets were cleaned using unsupervised Isolation Forests (IFs)84. Once anomalies had been removed, data augmentation was performed using two different unsupervised approaches; Generative Adversarial Networks (GANs)38,39,40,41 and Markov Chain Monte Carlo (MCMC) sampling44. Data augmentation was performed for two primary reasons; (1) the simulation of larger datasets to ensure supervised algorithms have enough information to train from, and (2) to balance datasets so each sample has the same size. Both MCMCs and GANs were trialed and tested using robust statistics to evaluate quality of augmented data41. Once the best model had been determined, each of the datasets were augmented so they had a total sample size of (n = 100). In the case of the Taxonomic Family dataset, augmentation was performed until all samples had the same size as the largest sample.Once augmented, samples were used for the training of supervised classification models. Two classification models were tried and tested; Support Vector Machines (SVM)85 and Neural Support Vector Machines (NSVM)86,87. NSVMs are an extension of SVM using Neural Networks (NNs)38 as feature extractors, in substituting the kernel functions typically used in SVMs. Hyperparameter optimization for both SVMs and NSVMs were performed using Bayesian Optimization Algorithms (BOAs)88.Supervised computational applications were performed in both the R and Python programming languages (Sup. Appendix 8). For full details on both unsupervised and supervised computational algorithms, consult the Extended Methods section of the Supplementary Materials.Evaluation of supervised learning algorithms took into account a wide array of different popular evaluation metrics in machine and deep learning. These included; Accuracy, Sensitivity, Specificity, Precision, Recall, Area Under the receiver operator characteristic Curve (AUC), the F-Measure (also known as the F1 Score), Cohen’s Kappa ((kappa )) statistic, and model Loss. Each of these metrics, with the exception of loss, are calculated using confusion matrices, measuring the ratio of correctly classified individuals (True Positive & True Negative) as well as miss-classified individuals (False Positive & False Negative). For more details see Supplementary Appendix 6.Accuracy is simply reported as either a decimal (left[ 0, 1right] ) or a percentage. Accuracy is a metric often misinterpreted, as explained in Supplementary Appendix 6, and should always be considered in combination with other values, such as Sensitivity or Specificity. Both Sensitivity and Specificity are values reported as decimals (left[ 0, 1right] ), and are used to evaluate the proportion of correct classifications and miss-classifications. AUC values are derived from receiver operator characteristic curves, a method used to balance and graphically represent the rate of correctly and incorrectly classified individuals. The closest the curve gets to reaching the top left corner of the graph, the better the classifier, while diagonal lines in the graph represent a random classifier (poor model). In order to quantify the curvature of the graph, the area under the curve can be calculated (AUC), with (AUC=1) being a perfect classifier and (AUC=0.5) being a random classifier. The (kappa ) statistic is a measure of observer reliability, usually employed to test the agreement between two systems. When applied to confusion matrix evaluations, (kappa ) can be used to assess the probability that a model will produce an output (hat{y}) that coincides with the real output y. (kappa ) values typically range between (left[ 0, 1right] ), with (kappa =1) meaning perfect agreement, (kappa =0) being random agreement, and (kappa =0.8) typically used as a threshold to define a near-perfect or perfect algorithm.While in the authors’ opinion, AUC, Sensitivity and Specificity values are the most reliable evaluation metrics for studies of this type (Supp. Appendix 6), for ease of comparison with other papers or authors who choose to use other metrics, we have also included Precision, Recall and F-Measure values. Precision and Recall values play a similar role to sensitivity and specificity, with recall being equivalent to sensitivity, and precision being the calculation of the number of correct positive predictions made. Precision and Recall, however, differ from their counterparts in being more robust to imbalance in datasets. F-Measures are a combined evaluation of these two measures. For more details consult Supplementary Appendix 6.Loss metrics were reported using the Mean Squared Error (Eq. 2);$$begin{aligned} MSE = frac{1}{n} sum _{i = 1}^{n} left( y_{i} – hat{y}_{i} right) ^{2} end{aligned}$$
(2)
Loss values are interpreted considering values closest to 0 as an indicator of greater confidence when using the model to make new predictions.Final evaluation metrics were reported when using algorithms to classify only the original samples, without augmented data. Augmented data was, therefore, solely used for training and validation. Finally, so as to assess the impact data augmentation has on supervised learning algorithms, algorithms were also trained on the raw data. This was performed using 70% of the raw data for training, while the remaining 30% was used as a test set. More

Study site and target foraminiferal speciesThe Maldivian archipelago is a partially drowned carbonate platform within the central, equatorial Indian Ocean. It consists of two rows of north–south orientated atolls, which encompass an Inner Sea. The lowermost neritic carbonate unit sits upon volcanic bedrock and has been dated back to the Eocene19 with continuous drift deposition, within the Inner Sea, starting ~ 12.9 Ma at the establishment of the modern South Asian Monsoon (SAM)35,36. This seasonally reversing, major climatic system has an impact on both the regional precipitation patterns as well as physiochemical oceanographic properties (Fig. 1). The summer southwest SAM brings warm, wet conditions to the Indian subcontinent, as well as higher saline surface waters from the Arabian Sea into the Maldives region. In comparison, the winter northeast SAM results in cool, dry continental conditions and transports lower salinity water from the Bay of Bengal into the central, equatorial Indian Ocean. As a result, the Maldives seasonal salinity depth profiles can vary significantly, yet due to its tropical location the seasonal sea water temperatures are relatively stable.Three symbiont-bearing foraminiferal species are used in this study: Amphistegina lessonii, Globigerinoides ruber (white) and Trilobatus sacculifer (with sac-like final chamber):Amphistegina lessonii is a larger benthic, symbiont-bearing (diatoms) foraminiferal species. It has a shallow depth range (0–50 m)37,38,39 and is globally abundant in tropical coral reef, benthic foraminiferal shoal and general carbonate shelf settings40. Similarly to corals, amphisteginids have been shown to bleach under high temperatures/high irradiance levels with the new development of the Amphistegina Bleaching Index (ABI) as an indicator of photo-inhibitory stress in coral reef settings41,42. From ~ 30 °C this species starts showing signs of thermal stress, with bleaching and mortality reported for temperatures  > 31 °C11,12.Globigerinoides ruber (w) hosts dinoflagellate endosymbionts and is the most common planktonic foraminiferal species in tropical-subtropical waters13 state that while G. ruber (w) is generally considered one of the shallowest-dwelling species, its depth distribution does vary in relation to regional ecological conditions. It has a particular relation to the nutricline depth in less turbid, oligotrophic conditions43 which has been confirmed for the Maldives28. It is omnivorous, however in comparison to other omnivorous, symbiont-bearing species, it has demonstrated an elevated adaptation for consuming phytoplankton protein over zooplankton protein13. From culture experiments, it has a broad temperature (14–31 °C) and salinity (22–49 PSU) tolerance, and has been reported as the most tolerant species to low sea surface salinity (SSS)13. This species occurs year-round and has a fortnightly reproduction13.Trilobatus sacculifer is a planktonic foraminiferal species abundant in tropical-subtropical surface waters and as such is extensively used in paleo-reconstructions. It hosts dinoflagellate endosymbionts yet is omnivorous, feeding predominantly on calanoid copepods13. It is a euryhaline species, with a broad salinity (24–47 PSU) and temperature (14–32 °C) tolerance. Similarly to G. ruber (w), this species occurs year-round and has a monthly reproduction on a synodic lunar cycle13. While a shallow dwelling species, it is generally reported to live slightly deeper in the water column, in comparison to G. ruber (w)28,30,44.SamplingAll planktonic foraminiferal specimens (G. ruber (w) and T. sacculifer (w/s)) for the geochemical analysis (δ18Oc and Mg/Ca) originate from the International Ocean Discovery Program (IODP) Expedition 359, Site U1467 (4° 51.0274′ N, 73° 17.0223′ E) drilled in 2015 within the Inner Sea of the Maldivian archipelago at a water depth of 487 m19. The age model for these samples was adopted from a previous study45 which is based on the correlation of their long-term (0–1800 kyr) Site 359-U1467 C. mabahethi and G. ruber (w) δ18Oc records to the stacked reference curve of46. Recent surface sediment samples (mudline A and B: representing the sample from the sediment/water interface), as well as three samples from the peak of MIS9e (U1467C, 2H6, 0–1 cm; U1467C, 2H6, 15–16 cm; U1467C, 2H6, 18–19 cm) and MIS11c (U1467B, 3H2, 147–148 cm; U1467B, 3H3, 9–10 cm; U1467B, 3H3, 12–13 cm) were analysed19,28 (sample locations are shown on Fig. 3). The mudline is identified as Recent, likely representing the last few hundred years, based on the presence of Rose Bengal (1 g/L) stained ostracods and benthic foraminifera. The study by45 has verified that diagenetic influences, within this shallow, carbonate environment, are not a concern for foraminiferal geochemical compositions over the investigated time-interval (MIS1-11).Rose Bengal stained A. lessonii specimens were obtained from modern surface rubble samples collected by hand, at 10 m water depth, during the 2015 International Union for Conservation of Nature (IUCN) REGENERATE cruise47 (Supplementary Table 6). Samples were collected from the reefs of two islands, Maayafushi and Rasdhoo, both located within the central part of the Maldivian archipelago. As the foraminifera shells were stained pink, it implies they were living at the time of collection. These specimens were used for stable isotopic analysis and their reconstructed temperatures represent modern (a cumulative signal encompassing their lifespan of four to twelve months48) conditions (Supplementary Tables 5–6). A full explanation of the Rose Bengal protein stain for foraminifera is detailed in49.δ18Oc stable isotopic analysisAll samples were initially washed using a 32 μm sieve to remove the finer clay and silt fractions. Subsequently, they were air dried and sieved into discrete sizes for foraminiferal picking. To ensure enough calcite for the measurements, all specimens for Individual Foraminifera Analysis (IFA) for both G. ruber (w) and T. sacculifer (w/s) (n = 632) were picked from the 355–400 μm size fraction. In addition, traditional whole-shell (pooled) measurements for G. ruber (w) (n = 24) were conducted on specimens from the 212–400 μm fraction (2–5 pooled specimens). Trilobatus sacculifer (w/s) traditional whole-shell analysis (n = 21) was measured on specimens (2 pooled specimens) from the 300–355 μm fraction. The majority of these pooled measurements are obtained from28,45,50,51 (Supplementary Table 1). Amphistegina lessonii measurements were run on single specimens  > 250 μm in size. Prior to stable isotopic analysis, all shells were briefly cleaned (1–2 s) by ultrasonication in Milli-Q water to remove any adhering particles. All stable isotopic measurements were conducted at the School of GeoSciences at the University of Edinburgh on a Thermo Electron Delta + Advantage mass spectrometer integrated with a Kiel carbonate III automated extraction line. Samples were reacted with 100% phosphoric acid (H3PO4) at 90 °C for 15 min, with the evolved CO2 gas collected in a liquid nitrogen coldfinger and analysed compared to a reference gas. All samples are corrected using an internal laboratory standard and expressed as parts per mil (‰) relative to Vienna Pee Dee Belemnite (VPDB). Replicate measurements of the standards give the instrument an analytical precision (1σ) of ~ 0.05 ‰ for δ18O and δ13C.Mg/Ca analysisThe Mg/Ca data is obtained from28,45,50,51 (Supplementary Table 1). Each G. ruber (w) Mg/Ca analysis (n = 17; 212–250 μm in size) was conducted on 30 pooled specimens by inductively coupled plasma optical emission spectrometry (ICP-OES) on a Thermoscientific iCap 6300 (dual viewing) at the Institute of Geosciences of the Goethe-University of Frankfurt. All samples were initially cleaned (1–2 s) by ultrasonication in Milli-Q water and then the standard oxidative cleaning protocol of52 followed to prevent clay mineral contamination. The final centrifuged sample solution was diluted with an yttrium solution (1 mg/l) prior to measurement to allow for the correction of matrix effects. In addition, before each analysis five calibration solutions were measured to allow for intensity ratio calibrations. All element/Ca measurements were standardized using an internal consistency standard (ECRM 752–1, 3.761 mmol/mol Mg/Ca). Furthermore, the elements Al, Fe, and Mn were screened and blanks periodically run to monitor for further signs of contamination during the analyses.Establishment of present and past seawater temperaturesPrior to temperature calculations, we test the IFA distributions for normality using the Shapiro‐Wilk test and the Fisher–Pearson coefficient of skewness with bootstrap confidence intervals, to define the skewness of the datasets53 (Supplementary Table 3). The Recent G. ruber (w) and T. sacculifer (w/s) and MIS11c T. sacculifer are normally distributed. In the case of both MIS9e datasets and the MIS11c G. ruber population, the null hypothesis that the data are normally distributed (p ≤ 0.05) is rejected (Supplementary Table 3). Considering bioturbation within the sediment record is a possibility, we use two methods to identify and remove outliers in the IFA datasets. Firstly, the inter-quartile range (IQR) is used for each δ18Oc dataset, which defines a measurement as an outlier if it falls outside the range [Q1 − 1.5 (Q3 − Q1), Q3 + 1.5 (Q3 − Q1)], with IQR = Q3 − Q1 and Q3 and Q1 representing the third and first quartile of the dataset20. But if there is considerable reworking, the IQR method would not necessarily identify reworked glacial measurements (highest δ18Oc values) within the interglacial samples. As such, the Recent IFA datasets, which are both normally distributed, are used to further set a rudimentary cut-off point for the highest δ18Oc (= lowest temperatures) value to expect during past interglacial minima periods for both G. ruber (w) and T. sacculifer (w/s) (this is discussed further in the Supplementary Materials, Supplementary Figs. 1–3).There are innumerable analytical techniques (e.g., traditional mass spectrometry, secondary-ion mass spectrometry, laser ablation inductively coupled plasma mass spectrometry), proxies (Mg/Ca, δ18O, clumped isotopes, TEX86, Uk’37) as well as target medians (e.g., calcitic shells of foraminifera, aragonitic coral skeletons, ice, lipids, alkenones) which are used in marine paleo-temperature reconstructions. Furthermore, different methods exist in the literature to calculate temperature estimates using both planktonic foraminiferal δ18Oc and Mg/Ca measurements with innumerable species-specific δ18O-temperature and Mg/Ca-temperature equations reported20,23,30,54,55,56. Moreover, due to the exponential nature of the Mg/Ca-temperature equations, if inappropriately applied, offsets in the upper temperature range are exacerbated. Additional considerations are species-specific offsets and differential geochemical compositions within the shell (e.g., high versus low Mg banding, gametogenic calcite). Trilobatus sacculifer gametogenic calcite has been reported to be significantly enriched in Mg in comparison to the rest of the shell57. As T. sacculifer specimens selected for use in this study underwent reproduction, indicated by the presence of a sac-like final chamber58, we can expect their Mg/Ca ratios to be biased. As such, to avoid overestimates we chose to use only G ruber (w, pooled) Mg/Ca and δ18Oc data to calculate representative δ18Osw values for each time interval, for use with both the G. ruber (w) and T. sacculifer (w/s) δ18Oc IFA datasets. Considering both planktonic species are considered as shallow-dwellers with similar living depths and an affinity for the DCM, the utilisation of common δ18Osw values is applicable13,28,30.The G. ruber Mg/Ca-temperature Eq. (1) from55 (temperature calibration range: ~ 22–27 °C), similarly applied in the Maldivian study of28, was used in this study:$$Mg/Ca=0.34left(pm 0.08right)mathrm{exp}(0.102left(pm 0.010right)*T)$$
(1)
The applied δ18O-temperature species-specific equations (Eqs. 2 and 3) were previously utilised in the local study by28. Both the G. ruber (Eq. 2) and T. sacculifer (Eq. 3) equations are from the Indian Ocean study of59 (temperature calibration range: ~ 20–31 °C):$$T=12.75-5({delta }^{18}{O}_{c}-{delta }^{18}{O}_{sw})$$
(2)
$$T=11.95-5.26({delta }^{18}{O}_{c}-{delta }^{18}{O}_{sw})$$
(3)
Using the above equations, the range in temperature estimates are obtained as follows (Fig. 4):

1.

The mean G. ruber (w) Mg/Ca measurements are used together with Eq. (1) to calculate a temperature estimate for each time point (Supplementary Table 1). Since the Mg/Ca calcification temperatures are based on 30 pooled specimens, they are considered to reflect mean calcification temperatures.

2.

The Mg/Ca derived temperature estimates are then used together with the mean traditional (pooled) G. ruber (w) δ18Oc data and Eq. (2) to calculate representative δ18Osw values for each time point (Supplementary Table 2). As these are calculated from pooled samples, they are considered to mirror mean δ18Osw values for both the Recent and fossil populations.

3.

The G. ruber (w) and T. sacculifer (w/s) IFA datasets are then used, together with the relevant species-specific δ18O-temperature equations and δ18Osw values, to calculate the spread in temperature estimates (Fig. 4, Supplementary Tables 3–4).

Trilobatus sacculifer (w/s) data from the glacial maxima of MIS12 are included in the study to illustrate the applicability of the IFA method, however, as they do not contribute to the discussion on bleaching thresholds, they are discussed further in the Supplementary Materials (Supplementary Figs. 1, 3).Finally, the temperature estimates for the shallow-dwelling symbiont-bearing benthic A. lessonii are obtained using the genus-specific δ18O-temperature equation of60 (Eq. 4) (Supplementary Tables 5–6).$$T=16.3-4.24({delta }^{18}{O}_{c}-{delta }^{18}{O}_{sw})$$
(4)
Considering the benthic specimens were deemed living at the time of collection (Rose Bengal stained), a mean regional surface (0 m) δ18Osw value (0.49 ‰) is used together with the δ18Oc data in the calculations (Supplementary Tables 5–6). More

1.EU. No 1143/2014 of the European Parliament and of the Council of 22 October 2014 on the prevention and management of the introduction and spread of invasive alien species. Off. J. Eur. Union 2014, 35–55 (2014).
Google Scholar 
2.Ludyanskiy, M. L., McDonald, D. & MacNeill, D. Impact of the Zebra Mussel, a Bivalve Invader. Bioscience 43, 533–544 (1993).Article 

Google Scholar 
3.Lalaguna, C. D. & Marco, A. A. The zebra mussel invasion in Spain and navigation rules. Aquat. Invasions 3, 315–324 (2008).Article 

Google Scholar 
4.Rajagopal, S. et al. Origin of spanish invasion by the zebra mussel, dreissena polymorpha (pallas, 1771) revealed by amplified fragment length polymorphism (AFLP) fingerprinting. Biol. Invasions 11, 2147–2159 (2009).Article 

Google Scholar 
5.CABI Invasive species Compendium: Dreissena polymorpha (zebra mussel) 2017. https://www.cabi.org/isc/datasheet/85295.6.Marescaux, J. & Van Doninck, K. Using DNA barcoding to differentiate invasive Dreissena species (Mollusca, Bivalvia). Zookeys 365, 235–244 (2013).Article 

Google Scholar 
7.Benson, A. J., Raikow, D., Larson, J., Fusaro, A., & Bogdanoff, A. K. Dreissena polymorpha (Pallas, 1771): U.S. Geological Survey, Nonindigenous Aquatic Species Database. Gainesville, FL (2017).8.Minchin, D., Lucy, F. & Sullivan, M. Zebra Mussel: Impacts and Spread. Invasive Aquat. Species Eur. Distrib. Impacts Manag. 135–146. © 2002 Kluwer Acad. Publ. Dordreicht, Netherlands. 135–148 (2002) doi:https://doi.org/10.1007/978-94-015-9956-6_15.9.Montero Melendez, J. Control of invasive alien species in Guadalquivir river basin. EURO-RIOB 2017. https://www.riob.org/es/node/404510.Molloy, D. P., Karatayev, A., Burlakova, L. E., Kurandina, D. P. & Laruelle, F. Natural enemies of zebra mussels: predators, parasites, and ecological competitors. Rev. Fish. Sci. 5, 27–97 (1997).Article 

Google Scholar 
11.Nalepa, T. F. & Schloesser, D. W. Zebra Mussels Biology, Impacts, and Control (Lewis Publishers, Boca Raton, 1993).
Google Scholar 
12.Birnbaum, C. NOBANIS – Invasive Alien Species Fact Sheet – Dreissena polymorpha. Accessed 2 Dec 2019. https://www.nobanis.org (Online Database of the European Network on Invasive Alien Species – NOBANIS, 2011).13.Lowe, S., Browne, M., Boudjelas, S., De Poorter, M. 100 of the World’s Worst Invasive Alien Species A selection from the Global Invasive Species
Database. (The Invasive Species Specialist Group (ISSG), 2000).14.Glomski, L. M. Zebra Mussel Chemical Control Guide. US Army Corps of Engineers: Waterways Experiment Station. https://erdclibrary.erdc.dren.mil/jspui/bitstream/11681/6966/1/ERDC-EL-TR-15-9.pdf (2015). 15.Boelman, S. F., Neilson, F. M., Dardeau, E. A. & Cross, T. Zebra mussel (Dreissena polymorpha) control handbook for facility operators, first edition . US Army Corps of Engineers: Waterways Experiment Station. https://hdl.handle.net/11681/2966 (1997).16.Durán, C., Lanao, M., Anadón, A. & Touyá, V. Management in practice management strategies for the zebra mussel invasion in the Ebro River basin. Aquat. Invasions 5, 309–316 (2010).Article 

Google Scholar 
17.Bij de Vaate, A. Rajagopal, S. & van der Velde, G. The zebra mussel in Europe: summary and synthesis. in The Zebra Mussel in Europe (ed van der Velde, G.
et al.) 415–421 (Backhuys Publishers, 2010).18.Herder, J. et al. Environmental DNA—a review of the possible applications for the detection of (invasive) species. Report 2013-104. Accessed 3 May 2019. https://www.researchgate.net/publication/283267157_Environmental_DNA_-_a_review_of_the_possible_applications_for_the_detection_of_invasive_species#fullTextFileContent (Stichting RAVON, Nijmegen, 2014).19.Xiong, W., Li, H. & Zhan, A. Early detection of invasive species in marine ecosystems using high-throughput sequencing: technical challenges and possible solutions. Mar. Biol. 163, 1–12 (2016).CAS 
Article 

Google Scholar 
20.Harvey, C. T., Qureshi, S. A. & MacIsaac, H. J. Detection of a colonizing, aquatic, non-indigenous species. Divers. Distrib. 15, 429–437 (2009).Article 

Google Scholar 
21.Jerde, C. L., Mahon, A. R., Chadderton, W. L. & Lodge, D. M. ‘Sight-unseen’ detection of rare aquatic species using environmental DNA. Conserv. Lett. 4, 150–157 (2011).Article 

Google Scholar 
22.Dejean, T. et al. Improved detection of an alien invasive species through environmental DNA barcoding: the example of the American bullfrog Lithobates catesbeianus. J. Appl. Ecol. 49, 953–959 (2012).Article 

Google Scholar 
23.Jerde, C. L. et al. Detection of Asian carp DNA as part of a Great Lakes basin-wide surveillance program. Can. J. Fish. Aquat. Sci. 70, 522–526 (2013).CAS 
Article 

Google Scholar 
24.Laramie, M. B., Pilliod, D. S. & Goldberg, C. S. Characterizing the distribution of an endangered salmonid using environmental DNA analysis. Biol. Conserv. 183, 29–37 (2015).Article 

Google Scholar 
25.Takahara, T., Minamoto, T. & Doi, H. Using environmental DNA to estimate the distribution of an invasive fish species in ponds. PLoS ONE 8, e56584 (2013).ADS 
CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
26.Gingera, T. D., Bajno, R., Docker, M. F. & Reist, J. D. Environmental DNA as a detection tool for zebra mussels Dreissena polymorpha (Pallas, 1771) at the forefront of an invasion event in Lake Winnipeg, Manitoba, Canada. Manag. Biol. Invasions 8, 287–300 (2017).Article 

Google Scholar 
27.Darling, J. A. & Mahon, A. R. From molecules to management: adopting DNA-based methods for monitoring biological invasions in aquatic environments. Environ. Res. 111, 978–988 (2011).CAS 
PubMed 
Article 

Google Scholar 
28.Kaprou, G. D. et al. Miniaturized devices for isothermal DNA amplification addressing DNA diagnostics. Microsyst. Technol. 22, 1529–1534 (2016).CAS 
Article 

Google Scholar 
29.Mori, Y. & Notomi, T. Loop-mediated isothermal amplification (LAMP): a rapid, accurate, and cost-effective diagnostic method for infectious diseases. J. Infect. Chemother. 15, 62–69 (2009).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
30.Fang, X., Liu, Y., Kong, J. & Jiang, X. Loop-mediated isothermal amplification integrated on microfluidic chips for point-of-care quantitative detection of pathogens. Anal. Chem. 82, 3002–3006 (2010).CAS 
PubMed 
Article 

Google Scholar 
31.Rafati, A. & Gill, P. Microfluidic method for rapid turbidimetric detection of the DNA of Mycobacterium tuberculosis using loop-mediated isothermal amplification in capillary tubes. Microchim. Acta 182, 523–530 (2014).Article 
CAS 

Google Scholar 
32.Notomi, T. et al. Loop-mediated isothermal amplification of DNA. Nucleic Acids Res. 28, e63 (2000).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
33.Mori, Y., Kitao, M., Tomita, N. & Notomi, T. Real-time turbidimetry of LAMP reaction for quantifying template DNA. J. Biochem. Biophys. Methods 59, 145–157 (2004).CAS 
PubMed 
Article 

Google Scholar 
34.Tomita, N., Mori, Y., Kanda, H. & Notomi, T. Loop-mediated isothermal amplification (LAMP) of gene sequences and simple visual detection of products. Nat. Protoc. 3, 877–882 (2008).CAS 
PubMed 
Article 

Google Scholar 
35.Garrido-Maestu, A., Fuciños, P., Azinheiro, S., Carvalho, J. & Prado, M. Systematic loop-mediated isothermal amplification assays for rapid detection and characterization of Salmonella spp., Enteritidis and Typhimurium in food samples. Food Control 80, 297–306 (2017).CAS 
Article 

Google Scholar 
36.Verkaar, E., Nijman, I., Boutaga, K. & Lenstra, J. Differentiation of cattle species in beef by PCR-RFLP of mitochondrial and satellite DNA. Meat Sci. 60, 365–369 (2002).CAS 
PubMed 
Article 

Google Scholar 
37.Wilson-Wilde, L., Norman, J. & Robertson, J. et al. Current issues in species identification for forensic science and the validity of using the cytochrome oxidase I
(COI) gene. Forensic Sci. Med. Pathol. 6, 233–241. https://doi.org/10.1007/s12024-010-9172-y (2010). CAS 
PubMed 
Article 

Google Scholar 
38.Hebert, P. D. N., Cywinska, A., Ball, S. L. & Jeremy, R. Biological identifications through DNA barcodes. Proc. Biol. Sci. 270(1512), 313–321 (2003).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
39.Staroscik, A. Copy number calculator for realtime PCR. http://www.scienceprimer.com/copy-number-calculator-for-realtime-pcr (2012).40.Ficetola, G. F., Miaud, C., Pompanon, F. & Taberlet, P. Species detection using environmental DNA from water samples. Biol. Lett. 4, 423–425 (2008).PubMed 
PubMed Central 
Article 

Google Scholar 
41.Zanoli, L. M. & Spoto, G. Isothermal amplification methods for the detection of nucleic acids in microfluidic devices. Biosensors 3, 18–43 (2013).CAS 
PubMed 
Article 

Google Scholar 
42.Egan, S. P. et al. Rapid molecular detection of invasive species in ballast and harbor water by integrating environmental DNA and light transmission spectroscopy. Environ. Sci. Technol. 49, 4113–4121 (2015).ADS 
CAS 
PubMed 
Article 

Google Scholar 
43.Xia, Z. et al. Early detection of a highly invasive bivalve based on environmental DNA (eDNA). Biol. Invasions 20, 437–447 (2018).Article 

Google Scholar 
44.Williams, M. R. et al. Isothermal amplification of environmental DNA (eDNA) for direct field-based monitoring and laboratory confirmation of Dreissena sp. PLoS ONE 12, e0186462 (2017).PubMed 
PubMed Central 
Article 
CAS 

Google Scholar 
45.Frackman, B. S., Kobs, G., Simpson, D., Storts, D. & Corporation, P. Betaine and DMSO: enhancing agents for PCR. Promega Notes 65, 9–12 (1998).
Google Scholar 
46.Wang, D.-G., Brewster, J., Paul, M. & Tomasula, P. Two methods for increased specificity and sensitivity in loop-mediated isothermal amplification. Molecules 20, 6048–6059 (2015).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
47.Jantz, B. & Neumann, D. Growth and reproductive cycle of the zebra mussel in the River Rhine as studied in a river bypass. Oecologia 114, 213–225 (1998).ADS 
CAS 
PubMed 
Article 

Google Scholar 
48.Grigorovich, I. A., Kelly, J. R., Darling, J. A. & West, C. W. The quagga mussel invades the Lake Superior Basin. J. Great Lakes Res. 34, 342–350 (2008).Article 

Google Scholar 
49.Mahon, A. R. et al. Molecular detection of invasive species in heterogeneous mixtures using a microfluidic carbon nanotube platform. PLoS ONE 6, 1–5 (2011).Article 
CAS 

Google Scholar 
50.PrimerExplorer. LAMP Primer Designing Software (Fujitsu Ltd, Tokyo, 2005).
Google Scholar 
51.Untergasser, A. et al. Primer3-new capabilities and interfaces. Nucleic Acids Res. 40, e115–e115 (2012).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar 
52.Altschul, S. F., Gish, W., Miller, W., Myers, E. W. & Lipman, D. J. Basic local alignment search tool. J. Mol. Biol. 215, 403–410 (1990).CAS 
PubMed 
PubMed Central 
Article 

Google Scholar  More

1.Smol, J. P. et al. Climate-driven regime shifts in the biological communities of arctic lakes. Proc. Natl Acad. Sci. USA 102, 4397–4402 (2005).CAS 
Article 

Google Scholar 
2.Wernberg, T. et al. Climate-driven regime shift of a temperate marine ecosystem. Science 353, 169–172 (2016).CAS 
Article 

Google Scholar 
3.Scheffer, M., Carpenter, S., Foley, J. A., Folke, C. & Walker, B. Catastrophic shifts in ecosystems. Nature 413, 591–596 (2001).CAS 
Article 

Google Scholar 
4.Staver, A. C., Archibald, S. & Levin, S. A. The global extent and determinants of savanna and forest as alternative biome states. Science 334, 230–232 (2011).CAS 
Article 

Google Scholar 
5.Su, H. et al. Long‐term empirical evidence, early warning signals and multiple drivers of regime shifts in a lake ecosystem. J. Ecol. https://doi.org/10.1111/1365-2745.13544 (2020).6.Barnosky, A. D. et al. Approaching a state shift in Earth’s biosphere. Nature 486, 52–58 (2012).CAS 
Article 

Google Scholar 
7.Steffen, W. et al. Trajectories of the Earth system in the Anthropocene. Proc. Natl Acad. Sci. 115, 8252–8259 (2018).CAS 
Article 

Google Scholar 
8.Holling, C. S. Resilience and stability of ecological systems. Ann. Rev. Ecol. Syst. 4, 1–23 (1973).Article 

Google Scholar 
9.Ratajczak, Z. et al. Abrupt change in ecological systems: inference and diagnosis. Trends Ecol. Evol. 33, 513–526 (2018).Article 

Google Scholar 
10.Pimm, S. L. The complexity and stability of ecosystems. Nature 307, 321–326 (1984).Article 

Google Scholar 
11.Holling, C. S. Engineering resilience versus ecological resilience. Eng. Ecol.Constraints 31, 32 (1996).
Google Scholar 
12.Li, X. et al. Temporal trade-off between gymnosperm resistance and resilience increases forest sensitivity to extreme drought. Nat. Ecol. Evol. 4, 1075–1083 (2020).Article 

Google Scholar 
13.Carpenter, S. R. & Brock, W. A. Rising variance: a leading indicator of ecological transition. Ecol. Lett. 9, 311–318 (2006).CAS 
Article 

Google Scholar 
14.Dakos, V. et al. Slowing down as an early warning signal for abrupt climate change. Proc. Natl Acad. Sci. USA 105, 14308–14312 (2008).CAS 
Article 

Google Scholar 
15.Guttal, V. & Jayaprakash, C. Changing skewness: an early warning signal of regime shifts in ecosystems. Ecol. Lett. 11, 450–460 (2008).Article 

Google Scholar 
16.Scheffer, M. et al. Early-warning signals for critical transitions. Nature 461, 53–59 (2009).CAS 
Article 

Google Scholar 
17.Drake, J. M. & Griffen, B. D. Early warning signals of extinction in deteriorating environments. Nature 467, 456 (2010).CAS 
Article 

Google Scholar 
18.Wang, R. et al. Flickering gives early warning signals of a critical transition to a eutrophic lake state. Nature 492, 419–422 (2012).Article 
CAS 

Google Scholar 
19.Clements, C. F. & Ozgul, A. Including trait-based early warning signals helps predict population collapse. Nat. Commun. 7, 10984 (2016).CAS 
Article 

Google Scholar 
20.Chevalier, M. & Grenouillet, G. Global assessment of early warning signs that temperature could undergo regime shifts. Sci. Rep. 8, 10058 (2018).Article 
CAS 

Google Scholar 
21.Cole, L. E., Bhagwat, S. A. & Willis, K. J. Recovery and resilience of tropical forests after disturbance. Nat. Commun. 5, 3906 (2014).CAS 
Article 

Google Scholar 
22.Willis, K. J., Jeffers, E. S. & Tovar, C. What makes a terrestrial ecosystem resilient? Science 359, 988–989 (2018).CAS 
Article 

Google Scholar 
23.Thomas, C. D. et al. Extinction risk from climate change. Nature 427, 145–148 (2004).CAS 
Article 

Google Scholar 
24.Seddon, A. W., Macias-Fauria, M., Long, P. R., Benz, D. & Willis, K. J. Sensitivity of global terrestrial ecosystems to climate variability. Nature 531, 229–232 (2016).CAS 
Article 

Google Scholar 
25.Ehleringer, J. R., Cerling, T. E. & Helliker, B. R. C4 photosynthesis, atmospheric CO2, and climate. Oecologia 112, 285–299 (1997).Article 

Google Scholar 
26.Higgins, S. I. & Scheiter, S. Atmospheric CO2 forces abrupt vegetation shifts locally, but not globally. Nature 488, 209 (2012).CAS 
Article 

Google Scholar 
27.Holmgren, M., Hirota, M., Van Nes, E. H. & Scheffer, M. Effects of interannual climate variability on tropical tree cover. Nat. Clim. Chang. 3, 755–758 (2013).Article 

Google Scholar 
28.Thornton, P. K., Ericksen, P. J., Herrero, M. & Challinor, A. J. Climate variability and vulnerability to climate change: a review. Glob. Change Biol. 20, 3313–3328 (2014).Article 

Google Scholar 
29.Ray, D. K., Gerber, J. S., MacDonald, G. K. & West, P. C. Climate variation explains a third of global crop yield variability. Nat. Commun. 6, 5989 (2015).CAS 
Article 

Google Scholar 
30.Jha, S., Das, J. & Goyal, M. K. Assessment of risk and resilience of terrestrial ecosystem productivity under the influence of extreme climatic conditions over India. Sci. Rep. 9, 18923 (2019).CAS 
Article 

Google Scholar 
31.Li, D., Wu, S., Liu, L., Zhang, Y. & Li, S. Vulnerability of the global terrestrial ecosystems to climate change. Glob. Change Biol. 24, 4095–4106 (2018).Article 

Google Scholar 
32.Gonzalez, P., Neilson, R. P., Lenihan, J. M. & Drapek, R. J. Global patterns in the vulnerability of ecosystems to vegetation shifts due to climate change. Glob. Ecol. Biogeogr. 19, 755–768 (2010).Article 

Google Scholar 
33.Wang, S. & Loreau, M. Ecosystem stability in space: α, β and γ variability. Ecol. Lett. 17, 891–901 (2014).Article 

Google Scholar 
34.Stenseth, N. C. et al. The effect of climatic forcing on population synchrony and genetic structuring of the Canadian lynx. Proc. Natl Acad. Sci. USA 101, 6056–6061 (2004).CAS 
Article 

Google Scholar 
35.Koenig, W. D. & Liebhold, A. M. Temporally increasing spatial synchrony of North American temperature and bird populations. Nat. Clim. Chang. 6, 614–617 (2016).Article 

Google Scholar 
36.Sheppard, L. W., Bell, J. R., Harrington, R. & Reuman, D. C. Changes in large-scale climate alter spatial synchrony of aphid pests. Nat. Clim. Chang. 6, 610–613 (2016).Article 

Google Scholar 
37.Dakos, V., van Nes, E. H., Donangelo, R., Fort, H. & Scheffer, M. Spatial correlation as leading indicator of catastrophic shifts. Theor. Ecol. 3, 163–174 (2010).Article 

Google Scholar 
38.Paruelo, J. M., Epstein, H. E., Lauenroth, W. K. & Burke, I. C. ANPP estimates from NDVI for the central grassland region of the United States. Ecology 78, 953–958 (1997).Article 

Google Scholar 
39.Piao, S., Fang, J., Zhou, L., Tan, K. & Tao, S. Changes in biomass carbon stocks in China’s grasslands between 1982 and 1999. Global Biogeochem. Cycles 21, 2 (2007).
Google Scholar 
40.Maurer, G. E., Hallmark, A. J., Brown, R. F., Sala, O. E. & Collins, S. L. Sensitivity of primary production to precipitation across the United States. Ecol. Lett. 23, 527–536 (2020).Article 

Google Scholar 
41.Brown, J. H. & Kodric-Brown, A. Turnover rates in insular biogeography: effect of immigration on extinction. Ecology 58, 445–449 (1977).Article 

Google Scholar 
42.Earn, D. J., Levin, S. A. & Rohani, P. Coherence and conservation. Science 290, 1360–1364 (2000).CAS 
Article 

Google Scholar 
43.Hodgson, D., McDonald, J. L. & Hosken, D. J. What do you mean,‘resilient’? Trends Ecol. Evol. 30, 503–506 (2015).Article 

Google Scholar 
44.Seidl, R. et al. Forest disturbances under climate change. Nat. Clim. Chang. 7, 395–402 (2017).Article 

Google Scholar 
45.Bernstein, L. et al. IPCC, 2007: Climate Change 2007: Synthesis Report. (IPCC, Geneva, 2008)46.Myers-Smith, I. H. et al. Shrub expansion in tundra ecosystems: dynamics, impacts and research priorities. Environ. Res. Lett. 6, 045509 (2011).Article 

Google Scholar 
47.Myers-Smith, I. H. et al. Climate sensitivity of shrub growth across the tundra biome. Nat. Clim. Chang. 5, 887–891 (2015).Article 

Google Scholar 
48.Thompson, I., Mackey, B., McNulty, S. & Mosseler, A. Forest resilience, biodiversity, and climate change. In Secretariat of the Convention on Biological Diversity, Montreal. Technical Series 43, 1–67 (2009).
Google Scholar 
49.Carpenter, S. R. et al. Early warnings of regime shifts: a whole-ecosystem experiment. Science 332, 1079–1082 (2011).CAS 
Article 

Google Scholar 
50.Gsell, A. S. et al. Evaluating early-warning indicators of critical transitions in natural aquatic ecosystems. Proc. Natl Acad. Sci. USA 113, E8089–E8095 (2016).CAS 
Article 

Google Scholar 
51.Clements, C. F., Blanchard, J. L., Nash, K. L., Hindell, M. A. & Ozgul, A. Body size shifts and early warning signals precede the historic collapse of whale stocks. Nat. Ecol. Evol. 1, 0188 (2017).Article 

Google Scholar 
52.Dakos, V., Carpenter, S. R., van Nes, E. H. & Scheffer, M. Resilience indicators: prospects and limitations for early warnings of regime shifts. Philos. Trans. R. Soc. B, Biol. Sci. 370, 20130263 (2015).Article 

Google Scholar 
53.Zemp, D. C. et al. Self-amplified Amazon forest loss due to vegetation-atmosphere feedbacks. Nat. Commun. 8, 14681 (2017).CAS 
Article 

Google Scholar 
54.Staal, A. et al. Forest-rainfall cascades buffer against drought across the Amazon. Nat. Clim. Chang. 8, 539–543 (2018).Article 

Google Scholar 
55.Poorter, L. et al. Biomass resilience of Neotropical secondary forests. Nature 530, 211–214 (2016).CAS 
Article 

Google Scholar 
56.Locosselli, G. M. et al. Global tree-ring analysis reveals rapid decrease in tropical tree longevity with temperature. Proc. Natl Acad. Sci. USA 117, 33358–33364 (2020).CAS 
Article 

Google Scholar 
57.Ruiz-Pérez, G. & Vico, G. Effects of temperature and water availability on Northern European boreal forests. Front. For. Glob.Change 3, 34 (2020).Article 

Google Scholar 
58.Kitzberger, T., Aráoz, E., Gowda, J. H., Mermoz, M. & Morales, J. M. Decreases in fire spread probability with forest age promotes alternative community states, reduced resilience to climate variability and large fire regime shifts. Ecosystems 15, 97–112 (2012).Article 

Google Scholar 
59.Scheffer, M., Hirota, M., Holmgren, M., Van Nes, E. H. & Chapin, F. S. Thresholds for boreal biome transitions. Proc. Natl Acad. Sci. USA 109, 21384–21389 (2012).CAS 
Article 

Google Scholar 
60.Newbold, T. et al. Climate and land-use change homogenise terrestrial biodiversity, with consequences for ecosystem functioning and human well-being. Emerg. Top. Life Sci. 3, 207–219 (2019).Article 

Google Scholar 
61.Senior, R. A., Hill, J. K., González del Pliego, P., Goode, L. K. & Edwards, D. P. A pantropical analysis of the impacts of forest degradation and conversion on local temperature. Ecol. Evol. 7, 7897–7908 (2017).Article 

Google Scholar 
62.Wang, S. et al. An invariability-area relationship sheds new light on the spatial scaling of ecological stability. Nat. Commun. 8, 1–8 (2017).Article 
CAS 

Google Scholar 
63.Mehrabi, Z. & Ramankutty, N. Synchronized failure of global crop production. Nat. Ecol. Evol. 3, 780–786 (2019).Article 

Google Scholar 
64.Post, E. & Forchhammer, M. C. Spatial synchrony of local populations has increased in association with the recent Northern Hemisphere climate trend. Proc. Natl Acad. Sci. 101, 9286–9290 (2004).CAS 
Article 

Google Scholar 
65.Ripa, J. Analysing the Moran effect and dispersal: their significance and interaction in synchronous population dynamics. Oikos 89, 175–187 (2000).Article 

Google Scholar 
66.Peterson, G., Allen, C. R. & Holling, C. S. Ecological resilience, biodiversity, and scale. Ecosystems 1, 6–18 (1998).Article 

Google Scholar 
67.Wang, S. & Loreau, M. Biodiversity and ecosystem stability across scales in metacommunities. Ecol. Lett. 19, 510–518 (2016).Article 

Google Scholar 
68.Dakos, V. et al. Methods for detecting early warnings of critical transitions in time series illustrated using simulated ecological data. PloS ONE 7, e41010 (2012).CAS 
Article 

Google Scholar 
69.R core team. R: a language and environment for statistical computing. R Foundation for Statistical Computing https://www.R-project.org/ (2019).70.Bivand, R., Keitt, T. & Rowlingson, B. rgdal: bindings for the ‘Geospatial’ Data Abstraction Library. R package version 1.5-16 https://CRAN.R-project.org/package=rgdal (2020).71.Tucker, C. J. et al. An extended AVHRR 8‐km NDVI dataset compatible with MODIS and SPOT vegetation NDVI data. Int. J. Remote Sens. 26, 4485–4498 (2005).Article 

Google Scholar 
72.Pinzon, J. E. & Tucker, C. J. A non-stationary 1981–2012 AVHRR NDVI3g time series. Remote Sens. 6, 6929–6960 (2014).Article 

Google Scholar 
73.Holben, B. N. Characteristics of maximum-value composite images from temporal AVHRR data. Int. J. Remote Sens. 7, 1417–1434 (1986).Article 

Google Scholar 
74.Piao, S. et al. Changes in vegetation net primary productivity from 1982 to 1999 in China. Global Biogeochem. Cycles 19, 2 (2005).Article 
CAS 

Google Scholar 
75.Olson, D. M. et al. Terrestrial ecoregions of the world: a new map of life on Earth. BioScience 51, 933–938 (2001).Article 

Google Scholar 
76.Harris, I., Osborn, T. J., Jones, P. & Lister, D. Version 4 of the CRU TS monthly high-resolution gridded multivariate climate dataset. Sci. Data 7, 1–18. (2020).Article 

Google Scholar 
77.Mitchell, A. The ESRI Guide to GIS Analysis: Spatial Measurements and Statistics (Environmental System Research Institute Press, 2005).78.Fang, J., Piao, S., He, J. & Ma, W. Increasing terrestrial vegetation activity in China, 1982–1999. Sci. China C Life Sci. 47, 229–240 (2004).
Google Scholar 
79.Peng, S. et al. Recent change of vegetation growth trend in China. Environ. Res. Lett. 6, 044027 (2011).Article 

Google Scholar 
80.Thenkabail, P. S. & Lyon, J. G. Hyperspectral Remote Sensing of Vegetation (CRC press, 2016).81.Feng, Y. et al. Changes in the trends of vegetation net primary productivity in China between 1982 and 2015. Environ. Res. Lett. 14, 124009 (2019).Article 

Google Scholar 
82.He, H. et al. Altered trends in carbon uptake in China’s terrestrial ecosystems under the enhanced summer monsoon and warming hiatus. Natl Sci. Rev. 6, 505–514 (2019).CAS 
Article 

Google Scholar  More