Philippa Kaur

1.
Strong, D. R., Lawton, J. H. & Southwood, S. R. Insects on Plants, Community Patterns and Mechanisms (Blackwell Scientific Publicatons, Hoboken, 1984).
Google Scholar 
2.
Daly, H. V., Doyen, J. T. & Purcell, A. H. Introduction to Insect Biology and Diversity (Oxford University Press, Oxford, 1998).
Google Scholar 

3.
Strong, D. R., Lawton, J. H. & Southwood, T. R. E. Insects on Plants: Community Patterns and Mechanisms (Harvard University Press, Cambridge, 1984).
Google Scholar 

4.
Mitter, C., Farrell, B. D. & Wiegmann, B. The phylogenetic study of adaptive zones: Has phytophagy promoted insect diversification?. Am. Nat 132, 107–128 (1988).
Article  Google Scholar 

5.
Farrell, B. D. “Inordinate fondness” explained: Why are there so many beetles?. Science 281, 555–559 (1998).
CAS  PubMed  Article  Google Scholar 

6.
Marvaldi, A. E., Sequeira, A. S., O’Brien, C. W. & Farrell, B. D. Molecular and morphological phylogenetics of weevils (Coleoptera, Curculionoidea): Do niche shifts accompany diversification?. Syst. Biol 51, 761–785 (2002).
PubMed  Article  Google Scholar 

7.
Futuyma, D. J. Evolution of Host Specificity in Herbivorous Insects: Genetic, Ecological, and Phylogenetic Aspects (Wiley, New York, 1991).
Google Scholar 

8.
Via, S., Bouck, A. C. & Skillman, S. Reproductive isolation between divergent races of pea aphids on two hosts. II. Selection against migrants and hybrids in the parental environments. Evolution 54, 1626–1637 (2000).
CAS  PubMed  Article  Google Scholar 

9.
Nosil, P. Transition rates between specialization and generalization in phytophagous insects. Evolution 56, 1701–1706 (2002).
CAS  PubMed  Article  Google Scholar 

10.
Woodward, F. I. Climate and Plant Distribution (Cambridge University Press, Cambridge, 1987).
Google Scholar 

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

12.
Soberón, J. & Nakamura, M. Niches and distributional areas: Concepts, methods, and assumptions. Proc. Natl. Acad. Sci 106, 19644–19650 (2009).
ADS  PubMed  Article  Google Scholar 

13.
Whittaker, R. H. Evolution and measurement of species diversity. Taxon 21, 213–251 (1972).
Article  Google Scholar 

14.
Vellend, M. Conceptual synthesis in community ecology. Q Rev. Biol 85, 183–206 (2010).
PubMed  Article  Google Scholar 

15.
Novotny, V. Beta diversity of plant-insect food webs in tropical forests: A conceptual framework. Insect. Conserv. Diver 2, 5–9 (2009).
Article  Google Scholar 

16.
Legendre, P., Borcard, D. & Peres-Neto, P. R. Analyzing beta diversity: Partitioning the spatial variation of community composition data. Ecol. Monogr 75, 435–450 (2005).
Article  Google Scholar 

17.
Price, P. W. Resource-driven terrestrial interaction webs. Ecol. Res 17, 241–247 (2002).
Article  Google Scholar 

18.
Tscharntke, T. & Hawkins, B. A. Multi-Trophic Level Interactions (Cambridge University Press, Cambridge, 2002).
Google Scholar 

19.
Pereira, M. L., Matos, M. A., Lewinsohn, T. M. & Almeida, N. M. Trophic level and host specialisation affect beta-diversity in plant–herbivore–parasitoid assemblages. Insect. Conserv. Diver 12, 404–413 (2019).
Article  Google Scholar 

20.
Novotny, V. & Weiblen, G. D. From communities to continents: Beta diversity of herbivorous insects. Ann. Zool. Fenn 42, 463–475 (2005).
Google Scholar 

21.
Lewinsohn, T. M. & Roslin, T. Four ways towards tropical herbivore megadiversity. Ecol. Lett 11, 398–416 (2008).
PubMed  Article  Google Scholar 

22.
Ødegaard, F. Host specificity, alpha- and beta-diversity of phytophagous beetles in two tropical forests in Panama. Biodivers. Conserv 15, 83–105 (2006).
Article  Google Scholar 

23.
Hawkins, B. A. & Porter, E. E. Does herbivore diversity depend on plant diversity? The case of California butterflies. Am. Nat. 161, 40–49 (2003).
PubMed  Article  Google Scholar 

24.
Kemp, J. E., Linder, H. P. & Ellis, A. G. Beta diversity of herbivorous insects is coupled to high species and phylogenetic turnover of plant communities across short spatial scales in the Cape Floristic Region. J. Biogeogr. 44, 1813–1823 (2017).
Article  Google Scholar 

25.
Lawrence, J. F. Synopsis and Classification of Living Organisms, ***Vol 2 (McGraw-Hill, New York, 1982).
Google Scholar 

26.
Hanks, L. M. Influence of the larval host plant on reproductive strategies of cerambycid beetles. Annu. Rev. Entomol. 44, 483–505 (1999).
CAS  PubMed  Article  Google Scholar 

27.
Myers, N., Mittermeier, R. A., Mittermeier, C. G., Fonseca, G. A. B. & Kent, J. Biodiversity hotspots for conservation priorities. Nature 403, 853–857 (2000).
ADS  CAS  PubMed  Article  Google Scholar 

28.
Corff, J. & Marquis, R. Difference between understory and canopy in herbivore community composition and leaf quantity for two oak species in Missouri. Ecol. Entomol. 24, 46–58 (2001).
Article  Google Scholar 

29.
Ulyshen, M. D. & Hanula, J. L. A comparison of the beetle (Coleoptera) fauna captured at two heights above the ground in a North American temperate deciduous forest. Am. Midl. Nat. 158, 260–278 (2007).
Article  Google Scholar 

30.
Groot, P. & Nott, R. Evaluation of traps of six different designs to capture pine sawyer beetles (Coleoptera: Cerambycidae). Agric. Forest. Entomol 3, 107–111 (2002).
Article  Google Scholar 

31.
Miller, D. et al. Ipsenol, ipsdienol, ethanol, and pinene: Trap lure blend for Cerambycidae and Buprestidae (Coleoptera) in Pine Forests of Eastern North America. J. Econ. Entomol. 20, 20 (2015).
Google Scholar 

32.
Simpson, E. H. Measurement of diversity. Nature 163, 688 (1949).
ADS  MATH  Article  Google Scholar 

33.
Baselga, A. Partitioning the turnover and nestedness components of beta diversity. Global. Ecol. Biogeogr 19, 134–143 (2010).
Article  Google Scholar 

34.
Baselga, A. The relationship between species replacement, dissimilarity derived from nestedness, and nestedness. Glob. Ecol. Biogeogr 21, 1223–1232 (2012).
Article  Google Scholar 

35.
Carvalho, J. C. et al. Determining the relative roles of species replacement and species richness differences in generating beta-diversity patterns. Glob. Ecol. Biogeogr. 21, 760–771 (2012).
Article  Google Scholar 

36.
Vavrek, M. J. fossil: Palaeoecological and palaeogeographical analysis tools. Palaeontol. Electron 14, 1 (2011).
Google Scholar 

37.
Zhang, J.-L. Plantlist: Looking Up the Status of Plant Scientific Names based on The Plant List Database (Version 0.3.7) (2018).

38.
Webb, C. O. & Donoghue, M. J. Phylomatic: Tree assembly for applied phylogenetics. Mol. Ecol. Notes 5, 181–183 (2005).
Article  Google Scholar 

39.
Davis, T. J. et al. Darwin’s abominable mystery: Insights from a supertree of the angiosperms. Proc. Natl. Acad. Sci 101, 1904–1909 (2004).
ADS  Article  CAS  Google Scholar 

40.
Kembel, S. W. et al. Picante: R tools for integrating phylogenies and ecology. Bioinformatics 26, 1463–1464 (2010).
CAS  PubMed  Article  Google Scholar 

41.
Whittaker, R. H. Vegetation of the Siskiyou Mountains, Oregon and California. Ecol. Monogr 30, 279–338 (1960).
Article  Google Scholar 

42.
Burnham, K. P. & Anderson, D. R. Model Selection and Multimodel Inference (Springer, Berlin, 2002).
Google Scholar 

43.
R Core Team R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna. https://www.R-project.org (2018).

44.
Chesson, P. Mechanisms of maintenance of species diversity. Annu. Rev. Ecol. Evol. S 31, 343–366 (2000).
Article  Google Scholar 

45.
Whittaker, R. J. Communities and Ecosystems (MacMillan, New York, 1975).
Google Scholar 

46.
Araújo, M. B. & Luoto, M. The importance of biotic interactions for modelling species distributions under climate change. Glob. Ecol. Biogeogr 16, 743–753 (2007).
Article  Google Scholar 

47.
Elton, C. S. Animal Ecology (Sedgwick & Jackson Ltd, London, 1927).
Google Scholar 

48.
Grinnell, J. The niche-relationships of the California Thrasher. Auk 34, 427–433 (1917).
Article  Google Scholar  More

Pot experiment of AG planting
As shown in Fig. 1, compared to CS, the survival rate of 10-year rotation AG decreased, indicating that 2-year-old AG survival rate in RS was lower than that of AG in CS. This confirmed the continued existence of AG continuous cropping obstacles in RS.
The decrease of physicochemical properties and enzyme activity
Plant growth requires water and nutrients. Because soil physicochemical properties influence water and nutrient availability, changes in soil physicochemical properties directly affect AG growth. In the present study, the water content of RS was significantly higher than that of CS under the same management conditions (Table 1). Shu et al.29 found that high soil water content induced root rot disease in AG when sandy loam water content exceeded 30% or that of clay exceeded 50%. Similarly, according to Wang et al.30, the incidence of rust rot positively correlated with soil moisture and rainfall. Therefore, high soil water content, caused by changes in soil physicochemical properties, may negatively affect AG replanting. Furthermore, the pH of RS was significantly lower than that of CS (Table 1). According to Rahman and Punja24, root rot severity at soil pH 5.05 was greater than that at pH 7.0, indicating that acidic conditions can negatively affect AG health. In addition, the available K content in RS was lower than that in CS (Table 1). Sun31 found that AG should be fertilized (N, P, K fertilizer) from emergence to early flowering, when its demand for potassium fertilizer is the highest, suggesting that AG has a high potassium requirement. The levels of ammonium N, nitrate N, available P, and available K, but not of total N and total C, were generally lower in RS than in CS (Table 1), indicating that the cultivation of AG may have long-term negative effects on these soil nutrients. The same trend was observed for soil enzyme activity. Urease, a nickel-containing enzyme, catalyzes the hydrolysis of urea into carbonate and ammonia. Here, urease activity was significantly higher in CS than in RS. Average phosphatase and sucrase activities were also higher in CS than those in RS, although these differences were not significant (Table 2). Yang32 also found that the activities of sucrase, urease, and phosphatase decreased during AG cultivation. In summary, compared to that of CS, RS had lower fertility, but higher soil water content and lower pH, two conditions which are conducive to AG disease, and that may, therefore, present obstacles to AG replanting.
The dual effects of phenolic acids
The results showed that the content of salicylic acid in RS was significantly higher than that in CS. Yang16 found that among the various phenolic acids tested, salicylic acid had the strongest inhibitory effect on AG radicle growth. In our study, higher salicylic acid content in RS may have posed direct autotoxicity to AG. As a major defense hormone, salicylic acid has the function of enhancing immune signals and reprogramming defense transcriptomes33. After planting AG, the soil salicylic acid content increased, which indicated that AG might release more salicylic acid in the growth process to improve immune response to the surrounding environment. Therefore, the role of salicylic acid in the continuous cropping obstacles to AG cultivation deserves further study.
In addition, we found that the content of most phenolic acids, such as p-coumaric, p-hydroxybenzoic, vanillic, caffeic, and cinnamic acid, decreased after AG cultivation, and had not returned to the levels in CS even after 10 years of subsequent crop rotation. AG requires a suitable environment for growth. Before germination in spring, the ginseng farmers’ association uses wheat straw to cover the soil, which not only maintains soil temperature and retains soil moisture, but also improves soil quality and promotes the growth of AG seedlings. Jia et al.34 detected the increase in ferulic, vanillic, cinnamic, and p-hydroxybenzoic acid in a wheat-corn rotation area. In addition, Zheng et al.35 found that straw return, a common method for soil improvement, also increased the concentration of phenolic acids in soil. In our study, the increased phenolic acid content in CS relative to RS may have been beneficial to the growth of AG. Similar to our research results, Jiao et al.36 also found that the content of phenolic acid substances such as syringic, vanillic, p-coumaric, and ferulic acid decreased by 49.1–81% after adding AG root residues (simulating the seasonal AG leaf and fibrous root senescence). Therefore, decreases in the soil contents of some phenolic acids after planting AG may underlie the decline of other soil properties, which is not conducive to the subsequent growth of AG.
As described above, some phenolic acids may be beneficial to the growth of AG; if so, by what mechanism do these beneficial phenolic acids exert their role? Phenolic acids are produced by plants under external stress37,38,39,40. They do have many beneficial functions, such as antibacterial, antioxidant and so on, which can alleviate the stress of plants41. However, with the increase of phenolic acid secretion, some phenolic acids will penetrate into the soil and affect the soil microorganisms. Li et al.42 found that cinnamic acid inhibits Cylindrocarpon destructans (a pathogen of ginseng) growth at high concentrations, while promoting it at low concentrations. Yang et al.43 found that vanillic acid promoted the growth of the pathogens Rhizoctonia solani and Fusarium solani at low concentrations, but inhibited it at high concentrations; many phenolic acid compounds can inhibit the proliferation of Phytophthora cactorum (a pathogenic bacterium that causes AG phytophthora disease) at high concentrations. In addition, Yuan et al.44 found that p-coumaric acid strongly suppressed the in vitro growth of fungi, significantly reducing the decay caused by Alternaria alternata. Therefore, it can be seen that phenolic acids have inhibitory effects on pathogens at higher concentrations. With a decrease in soil phenolic acid content, this inhibitory effect on pathogenic bacteria will be weakened, resulting in an imbalance in the soil microbial composition that affects AG growth performance. Overall, soil phenolic acid content may indirectly affect AG growth performance by affecting soil microorganisms.
The change in the relative abundance of key bacteria
Our results showed that there was no significant difference in bacterial α-diversity between 10-year post-ginseng RS and CS, but there were differences in β-diversity, which reflects community composition and structure, between CS and RS. In other words, there were significant differences in the relative abundance of key bacteria in the bacterial community, such as Chlamydiae (phylum level, RS: 0.28%, CS: 0.10%, P = 0.035), within this phylum, the c_Chlamydiae, o_Chlamydiales, f_Simkaniaceae, and g_uncultured; Acidothermus (genus level, RS: 2.40%, CS: 5.40%, P = 0.030); Sphingomonadales (order level, CS: 2.98%, RS: 1.68%, P = 0.002), Sphingomonadaceae (family level, CS: 2.88%, RS: 1.48%, P = 0.004), genera Novosphingobium (CS: 0.03%, RS: 0.20%, P = 0.035) and Sphingomonas (CS: 2.83%, RS: 1.10%, P = 0.000); Rhodanobacter (CS: 0.38%, RS: 3.45%, P = 0.050); Arthrobacter (CS: 0.03%, RS: 0.43%, P = 0.001); Mizugakiibacter (CS: 0.63%, RS: 2.28%, P = 0.048); Jatrophihabitans (CS: 1.15%, RS: 0.75%, P = 0.048); Pseudomonas (RS: 0.15%, CS: 0.03%, P = 0.029) among others (Fig. 4, see Supplementary Table S2).
There was no difference in soil bacterial α-diversity between RS and CS, which may be due to the recovery of soil bacterial diversity after 10 years of rotation. However, the results of the pot experiment showed that RS still presented continuous cropping obstacles, which indicated that restoring soil microbial α-diversity does not alleviate continuous cropping obstacles for AG. Instead, differences in microbial community composition (i.e., β-diversity), particularly the abundances of bacterial taxa that play key roles, may explain the persistence of AG continuous cropping obstacles in RS after 10 years.
Among the differences in microbial community composition, CS had higher relative abundances of some bacterial genera that may be beneficial bacteria. The genus Acidothermus had the highest abundance, and it contained a single species, A. cellulolyticus, which is thermophilic, acidophilic, and has the ability to produce many thermostable cellulose-degrading enzymes45. Therefore, higher cellulose-degrading capacity might exist in CS than that in RS. Sphingomonas, a bacterium with the ability to decompose mono- and polycyclic aromatic compounds, as well as heterocyclic compounds, was more abundant in CS than RS, suggesting that bacterial decay of recalcitrant plant compounds was also higher in CS than RS. In addition, Sphingomonas not only decomposes monoaromatic phenolic acids but also improves plant stress resistance, and it is considered a plant probiotic46. Similar to our results, Li and Jiang23 found that Jatrophihabitans relative abundance in soil used for AG for 4 years was significantly (P  root rot group  > control group; in addition, compared with CS, there was a higher abundance of Rhodanobacter in the soil in which Korean ginseng (Panax ginseng) was grown49. We also found that this genus might be increased by the influence of Panax plants, which warrants further study. Our results showed that Arthrobacter was higher in the RS group, and Jiang et al.48 also found that the relative abundance of Arthrobacter in the root rot group was higher than that in the healthy root group; therefore, we speculate that Arthrobacter might be a factor causing root rot of P. quinquefolius, leading to a continuous cropping obstacle to AG growth. Our results showed that the abundance of Pseudomonas sp. in RS was higher than that in CS (RS: 0.15%, CS: 0.03%, P = 0.029, see Supplementary Table S2). Tan et al.50 showed that Pseudomonas sp. was the main pathogen causing root rot disease in P. notoginseng. In addition, Jiang et al.48 also found that Pseudomonas is abundant in the rhizosphere soils of diseased ginseng roots. Therefore, it is necessary to further study the effects of Pseudomonas species on AG growth. To sum up, the relative abundances of a large number of bacteria that are either confirmed or potentially harmful to other plants increased in RS, which may be an important factor leading to the occurrence of continuous cropping obstacles in the 10-year post-ginseng rotation soil.
As shown in Fig. 6, there are many correlations among the three factors. The abundances of Acidothermus, Sphingomonas, Jatrophihabitans, and Actinospica were each positively correlated with that of available K, caffeic acid, and cinnamic acid, but negatively correlated with that of salicylic acid. Therefore, the interactions among phenolic acids, microorganisms, and soil nutrients evidenced possible “synergistic” or “antagonistic” effects within the microecosystem. Overall, these complex relationships are the main reason for AG continuous cropping obstacles, but it is still unknown which of these factors plays the primary role. Finally, Nitrobacter, Actinospica, Clostridium sensu stricto 1, Thermosporothrix, Holophaga, and Peptoclostridium, also showed significant differences in abundance between RS and CS (Fig. 4), which also should receive more attention. More

No statistical methods were used to predetermine sample size. The experiments were not randomized, and investigators were not blinded to allocation during experiments and outcome assessment.
Etymology
The designation ‘Azoamicus’ combines the prefix azo- (New Latin, pertaining to nitrogen) with amicus (Latin, masculine noun, friend); thus giving azoamicus (‘friend that pertains to nitrogen’), alluding to its role as denitrifying endosymbiont. ‘Ciliaticola’ combines ciliate (referring to a group of ciliated protozoa) with the suffix -cola (derived from the Latin masculine noun incola, dweller, inhabitant), thus meaning ‘dwelling within a ciliate’.
Geochemical profiling
We carried out sampling for geochemical profiling in September 2016 and October 2018 at a single station located in the deep, southern lake basin of Lake Zug (about 197-m water depth) (47° 06′ 00.8′′ N, 8° 29′ 35.0′′ E). In September 2016, we used a multi-parameter probe to measure conductivity, turbidity, depth (pressure), temperature and pH (XRX 620, RBR). Dissolved oxygen was monitored online with normal and trace micro-optodes (types PSt1 and TOS7, Presens) with detection limits of 125 and 20 nM, respectively, and a response time of 7 s. In October 2018, we used a CTD (CTD60, Sea&Sun Technology) equipped with a Clark-type oxygen sensor (accuracy ± 3%, resolution 0.1%) to record oxygen profiles.
Sample collection
Water for bulk DNA and RNA analyses was collected in September 2016 and October 2018. Sample collection for DNA and RNA extraction in September 2016 has previously been described46. In October 2018, water was sampled with a Niskin bottle (Hydro-Bios) from 160 m, 170 m and 180 m. For each depth, 2 l of lake water was directly filtered on board our boat onto 0.22-μm Sterivex filter cartridges (Merck Millipore) using a peristaltic pump, subsequently purged with RNAlater preservation solution (Life Technologies) and stored at −20 °C until further processing. For fluorescence in situ hybridization (FISH) analyses, water from the same depths was fixed on board the boat with formaldehyde (1.5% final concentration; Electron Microscopy Sciences) and incubated in a chilled cool box for about 6 h before filtration onto 3-μm polycarbonate filters (Merck Millipore). Additional FISH samples using the same approach were collected in May 2019 from 189 m water depth.
Water for incubation experiments and single-ciliate PCR was sampled in May 2019 from 189-m water depth using a 10 l Go-Flo bottle (General Oceanics), filled into 2.5-l glass bottles without headspace, closed with butyl rubber stoppers and kept cold (at about 4 °C) and dark until further handling. During sampling, oxygen contamination was minimized by overflowing the bottle with anoxic lake water.
For combined FISH and differential interference contrast microscopy analyses, individual live ciliates were picked from lake water (from 186-m depth, collected February 2020) and directly fixed on microscope slides. In brief, microscope slides were treated with 0.1 mg ml−1 poly-l-lysine for 10 min at room temperature, washed with MilliQ water and dried. Ciliates were pre-enriched by gravity flow of bulk lake water through a 5-μm membrane filter and picked using a glass capillary under a binocular microscope. Picked ciliates were transferred into a droplet of formaldehyde (2% in 0.1× PBS, pH 7.6) on poly-l-lysine-coated microscope slides, incubated (for 1 h at room temperature) and washed with MilliQ water. FISH was performed as described in ‘Double-labelled oligonucleotide probe fluorescence in situ hybridization and microscopy’.
Nutrient measurements
Water samples for measurements of nutrients (ammonium, NOx and nitrite) were retrieved with a syringe sampler from 15 discrete depths at and below the base of the oxycline. Forty ml of water was directly injected into a 50-ml Falcon tube containing 10 ml of OPA reagent for fluorometric ammonium quantification47. In 2018, ammonium concentration was determined using the same method, except that the lake water was immediately sterile-filtered after sampling and frozen at −20 °C until further processing. For NOx quantification, 10 ml of water was sterile-filtered into a 15-ml Falcon tube and combined nitrate and nitrite concentration was determined by a commercial QuAAtro Segmented Flow Analyzer (SEAL Analytical).
Clone library construction and Sanger sequencing
‘Candidatus A. ciliaticola’-specific 16S rRNA gene primers were designed on the basis of the ‘Ca. A. ciliaticola’ circular metagenome-assembled genome sequence. Primers targeted the intergenic spacer regions about 50 bp up- and downstream of the 16S rRNA gene, resulting in a 1,568-bp-long PCR product. For clone library construction, the ‘Ca. A. ciliaticola’ 16S rRNA gene was amplified by a nested PCR approach from the same DNA extract used for metagenome sequencing obtained in September 2016 from 160-m water depth using the newly designed ‘Ca. A. ciliaticola’-specific primers (eub62A3_29F and eub62A3_1547R) followed by PCR amplification with general bacterial 16S rRNA gene primers (8F and 1492R) (Supplementary Table 2). Cloning and construction of the clone library is described in more detail in Supplementary Methods. Inserts of purified plasmids from five clones were sequenced by Sanger sequencing using the BigDye Terminator v.3.1 sequencing kit (Thermo Fisher Scientific) and primers M13f or M13r. The sequencing PCR contained 3 μl purified plasmid, 0.5 μl 10× sequencing buffer, 0.5 μl primer (10 μM) and 1 μl BigDye reagent. The PCR reactions were performed as follows: 99 cycles (1 °C s−1 ramp) of denaturation (10 s at 96 °C), annealing (5 s at 60 °C) and elongation (4 min at 60 °C). The PCR products were purified using gel filtration (Sephadex G-50 Superfine, Amersham Bioscience) followed by Sanger sequencing (3130xl genetic analyser, Applied Biosystems). The Sanger sequences were quality-trimmed and assembled using Sequencher v.5.4.6 and standard settings before trimming vector and primer sequences.
Probe design for fluorescence in situ hybridization
To visualize ‘Ca. A. ciliaticola’ cells in the environment, we designed a specific FISH probe on the basis of the ‘Ca. A. ciliaticola’ circular metagenome-assembled genome 16S rRNA gene sequence and closely related sequences within the clade eub62A3. The ‘Ca. A. ciliaticola’ 16S rRNA gene sequence was imported into Arb48 v.6.1 and aligned to the SILVA SSU Ref NR 99 132 database using the SINA-Aligner49. A FISH probe specific for ‘Ca. A. ciliaticola’ and most members of clade eub62A3 was designed (probe eub62A3_813 5′ CTAACAGCAAGTTTTCATCGTTTA 3′) (Supplementary Table 3) using the probe design tool implemented in Arb, and further manually refined and evaluated in silico using MathFISH50. The newly designed probe eub62A3_813 targets ‘Ca. A. ciliaticola’ and 78% of the ‘Ca. Azoamicus’ subgroup A and B sequences included in SILVA SSU Ref NR 99 138 (7 out of 9; the 2 sequences that are not targeted belong to ‘Ca. Azoamicus’ subgroup B), and shows no nontarget hits. Some sequences in the database had only 1 or 2 weak mismatches (98% identity with the 16S rRNA gene sequence of ‘Ca. A. ciliaticola’. The remaining reads shared >94% identity with ‘Ca. A. ciliaticola’ and all shared as top hits sequences belonging to ‘Ca. Azoamicus’ subgroup A when blasted against the NCBI nr database (accessed June 2020).
Clone-FISH
Clone-FISH was performed as previously described51. In brief, a purified plasmid containing the ‘Ca. A. ciliaticola’ 16S rRNA gene sequence (as described in ‘Clone library construction and Sanger sequencing’) in the correct orientation (confirmed by PCR using M13F and 1492R primers) was transformed into electrocompetent Escherichia coli JM109(DE3) cells (Promega) by electroporation using the Cell Porator and Voltage Booster System (Gibco) with settings 350 V, 330 μF capacitance, low ohm impedance, fast charge rate and 4 kΩ resistance (Voltage Booster). After electroporation, cells were transferred into SOC medium (Sigma Aldrich), incubated for 1 h at 37 °C and plated onto LB plates containing 100 mg l−1 kanamycin. After incubation overnight at 37 °C, 4 clones were picked and the presence of the insert was checked with PCR (primers M13F and 1492R) as described in ‘Clone library construction and Sanger sequencing’, followed by gel electrophoresis. An insert-positive clone was selected at random and grown in LB medium containing 100 mg l−1 kanamycin at 37 °C until optical density at 600 nm reached 0.37. Transcription of the plasmid insert was induced using isopropyl β-d-1-thiogalactopyranoside (IPTG) (1 mM final concentration). After addition of IPTG, cells were incubated for 1 h at 37 °C followed by addition of 170 mg l−1 chloramphenicol and subsequent incubation for 4 h. Cells were collected by centrifugation, fixed in 2% formaldehyde solution for 1 h at room temperature, washed and stored at 4 °C in phosphate-buffered saline (pH 7.4) containing 50% ethanol until further processing. Formamide melting curves52 were carried out using a ‘Ca. Azoamicus’-specific, HRP-labelled probe (eub62A3_813) (Extended Data Fig. 5). In brief, cells were applied to glass slides. Permeabilization with lysozyme, peroxidase inactivation, hybridization (10%, 30%, 35%, 40%, 45% and 50% formamide) and tyramide signal amplification (Oregon Green 488) was performed as previously described53. For each formamide concentration, images of three fields of view were recorded using the same exposure time for all formamide concentrations, which was optimized at 10% formamide all same settings using an Axio Imager 2 microscope (Zeiss) and analysed using Daime 2.154.
Double-labelled oligonucleotide probe fluorescence in situ hybridization and microscopy
Hybridization with double-labelled oligonucleotide probes (terminally double-labelled with Atto488 dye; details of the probe are in Supplementary Table 3) (Biomers) and counterstaining with DAPI was performed as previously described55. In brief, samples (either cut filter sections or ciliates picked and fixed on a glass slide, as described in ‘Sample collection’) were incubated in hybridization buffer containing 25% formamide and 5 ng DNA μl−1 probe (the same concentration was used for eub62A3_813 competitor 1 and 2) for 3 h at 46 °C, and subsequently washed in prewarmed washing buffer (5 mM EDTA, 159 mM NaCl) for 30 min at 48 °C. After a brief MilliQ water wash, samples were incubated for 5 min at room temperature in DAPI solution (1 μg ml−1), briefly washed in MilliQ water and air-dried. Filter sections were mounted onto glass slides. Samples were embedded in Prolong Gold Antifade Mountant (Thermo Fisher Scientific), and left to cure for 24 h. Confocal laser scanning and differential interference contrast microscopy were performed on a Zeiss LSM 780 (Zeiss, 63× oil objective, 1.4 numerical aperture, with differential interference contrast prism). Z-stack images were obtained to capture entire ciliate cells and fluorescence images of FISH probe and DAPI signals were projected into 2D for visualization. Cell counting was performed using a Axio Imager 2 microscope (Zeiss) in randomly selected fields of view (40× objective, grid length = 312.5 μm) on polycarbonate filters (3-μm pore size, 32-mm effective filter diameter; Merck Millipore) onto which 0.5 l PFA-fixed lake water was filtered.
For light microscopy, live ciliates were picked from Lake Zug water (May 2019, 189 m) and prepared for live ciliate imaging using light microscopy as previously described56 on an Axio Imager 2 microscope (Zeiss). For image acquisition and processing, Zeiss ZEN (blue edition) 2.3 was used.
Scanning electron microscopy
Ciliates sampled in February 2020 were individually picked under a binocular microscope, and washed in droplets of sterile-filtered lake water. Washed ciliates were then transferred into approximately 200 μl fixative on top of a polylysine-coated silicon wafer (0.1 mg ml−1 poly-l-lysine for 10 min at room temperature) and fixed for 1 h at room temperature. The fixative contained 2.5% glutaraldehyde (v/v, electron microscopy grade) in PHEM buffer57 (pH 7.4). Fixed ciliates attached to the silicon wafer were dehydrated in an ethanol series (30%, 50%, 70%, 80%, 96% and 100%) before automated critical-point drying (EM CPD300, Leica). Scanning electron microscopy was performed on a Quanta FEG 250 (FEI). Images were obtained using FEI xTM v.6.3.6.3123 at an acceleration voltage of 2 kV under high vacuum conditions and were captured using an Everhart–Thornley secondary electron detector. The image represents an integrated and drift corrected array of 128 images captured with a dwell time of 50 ns.
Single-ciliate PCR
Ciliates were picked from Lake Zug water (189 m, 2019) under the binocular microscope with a glass micropipette, and subsequently washed twice in drops of sterile nuclease-free water (Ambion) before being transferred into lysis buffer. DNA was extracted using MasterPure DNA purification kit (Ambion) following the manufacturer’s instructions with a final elution in 1× TE buffer (25 μl). Overall, DNA was extracted from four individual ciliates (S1–S4), five (C5) and ten pooled ciliates (C10) as well as no ciliate (control). 16S (‘Ca. A. ciliaticola’) and 18S rRNA gene sequences (ciliate host) were then separately amplified by PCR using primer pairs eub62A3_29F/_1547R and cil_384F/_1147R. PCR reactions (20 μl) with ‘Ca. A. ciliaticola’-specific primers (eub62A3_29F/_1547R) were performed as described in ‘Clone library construction and Sanger sequencing’ with following modifications: 5 μl template, 58 °C annealing temperature and 40 amplification cycles. PCR with ciliate-specific primers (cil_384F/_1147R) was performed analogously with following modifications: 55 °C annealing temperature, 50 s elongation and 35 amplification cycles. The PCR reactions with primer pairs eub62A3_29F/_1547R were further amplified in a second round of PCR (under the same conditions) using 2 μl PCR reaction from the first round. For each PCR step, successful amplification of products was checked using gel electrophoresis as described in Supplementary Methods. PCR reactions were subsequently purified using QIAquick PCR purification Kit (Qiagen) according to the manufacturer’s instructions with a final elution in sterile nuclease-free water (25 μl) (Ambion). Purified PCR reactions were subsequently sequenced using Sanger sequencing and processed as described in ‘Clone library construction and Sanger sequencing’ with the following modifications: primers eub62A3_29F, eub62A3_1547R (annealing temperature 58 °C) or cil_384F, cil_1147R (annealing temperature 55 °C). Two of the single-ciliate DNA extracts amplified with the endosymbiont-specific primers either showed a faint (S4) or no (S2) band and were not sequenced.
Nucleic acid extraction from lake water
Bulk DNA and RNA extraction as well as metagenome and metatranscriptome sequencing of lake water samples collected in 2016 have previously been described46. For samples from 2018, filters were purged of RNAlater, briefly rinsed with nuclease-free water and removed from the Sterivex cartridge. RNA and DNA was then extracted from separate filters using PowerWater RNA or DNA isolation kits (MoBio Laboratories) according to the manufacturer’s instructions.
Metagenome and bulk metatranscriptome sequencing
For metagenomic sequencing, DNA libraries were prepared as recommended by the NEBNext Ultra II FS DNA Library Prep Kit for Illumina (New England Biolabs). Sequencing-by-synthesis was performed on the Illumina HiSeq2500 sequencer (Illumina Inc.) with the 2 × 250-bp read mode. For metatranscriptomic sequencing, rRNA was removed (Ribo-Zero rRNA Removal Kit for bacteria (Illumina)) and an RNA-sequencing library was prepared according to the NEBNext Ultra II Directional RNA Library Prep Kit for Illumina (New England Biolabs). Sequencing-by-synthesis was performed on the Illumina HiSeq3000 sequencer (Illumina) with the 1 × 150-bp read mode. Library preparation and sequencing was performed by the Max Planck Genome Centre Cologne (http://mpgc.mpipz.mpg.de/home/). Detailed information for each metagenomic and metatranscriptomic dataset can be found in Supplementary Table 4.
Genome assembly and finishing
The genome of ‘Ca. A. ciliaticola’ was reconstructed from a metagenomic dataset sampled in 2018, as follows: metagenomic reads (MG_18_C) were trimmed using Trimmomatic58 v.0.32 as previously described46 and assembled using metaSPAdes59 v.3.13.0 and k-mer lengths of 21, 33, 55, 77, 99 and 127. From this assembly, contigs with high similarity to the previously reconstructed ‘Ca. A. ciliaticola’ genome from the metagenome samples from 2016 (further details are provided in Supplementary Methods) were identified by blastn (identity >95%) and metagenomic reads were mapped back to the contigs using BBmap60 v.35.43 (minid = 0.98). Mapped reads were subsequently reassembled using SPAdes v.3.13.0 with mismatch corrector and coverage threshold enabled (–careful –cov-cutoff 60), resulting in the assembly of a single contig (292,647 bp) that was circularized by trimming the identical overlapping ends (127 bp) giving rise to the closed genome (292,520 bp). The start position was set in an intergenic spacer region near the maximum of the GC disparity curve generated with oriFinder61 v.1.0. The two independently assembled circular metagenome-assembled genomes (from 2016 and 2018 metagenomes (Supplementary Methods)) shared 99.99% identity. For all subsequent analyses, the genome reconstructed from the 2018 dataset was used owing to higher coverage.
Genome annotation and comparative analyses
Genome annotation was performed using a modified version of Prokka62 v.1.13.3 to allow annotation of genes that overlap with tRNA genes. The annotation of key metabolic genes was manually inspected and refined using searches against NCBI non-redundant protein or conserved domain database63. Transmembrane transporters were predicted and classified using the Transporter Automatic Annotation Pipeline web service64 and the Transporter Classification Database65. Pseudogenes were predicted using pseudo finder66 v.0.11 and standard settings. Circular genome maps were created using DNAplotter67 v.18.1.0.
For comparative analyses, protein-coding CDS encoded in the genomes of insect endosymbionts (C. ruddii PV, AP009180.1 and B. aphidicola BCc, CP000263.1), mitochondrion of J. libera (NC_021127) and a free-living relative of ‘Ca. A. ciliaticola’ (L. clemsonensis, CP016397.1) were downloaded from NCBI GenBank. Additionally, protein-coding CDS of the ciliate endosymbiont C. taeniospiralis (PGGB00000000.1) were obtained using Prokka annotation. Classification of functional categories was performed using the eggNOG-mapper v.1 web service68 with mapping mode DIAMOND and standard settings. The classification of the functional category C (energy production and conversion) for ‘Ca. A. ciliaticola’ was modified to also include norB and nirK, which were grouped by eggNOG into different categories (Q and P, respectively).
Multiple sequence alignment of ‘Ca. A. ciliaticola’ and other plastidic and bacterial ATP/ADP translocases was generated using MuscleWS69 v.3.8.31 with default settings implemented in Jalview70 v.2.11.1.0.
Metatranscriptomic analyses of bulk water samples
Raw metatranscriptomic Illumina reads trimming and removal of rRNA sequences was performed as previously described46. Non-rRNA reads were then mapped to the genome of ‘Ca. A. ciliaticola’ using Bowtie271 v.2.2.1.0 and standard parameters. Sorted and indexed BAM files were generated using samtools72 v.0.1.19 and transcripts per feature (based on the Prokka annotation) were quantified using EDGE-pro73 v.1.3.1 and standard settings. Gene transcription was subsequently quantified as transcripts per million74 (TPM) using the formula:

$${{rm{T}}{rm{P}}{rm{M}}}_{{rm{i}}}=frac{{c}_{i}}{{l}_{i}}times frac{1}{{sum }_{j}frac{{c}_{i}}{{l}_{i}}}times {10}^{6}$$

to assign each feature (i) a TPM value, in which c = feature count, l = length (in kb) and j = all features.
Read coverage visualization and plotting was performed using pyGenomeTracks75 (average coverage over 100-bp bins) implemented in deepTools276 v.3.2.0.
Phylogenetic analyses
The full-length 16S rRNA gene sequence was retrieved from the circular metagenome-assembled genome of ‘Ca. A. ciliaticola’ using RNAmmer77 v.1.5, aligned using the SILVA incremental aligner49 (SINA 1.2.11) and imported to the SILVA SSU NR99 database45 (release 132) using ARB48 v.6.1. Additional closely related 16S rRNA gene sequences were identified by BLASTN in the NCBI non-redundant nucleotide database and JGI IMG/M 16S rRNA Public Assembled Metagenomes (retrieved July 2018) and also imported into ARB. A maximum-likelihood phylogenetic tree of 16S rRNA gene sequences was calculated using RAxML78 v.8.2.8 integrated in ARB with the GAMMA model of rate heterogeneity and the GTR substitution model with 100 bootstraps. The alignment was not constrained by a weighting mask or filter. For the complete tree shown in Extended Data Fig. 4, additional ‘Ca. A. ciliaticola’ sequences obtained from a clone library and individual single ciliates were added to the tree using the Parsimony ‘Quick add’ algorithm implemented in ARB.
For the ciliate phylogeny, sequences obtained from Sanger sequencing of picked ciliates were added to the EukRef-Ciliphora30 Plagiopylea subgroup alignment using MAFFT79 online service version 7 (argument:–addfragments). An additional metagenome-assembled full-length 18S rRNA gene sequence assigned to Plagiopylea was obtained using phyloFlash80 v.3.0 from one of the Lake Zug metagenomes (MG_18_C) and also added to the alignment (argument:–add). A phylogenetic tree was calculated using IQ-TREE webserver (http://iqtree.cibiv.univie.ac.at) running IQ-TREE81 1.6.11 with default settings and automatic substitution model selection (best-fit model: TIM2+F+I+G4). Phylogenetic trees were visualized using the Interactive Tree of Life v.4 web service82.
For the ATP/ADP translocase phylogenetic tree, amino acid sequences were retrieved from NCBI RefSeq (250 top hits) and NCBI nr (15 top hits) (both accessed June 2019) using NCBI blastp web-service83 with the amino acid sequence of ATP/ADP translocase sequence of ‘Ca. A. ciliaticola’ (ESZ_00147) as query. Additional amino acid sequences of characterized nucleotide transporters listed in Supplementary Table 8 were also added. After removal of duplicates, sequences were clustered at 95% identity using usearch84 v.8.0.1623 and aligned using MUSCLE69 v.3.8.31. Phylogenetic tree construction using IQ-TREE (best-fit model: LG+F+I+G4) and visualization was performed as described for the 18S rRNA gene phylogenetic tree.
Incubation experiments
Incubation experiments were performed to provide experimental evidence for the denitrifying activity linked to the ciliate host. Three incubations were set up that contained (a) no ciliates (filtered fraction), (b) lake water that was enriched in ciliates (enriched fraction) and (c) bulk lake water. For these experiments, lake water was size-fractionated using a 10-μm polycarbonate filter (Merck Millipore) under N2 atmosphere in a glove bag at 12 °C. Enriched and filtered fractions were obtained by gravity filtration of 0.5 l water until 0.25 l supernatant was left. Thus, in the enriched fraction, ciliates from 0.5 l lake water were concentrated in 0.25 l lake water. In the filtered fraction, organisms larger than 10 μm (including ciliates) were filtered out. Both the enriched water (plus filter) and the filtered water were transferred to separate serum bottles (no headspace) and closed with butyl rubber stoppers. For bulk incubations, unfiltered water was directly filled into 250 ml serum bottles under N2 atmosphere.
Denitrification potential was assessed by measuring the production of 30N2 over time in 15N-nitrite and 15N-nitrate amended incubations by isotope ratio mass spectrometry (Isoprime Precision running ionOS v.4.04, Elementar). A 30 ml helium headspace was set for the 250 ml serum bottles and the water was degassed with helium for 10 min to ensure anoxic conditions and reduce N2 background. A 15N-labelled mixture of nitrate and nitrite (20 μM and 5 μM final concentration, respectively; Sigma Aldrich) was supplied at a 99% labelling percentage and the water was incubated for a total of 30 h at 4 °C in the dark. Subsamples of the headspace were taken at regular time intervals during the incubation by withdrawing 3 ml of the gaseous headspace and simultaneously replacing the removed volume by helium. This gas sample was transferred into 12 ml Exetainers (LabCo) that were pre-filled with helium-degassed water and stored until analysis. 30N2 in the gas samples was measured on an isotope ratio mass spectrometer, and denitrification rates were calculated from the slope of the linear increase of 30N2 in the headspace over the time course of the experiments. The rate of 30N2 production was corrected for dilution of the headspace introduced by subsampling and by the measured total 15N labelling percentage. ‘Ca. A. ciliaticola’-containing ciliate abundance in the different incubation bottles was assessed by microscopic counts after cell fixation, FISH hybridization (eub62A3_813) and DAPI staining as described in ‘Double-labelled oligonucleotide probe fluorescence in situ hybridization and microscopy’.
Statistics and reproducibility
No statistical methods were used to predetermine sample size and experiments were not randomized. The investigators were not blinded to allocation during experiments and outcome assessment.
In Fig. 1a, the scanning electron microscopy image is a representative of n = 6 recorded images that were obtained from 1 experiment. In Fig 1b, the differential interference contrast image is a representative of n = 6 recorded images that were obtained from 1 experiment.
In Fig. 2c, Extended Data Fig. 2i. FISH fluorescence images (eub62A3_813 probe) are representative of n = 33 recorded images that were obtained from 5 independent experiments of 3 biological replicate samples.
In Fig. 2c, Extended Data Fig. 2f, h, j. DAPI fluorescence images are representative of n = 21 recorded images that were obtained from 5 independent experiments of 3 biological replicate samples.
In Extended Data Fig. 2a, the FISH fluorescence image (Arch915 probe) is representative of n = 6 images that were obtained from 1 experiment. In Extended Data Fig. 2b, d, F420 autofluorescence images are representative of n = 11 recorded images that were obtained from 3 independent experiments of 1 sample. In Extended Data Fig. 2c, g, FISH fluorescence images (EUB-I probe) are representative of n = 15 images obtained from 3 independent experiments of 2 biological replicate samples. In Extended Data Fig. 2e, the FISH fluorescence image (NON338 probe) is representative of n = 15 recorded images that were obtained from 3 independent experiments of 2 biological replicate samples.
In Extended Data Fig. 5a, each of the 6 FISH fluorescence images (eub62A3_813 probe) is representative of n = 3 images from 1 experiment.
For the fluorescence images shown, the number of analysed cells was typically much higher (n  > 100) than the ones that were eventually recorded.
Reporting summary
Further information on research design is available in the Nature Research Reporting Summary linked to this paper. More

1.
Jost C, Lawrence CA, Campolongo F, Van De Bund W, Hill S, DeAngelis DL. The effects of mixotrophy on the stability and dynamics of a simple planktonic food web model. Theor Popul Biol. 2004;66:37–51.
PubMed  Article  Google Scholar 
2.
Tittel J, Bissinger V, Zippel B, Gaedke U, Bell E, Lorke A, et al. Mixotrophs combine resource use to outcompete specialists: Implications for aquatic food webs. Proc Natl Acad Sci. 2011;100:12776–81.
Article  CAS  Google Scholar 

3.
Ward BA, Follows MJ. Marine mixotrophy increases trophic transfer efficiency, mean organism size, and vertical carbon flux. Proc Natl Acad Sci. 2016;113:2958–63.
CAS  PubMed  Article  Google Scholar 

4.
Hansen PJ, Tillmann U. Mixotrophy among dinoflagellates—prey selection, physiology and ecological imporance. In: Subba Rao DV, editor. Dinoflagellates: classification, evolution, physiology and ecological significance. Hauppauge, NY, USA: Nova; 2020;201–60.

5.
Unrein F, Gasol JM, Not F, Forn I, Massana R. Mixotrophic haptophytes are key bacterial grazers in oligotrophic coastal waters. ISME J. 2014;8:164–76.
CAS  PubMed  Article  Google Scholar 

6.
Anderson R, Charvet S, Hansen P. Mixotrophy in chlorophytes and haptophytes – effect of irradiance, macronutrient, micronutrient and vitamin limitation. Front Microbiol. 2018;9:1704.
PubMed  PubMed Central  Article  Google Scholar 

7.
Lewitus AJ, Caron DA, Miller KR. Effect of light and glycerol on the organization of the photosynthetic apparatus in the facultative heterotroph Pyrenomonas salina (cryptophyceae). J Phycol. 1991;27:578–87.
Article  Google Scholar 

8.
Du YooY, Seong KA, Jeong HJ, Yih W, Rho J-R, Nam SW, et al. Mixotrophy in the marine red-tide cryptophyte Teleaulax amphioxeia and ingestion and grazing impact of cryptophytes on natural populations of bacteria in Korean coastal waters. Harmful Algae. 2017;68:105–17.
Article  Google Scholar 

9.
Caron DA, Porter KG, Sanders RW. Carbon, nitrogen, and phosphorus budgets for the mixotrophic phytoflagellate Poterioochromonas malhamensis (Chrysophyceae) during bacterial ingestion. Limnol Oceanogr. 1990;35:433–43.
CAS  Article  Google Scholar 

10.
Holen DA, Boraas ME. Mixotrophy in chrysophytes. Chrysophyte algae. Cambridge, UK: Cambridge University Press; 1995;119–40.

11.
Fenchel T. Ecology of heterotrophic microflagellates. II. Bioenerg growth Mar Ecol Prog Ser. 1982;8:225–31.
Article  Google Scholar 

12.
Rottberger J, Gruber A, Boenigk J, Kroth P. Influence of nutrients and light on autotrophic, mixotrophic and heterotrophic freshwater chrysophytes. Aquat Micro Ecol. 2013;71:179–91.
Article  Google Scholar 

13.
Bell EM, Laybourn-Parry J. Mixotrophy in the antarctic phytoflagellate, Pyramimonas gelidicola (Chlorophyta: Prasinophyceae). J Phycol. 2003;39:644–9.
Article  Google Scholar 

14.
McKie-Krisberg ZM, Sanders RW. Phagotrophy by the picoeukaryotic green alga Micromonas: implications for Arctic Oceans. ISME J. 2014;8:1953–61.
CAS  PubMed  PubMed Central  Article  Google Scholar 

15.
McKie-Krisberg ZM, Gast RJ, Sanders RW. Physiological responses of three species of Antarctic mxotrophic phytoflagellates to changes in light and dissolved nutrients. Micro Ecol. 2015;70:21–29.
CAS  Article  Google Scholar 

16.
Paasch A. Physiological and genomic characterization of phagocytosis in green algae. New York, NY, USA: American Museum of Natural History; 2017.

17.
Not F, Latasa M, Scharek R, Viprey M, Karleskind P, Balagué V, et al. Protistan assemblages across the Indian Ocean, with a specific emphasis on the picoeukaryotes. Deep Res Part I Oceanogr Res Pap. 2008;55:1456–73.
Article  Google Scholar 

18.
Shi XL, Marie D, Jardillier L, Scanlan DJ, Vaulot D. Groups without cultured representatives dominate eukaryotic picophytoplankton in the oligotrophic South East Pacific Ocean. PLoS ONE. 2009;4:e7657.
PubMed  PubMed Central  Article  CAS  Google Scholar 

19.
Rii YM, Duhamel S, Bidigare RR, Karl DM, Repeta DJ, Church MJ. Diversity and productivity of photosynthetic picoeukaryotes in biogeochemically distinct regions of the South East Pacific Ocean. Limnol Oceanogr. 2016;61:806–24.
Article  Google Scholar 

20.
Maruyama S, Kim E. A modern descendant of early green algal phagotrophs. Curr Biol. 2013;23:1081–4.
CAS  PubMed  Article  Google Scholar 

21.
O’Kelly C. Flagellar apparatus architecture and the phylogeny of ‘green’ algae: Chlorophytes, Euglenoids, Glaucophytes. In: Menzel D, editor. The cytoskeleton of the algae. Boca Raton: CRC Press; 1992. p. 315–41.
Google Scholar 

22.
Burns JA, Pittis AA, Kim E. Gene-based predictive models of trophic modes suggest Asgard archaea are not phagocytotic. Nat Ecol Evol. 2018;2:697–704.
PubMed  Article  Google Scholar 

23.
Wilken S, Yung CCM, Hamilton M, Hoadley K, Nzongo J, Eckmann C, et al. The need to account for cell biology in characterizing predatory mixotrophs in aquatic environments. Philos Trans R Soc B Biol Sci. 2019;374:20190090.
CAS  Article  Google Scholar 

24.
Inouye I, Hori T, Chihara M. Absolute configuration analysis of the flagellar apparatus of Pterosperma Cristatum (Prasinophyceae) and consideration of Its phylogenetic position. J Phycol. 1990;26:329–44.
Article  Google Scholar 

25.
Bhuiyan MAH, Faria DG, Horiguchi T, Sym SD, Suda S. Taxonomy and phylogeny of Pyramimonas vacuolata sp. nov. (Pyramimonadales, Chlorophyta). Phycologia. 2015;54:323–32.
CAS  Article  Google Scholar 

26.
Adl SM, Bass D, Lane CE, Lukeš J, Schoch CL, Smirnov A, et al. Revisions to the classification, nomenclature, and diversity of eukaryotes. J Eukaryot Microbiol. 2019;66:4–119.
PubMed  PubMed Central  Article  Google Scholar 

27.
Burns JA, Paasch A, Narechania A, Kim E. Comparative genomics of a bacterivorous green alga reveals evolutionary causalities and consequences of phago-mixotrophic mode of nutrition. Genome Biol Evol. 2015;7:3047–61.
CAS  PubMed  PubMed Central  Article  Google Scholar 

28.
Guillard R. Culture of phytoplankton for feeding marine invertebrates. In: Smith WL, Chanley WH, editors. Culture of marine invertebrate animals. 1975. New York: Plenum Press; 1975. p. 22–60.

29.
Cho J-C, Giovannoni SJ. Pelagibaca bermudensis gen. nov., sp. nov., a novel marine bacterium within the Roseobacter clade in the order Rhodobacterales. Int J Syst Evol Microbiol. 2006;56:855–9.
CAS  PubMed  Article  Google Scholar 

30.
Thrash JC, Cho J-C, Ferriera S, Johnson J, Vergin KL, Giovannoni SJ. Genome sequences of Pelagibaca bermudensis HTCC2601T and Maritimibacter alkaliphilus HTCC2654T, the type strains of two marine Roseobacter genera. J Bacteriol. 2010;192:5552–3.
CAS  PubMed  PubMed Central  Article  Google Scholar 

31.
First MR, Park NY, Berrang ME, Meinersmann RJ, Bernhard JM, Gast RJ, et al. Ciliate ingestion and digestion: Flow cytometric measurements and regrowth of a digestion-resistant Campylobacter jejuni. J Eukaryot Microbiol. 2012;59:12–19.
PubMed  Article  Google Scholar 

32.
Sherr BF, Sherr EB, Fallon RD. Use of monodispersed, fluorescently labeled bacteria to estimate in situ protozoan bacterivory. Appl Environ Microbiol. 1987;53:958–65.
CAS  PubMed  PubMed Central  Article  Google Scholar 

33.
Vazquez-Dominguez E, Peters F, Gasol JM, Vaqué D. Measuring the grazing losses of picoplankton: methodological improvements in the use of fluorescently labeled tracers combined with flow cytometry. Aquat Micro Ecol. 1999;20:119–28.
Article  Google Scholar 

34.
Leebens-Mack J, Barker M, Carpenter EJ. One thousand plant transcriptomes and the phylogenomics of green plants. Nature. 2019;574:679–85.
Article  CAS  Google Scholar 

35.
Wincker P. A thousand plants’ phylogeny. Nat Plants. 2019;5:1106–7.
PubMed  Article  PubMed Central  Google Scholar 

36.
Keeling PJ, Burki F, Wilcox HM, Allam B, Allen EE, Amaral-Zettler LA, et al. The marine microbial eukaryote transcriptome sequencing project (MMETSP): illuminating the functional diversity of eukaryotic life in the oceans through transcriptome sequencing. PLOS Biol. 2014;12:e1001889.
PubMed  PubMed Central  Article  CAS  Google Scholar 

37.
Johnson LK, Alexander H, Brown CT. Re-assembly, quality evaluation, and annotation of 678 microbial eukaryotic reference transcriptomes. Gigascience. 2019;8:1–12.
Google Scholar 

38.
Besemer J. GeneMarkS: a self-training method for prediction of gene starts in microbial genomes. Implications for finding sequence motifs in regulatory regions. Nucleic Acids Res. 2001;29:2607–18.
CAS  PubMed  PubMed Central  Article  Google Scholar 

39.
Simão FA, Waterhouse RM, Ioannidis P, Kriventseva EV, Zdobnov EM. BUSCO: assessing genome assembly and annotation completeness with single-copy orthologs. Bioinformatics. 2015;31:3210–2.
PubMed  PubMed Central  Google Scholar 

40.
Gawryluk RMR, Tikhonenkov DV, Hehenberger E, Husnik F, Mylnikov AP, Keeling PJ. Non-photosynthetic predators are sister to red algae. Nature. 2019;572:240–3.
CAS  PubMed  Article  Google Scholar 

41.
Newcombe RG. Interval estimation for the difference between independent proportions: comparison of eleven methods. Stat Med. 1998;17:873–90.
CAS  PubMed  Article  Google Scholar 

42.
Kursa M, Rudnicki W. Feature selection with the Boruta package. J Stat Softw. 2010;36:1–13.
Article  Google Scholar 

43.
Chasset PO. Probabilistic neural network for the R statistical language. https://github.com/chasset/pnn. Github. 2013.

44.
Maia R, Eliason CM, Bitton P-P, Doucet SM, Shawkey MD. pavo: an R package for the analysis, visualization and organization of spectral data. Methods Ecol Evol. 2013;4:906–13.
Google Scholar 

45.
Jimenez V, Burns J, Le Gall F, Not F, Vaulot D. No evidence of phago-mixotropy in Micromonas polaris, the dominant picophytoplankton species in the Arctic. J Phycol. 2021. https://doi.org/10.1111/jpy.13125.

46.
R Core Team. R development core team. R A Lang Environ Stat Comput. Vienna: R Foundation for Statistical Computing; 2016.

47.
Figueroa-Martinez F, Nedelcu AM, Smith DR, Reyes-Prieto A. When the lights go out: the evolutionary fate of free-living colorless green algae. N. Phytol. 2015;206:972–82.
Article  Google Scholar 

48.
Nakada T, Misawa K, Nozaki H. Molecular systematics of Volvocales (Chlorophyceae, Chlorophyta) based on exhaustive 18S rRNA phylogenetic analyses. Mol Phylogenet Evol. 2008;48:281–91.
CAS  PubMed  Article  Google Scholar 

49.
Johnson I. The molecular probes handbook: a guide to fluorescent probes and labeling technologies. 11th ed. Waltham, MA, USA: Life Technologies Corporation; 2010.

50.
Leliaert F, Smith DR, Moreau H, Herron MD, Verbruggen H, Delwiche CF, et al. Phylogeny and molecular evolution of the green algae. CRC Crit Rev Plant Sci. 2012;31:1–46.
Article  Google Scholar 

51.
Leliaert F. Green algae: chlorophyta and streptophyta. Reference module in life sciences. Amsterdam, DK: Elsevier; 2019.

52.
Parke M, Adams I. The Pyramimonas-like motile stage of Halosphaera viridis Schmitz. Bull Res Counc Isr. 1961.

53.
Thorndsen J. Cymbomonas Schiller (Prasinophyceae) reinvestigated by light and electron microscopy. Arch fur Protistenkd. 1988;136:327–36.
Article  Google Scholar 

54.
González JM, Sherr BF, Sherr EB. Digestive enzyme activity as a quantitative measure of protistan grazing: the acid lysozyme assay for bacterivory. Mar Ecol Prog Ser. 1993;100:197–206.
Article  Google Scholar 

55.
Moestrup Ø, Inouye I, Hori T. Ultrastructural studies on Cymbomonas tetramitiformis (Prasinophyceae). I. General structure, scale microstructure, and ontogeny. Can J Bot. 2003;81:657–71.
Article  Google Scholar 

56.
Turmel M, Lopes dos Santos A, Otis C, Sergerie R, Lemieux C. Tracing the evolution of the plastome and mitogenome in the Chloropicophyceae uncovered convergent tRNA gene losses and a variant plastid genetic code. Genome Biol Evol. 2019;11:1275–92.
CAS  PubMed  PubMed Central  Article  Google Scholar 

57.
Lopes dos Santos A, Gourvil P, Tragin M, Noël M, Decelle J, Romac S, et al. Diversity and oceanic distribution of prasinophytes clade VII, the dominant group of green algae in oceanic waters. ISME J. 2017;11:512–28.
PubMed  Article  Google Scholar 

58.
Lemieux C, Turmel M, Otis C, Pombert J-F. A streamlined and predominantly diploid genome in the tiny marine green alga Chloropicon primus. Nat Commun. 2019;10:4061.
PubMed  PubMed Central  Article  Google Scholar 

59.
Zingone A, Borra M, Brunet C, Forlani G. Kooistra WHCF, Procaccini G. Phylogenetic position of Crustomatix stigmatica sp. nov. and Dolichomastix tenuilepis in relation to the mamiellales (Prasinophyceae, Chlorophyta). J Phycol. 2002;38:1024–39.
CAS  Article  Google Scholar 

60.
Liang Z, Geng Y, Ji C, Du H, Wong CE, Zhang Q, et al. Mesostigma viride genome and transcriptome provide insights into the origin and evolution of Streptophyta. Adv Sci. 2020;7:1901850.
CAS  Article  Google Scholar 

61.
Buckley CM, Gopaldass N, Bosmani C, Johnston SA, Soldati T, Insall RH, et al. WASH drives early recycling from macropinosomes and phagosomes to maintain surface phagocytic receptors. Proc Natl Acad Sci. 2016;113:E5906–15.
CAS  PubMed  Article  Google Scholar 

62.
Shpak M, Kugelman JR, Varela-Ramirez A, Aguilera RJ. The phylogeny and evolution of deoxyribonuclease II: An enzyme essential for lysosomal DNA degradation. Mol Phylogenet Evol. 2008;47:841–54.
CAS  PubMed  Article  Google Scholar 

63.
Gast RJ, McKie-Krisberg ZM, Fay SA, Rose JM, Sanders RW. Antarctic mixotrophic protist abundances by microscopy and molecular methods. FEMS Microbiol Ecol. 2014;89:388–401.
CAS  PubMed  Article  Google Scholar 

64.
Mitra A, Flynn KJ, Tillmann U, Raven JA, Caron D, Stoecker DK, et al. Defining planktonic protist functional groups on mechanisms for energy and nutrient acquisition: Incorporation of diverse mixotrophic strategies. Protist. 2016;167:106–20.
CAS  Article  Google Scholar 

65.
Kirkham AR, Lepère C, Jardillier LE, Not F, Bouman H, Mead A, et al. A global perspective on marine photosynthetic picoeukaryote community structure. ISME J. 2013;7:922–36.
CAS  PubMed  PubMed Central  Article  Google Scholar 

66.
Moon-van Der Staay SY, Wachter R De, Vaulot D. Oceanic 18S rDNA sequences from picoplankton reveal unsuspected eukaryotic diversity. Nature. 2001;409:607–10.
CAS  PubMed  Article  Google Scholar 

67.
Worden A. Picoeukaryote diversity in coastal waters of the Pacific Ocean. Aquat Micro Ecol. 2006;43:165–75.
Article  Google Scholar 

68.
Van Hannen EJ, Veninga M, Bloem J, Gons HJ, Laanbroek HJ. Genetic changes in the bacterial community structure associated with protistan grazers. Fundam Appl Limnol. 1999;145:25–38.
Article  Google Scholar 

69.
Jürgens K, Güde H. The potential importance of grazing-resistant bacteria in planktonic systems. Mar Ecol Prog Ser. 1994;112:169–88.
Article  Google Scholar 

70.
Jürgens K, Pernthaler J, Schalla S, Amann R. Morphological and compositional changes in a planktonic bacterial community in response to enhanced protozoan grazing. Appl Environ Microbiol. 1999;65:1241–50.
PubMed  PubMed Central  Article  Google Scholar 

71.
Suzuki M. Effect of protistan bacterivory on coastal bacterioplankton diversity. Aquat Micro Ecol. 1999;20:261–72.
Article  Google Scholar 

72.
Sherr EB, Sherr BF. Significance of predation by protists in aquatic microbial food webs. Antonie van Leeuwenhoek2. 2002;81:293–308.
CAS  Article  Google Scholar 

73.
González J, Sherr EB, Sherr BF. Size-selective grazing on bacteria by natural assemblages of estuarine flagellates and ciliates. Appl Environ Microbiol. 1990;56:583–9.
PubMed  PubMed Central  Article  Google Scholar 

74.
Sherr BF, Sherr EB, McDaniel J. Effect of protistan grazing on the frequency of dividing cells in bacterioplankton assemblages. Appl Environ Microbiol. 1992;58:2381–5.
CAS  PubMed  PubMed Central  Article  Google Scholar 

75.
González J, Sherr EB, Sherr BF. Differential feeding by marine flagellates on growing vs starving bacteria, and on motile vs non-motile bacteria. Mar Ecol Prog Ser. 1993;102:257–67.
Article  Google Scholar 

76.
del Giorgio PA, Gasol JM, Vaqué D, Mura P, Agustí S, Duarte CM. Bacterioplankton community structure: Protists control net production and the proportion of active bacteria in a coastal marine community. Limnol Oceanogr. 1996;41:1169–79.
Article  Google Scholar 

77.
Andersen OK, Goldman JC, Caron DA, Dennett MR. Nutrient cycling in a microflagellate food chain: III. Phosphorus dynamics. Mar Ecol Prog Ser. 1986;31:47–55.
CAS  Article  Google Scholar 

78.
Fenchel T. Protistan filter feeding. Prog Protistol. 1986;1:65–113.
Google Scholar 

79.
Epstein S, Shiaris M. Size selective grazing of coastal bacterioplankton by natural assemblages of pigmented flagellates, colourless flagellates and ciliates. Micro Ecol. 1992;23:211–25.
CAS  Article  Google Scholar 

80.
Montagnes D, Barbosa A, Boenigk J, Davidson K, Jurgens K, Macek M, et al. Selective feeding behaviour of key free-living protists: avenues for continued study. Aquat Micro Ecol. 2008;53:83–98.
Article  Google Scholar 

81.
Pfister G, Arndt H. Food selectivity and feeding behaviour in omnivorous filter-feeding ciliates: a case study for Stylonychia. Eur J Protistol. 1998;34:446–57.
Article  Google Scholar 

82.
Boenigk J, Arndt H. Bacterivory by heterotrophic flagellates: community structure and feeding strategies. Antonie Van Leeuwenhoek. 2002;81:465–80.
PubMed  Article  Google Scholar 

83.
Pickup ZL, Pickup R, Parry JD. Growth of Acanthamoeba castellanii and Hartmannella vermiformis on live, heat-killed and DTAF-stained bacterial prey. FEMS Microbiol Ecol. 2007;61:264–72.
CAS  PubMed  Article  Google Scholar 

84.
Legrand C, Johansson N, Johnsen G, Borsheim K, Graneli E. Phagotrophy and toxicity variation in mixotrophic Prymnesium patelliferum (Haptophyceae). Limnol Oceanogr. 2001;46:1208–14.
Article  Google Scholar 

85.
Caron DA, Sanders RW, Lim EL, Marrasé C, Amaral LA, Whitney S, et al. Light-depend phagotrophy freshwater mixotrophic chrysophyte Dinobryon cylindricum. Micro. Ecol. 1993;25:93–111.
CAS  Article  Google Scholar 

86.
Fenchel T. The microbial loop – 25 years later. J Exp Mar Bio Ecol. 2008;366:99–103.
Article  Google Scholar 

87.
Tittel J, Bissinger V, Zippel B, Gaedke U, Bell E, Lorke A, et al. Mixotrophs combine resource use to outcompete specialists: Implications for aquatic food webs. Proc Natl Acad Sci. 2003;100:12776–81.
CAS  PubMed  Article  Google Scholar 

88.
Moorthi S, Ptacnik R, Sanders R, Fischer R, Busch M, Hillebrand H. The functional role of planktonic mixotrophs in altering seston stoichiometry. Aquat Micro Ecol. 2017;79:235–45.
Article  Google Scholar 

89.
Katechakis A, Haseneder T, Kling R, Stibor H. Mixotrophic versus photoautotrophic specialist algae as food for zooplankton: The light: nutrient hypothesis might not hold for mixotrophs. Limnol Oceanogr. 2005;50:1290–9.
CAS  Article  Google Scholar 

90.
Weisse T, Anderson R, Arndt H, Calbet A, Hansen PJ, Montagnes D. Functional ecology of aquatic phagotrophic protists – concepts, limitations, and perspectives. Eur J Protistol. 2016;55:50–74.
PubMed  Article  Google Scholar 

91.
Graham LE, Graham JM, Wilcox WL, Cook ME. Algae. 3rd ed. Madison, WI, USA: LJLM Press; 2016.

92.
Guillou L, Eikrem W, Chrétiennot-Dinet M-J, Le Gall F, Massana R, Romari K, et al. Diversity of picoplanktonic prasinophytes assessed by direct nuclear SSU rDNA sequencing of environmental samples and novel isolates retrieved from oceanic and coastal marine ecosystems. Protist. 2004;155:193–214.
CAS  PubMed  Article  Google Scholar 

93.
Lemieux C, Otis C, Turmel M. Six newly sequenced chloroplast genomes from prasinophyte green algae provide insights into the relationships among prasinophyte lineages and the diversity of streamlined genome architecture in picoplanktonic species. BMC Genom. 2014;15:857.
Article  CAS  Google Scholar  More

Humerus
Left humeri display a non-stationary (i.e., spatially variable in intensity) distribution of tooth marks, with inhomogeneous intensity and with a clustering trend, especially on the proximal end. Tooth marks cluster on the proximal epiphyses, both on the tubercles as well as the articular surface area and proximal metadiaphyses. They also occur in the vicinity of the deltoid crest. Shafts present more abundant modifications on the caudal and medial sides. The cranial and lateral distal shafts show very few tooth marks in comparison. This distribution shows a connection between tooth mark occurrence and areas of muscle and ligament insertions. Tooth marks were probably created during defleshing and limb detachment from the trunk. They are most abundant on the neck junction between the articular head and the proximal metadiaphysis (Fig. 2).
Figure 2

Examples of three-dimensional tooth mark distribution from the lion-consumed carcass sample on each of the four long bones. Distribution of marks is shown on bilateral representation.

Full size image

Both, the K function and the pair-correlation function indicate an overall trend of clustering. This is nuanced by the other functions. The near-neighbour G function shows a slight clustering trend in short distances and a general asymptotic trend of dispersal in longer distances. The empty-space F function suggests a trend towards clustering within an overall CSR pattern (Fig. 3).
Figure 3

Three-dimensional plot of the distribution of tooth marks on the left humerus. (A) K-function plot. (B) G near-neighbour function plot. (C) F empty space function. (C) Pair-correlation function. The F function suggests a pattern non-differentiable from CSR. The other three functions suggest a mild clustering trend in short distances. Key to (A,B,D) Dotted red line shows the Poisson Complete Spatial Random (CSR) process and the gray band shows its confidence envelope. Black line shows the point process of the target sample (here, tooth marks on humerus). When above the CSR Poisson process, it indicates a clustering trend. When below, it indicates a regular scattering trend. The interpretation is reverse for (C) (F empty space function). Same interpretation applies to equivalent figures in the Supplementary Information.

Full size image

Right humeri show a similar tendency of mark clustering around the neck under the articular surface with both tubercles impacted. Marks on the medial shaft are slightly more abundant than on the lateral shaft, while the latter shows higher concentrations around the deltoid crest. Interestingly, marks on the shaft cluster on the proximal and distal portions and the mid-shaft is mostly devoid of marks, regardless of orientation. The cranial side, especially the shaft, is again the least impacted by lions. All the functions show a moderate tendency to clustering; so much so in the K and pair-correlation functions because the latter is a modified version of the former (by using rings within the distance radius). The G function shows a very slight clustering trend in short distances, which in the F function is barely outside the CSR envelope (Fig. S4).
In sum, the slight clustering trends in both sides indicate a redundant pattern of tooth mark location. This shows that mark distribution in humeri is not random, since it is repeated across all the carcasses studied.
Femur
Left femora do not show a more widespread distribution of tooth marks than documented in both humeri. Most tooth marks also occur on the proximal half of the element. Most distal tooth marks appear concentrated on the epiphyses. They occur mostly on the medial condyle (on its medial facet) and on the medial portion of the trochlea. Marks on the proximal end occur on the trochanters and also on the spiral line of the neck. The lateral sides of the shaft are the least modified, followed by the caudal distal shaft. Tooth marks on the caudal shaft occur on both sides of the line aspera. As was the case with the humerus, a large portion of marks appear at or near muscle insertion areas. All the functions show a slight clustering trend in short distances and a CSR pattern in longer distances (Fig S5).
Right femora appear substantially more toothmarked than the left ones. Again, the proximal and distal ends exhibit the highest amount of marks. Both trochanters and the proximal metadiaphysis contain large numbers of modifications. Marks on the distal epiphysis occur both on the medial facet of the throclea and on both condyles. Marks on the caudal shaft, along the linea aspera, are more abundant than on the cranial shaft. All functions coincide in finding a moderate clustering trend, which indicates that BSM are not following a CSR pattern (Fig. S6).
As was the case for humeri, the non-random and moderately clustered pattern shows that there are locations, mostly coinciding with tendon and muscle insertions, that are more prone to be impacted by lions during carcass consumption than others.
Radius-ulna
Radii from carcasses consumed by lions are generally left unmodified12,26,27. Most of the damage concentrates on the olecranon of the ulna (Fig. S7). Only a few tooth marks have been documented scattered on the proximal metadiaphysis, some under the articular facet of the lateral epiphysis. The rest occur mostly in the form of isolated marks, without any specific preference for clustering or side. The left radius shows this distribution. Marks outside the ulna are very few and occur on the cranial and lateral sides of the proximal metadiaphysis, in proximity to the articular facet. Scattered marks can be observed on the distal end. In contrast with the stylopodials, the left radius-ulna shows more intense clustering of tooth marks, as denoted by the K,G, F and pair-correlation functions (Fig. S4). This may be the effect of the intense damage on the olecranon.
The right radius appears also very slightly toothmarked, despite the large number of carcasses involved. Most tooth marks concentrate on the ulnar olecranon, with very few scattered along the ulnar shaft and even less so on the radial shaft. The few tooth marks documented on the shaft appear on the uppermost cranial shaft and a couple on the lower caudal shaft. As was the case with the left radius, the second-order functions indicate a clear clustering of tooth marks in slightly longer distances than documented in the stylopods (Fig. S8).
In sum, marks in radii are few and mostly clustered on the ulna. Those on the radial shaft are scattered but also seem to be in connection with damage on the proximal end imparted during defleshing by lions.
Tibia
The left tibia shows a concentration of tooth marks on the proximal end, more specifically, on the epiphysis and, especially, on the crest. Marks on the shaft are not common and they cluster mostly on the lateral and medial sides and on the lateral portion of the caudal side. Marks in the lower half of the shaft are uncommon regardless of orientation (Fig. S9). This element exhibits the lowest frequency of marks of the whole long bone set. The second-order functions indicate a very minor clustering trend in short distances, probably caused by redundancy in damage in the proximal portion of the element, but most of the shaft, where the few scattered marks occur, seems very similar to a Poisson process. This suggests that damage to the tibia (with the exception of the crest and proximal end) is more stochastic than on the other elements.
Right tibiae are only slightly more toothmarked than the left tibiae. Given its overall greater length than other long bones, its low toothmarking frequencies renders them the least impacted elements in number of tooth marks. Most marks cluster on the proximal end, more specifically, on the tibial crest. The lateral side is more damaged than the other sides. In the whole collection, only one tooth mark was found in the distal half of the shaft (Fig. 2). Again, most of the damage on the caudal side was concentrated on the proximal lateral side, coinciding with the more intensive damage on the lateral portion of the cranial side. The second-order functions suggest also a very minor clustering trend, slightly more marked than on left tibiae, probably because all tooth marks documented concentrate on the proximal half of the element (Fig. S10).
In sum, tibiae show some of the least intense point processes resulting from toothmarking by lions on long bones. Marks occurring on the shaft are usually isolated and more random than on other elements, where they are more spatially recurrent.
Bilateral element comparison
Left and right humeri display a similar pattern in the location of most damage as the three-dimensional coordinates of the PCA show (Fig. 4). This is reinforced by the bivariate wavelet analysis, which shows that both sides of humeri show a strong correlation ( > 0.8) in the location of most tooth marks in specific locations (Fig. 5). Both humeri display high frequencies of tooth marks (remember, the lower the frequency, the higher the scale) and a clear clustering on the proximal epiphysis and proximal metadiaphysis, as well as on the distal shaft. Most of the mid-shaft shows almost no tooth marks and when they do, they occur in very low frequencies. High frequency marks have only been documented on the proximal epiphyseal portion (Fig. 5). The frequency distribution also shows that the medial and caudal sides bear more marks than the lateral and cranial sides. Most cranial marks are concentrated in the articular surface, tubercles and metadiaphyseal portion of the proximal end.
Figure 4

Principal component analysis (PCA) of each of the four long bones (humerus, femur, radius-ulna and tibia) according to side (left–right) showing point distribution according to components generated by compressing the three-dimensional coordinates. A 95% confidence ellipse per side shows variation and similarity of toothmark patterns in each of the bones. Percentages shown are for the first and second component respectively.

Full size image

Figure 5

Bivariate wavelet coherence plot showing the correlation of most tooth mark damage on the proximal and distal sections of left and right humeri in low frequencies. Arrows indicate that in these two high-correlation areas, both humeral sides are in phase (i.e., the covary together in the same direction). In the distal area, the right humerus is leading (arrows pointing to the right-down or left-up) and in the proximal area, the left humerus leads (arrows pointing to the right-up or left-down). Binning of histograms is described in Table 3. (A) frequency of marks from distal end (left) to proximal (right) end; (B) frequency of marks from lateral (left) to medial (right), and, (C) frequency of marks on caudal (left) to cranial (right).

Full size image

Table 3 Binning of histograms according to bone length.
Full size table

Left and right femora also display a similar toothmarking pattern (Fig. 4). Both 95% confidence PCA ellipses overlap in most of their areas. The wavelet coherence analysis shows that both sides display a high correlation ( > 0.8) in toothmarking on proximal and distal ends as well as on the shaft when the frequency of marks is low or moderate. Most marks occur on the proximal portion of the element, with a higher impact on the cranial side and more medial for the left femur and more lateral for the right one (Fig. 6). Femoral mid-shafts, thus, appear more highly toothmarked than humeral shafts. Interestingly, the wavelet analysis also shows that when modifications are abundant, there is correspondence between left and right sides only at the distal end. This seems to respond to bone and muscle insertions and ways in which lions deflesh carcasses at this part of the limb. A moderate correlation ( > 0.6) between both sides of the element can be found at the level of the proximal articular neck (metadiaphysis) and surrounding the trochanter section (see yellow islands at the level of the 20th-23rd bins in Fig. 6).
Figure 6

Bivariate wavelet coherence plot showing the correlation of tooth mark damage on the proximal and distal sections of left and right femora in moderate frequencies. Arrows indicate that in these two high-correlation areas, both femoral sides are in phase (i.e., the covary together in the same direction). The right femur is leading (arrows pointing to the right-down or left-up). Binning of histograms is described in Table 3. (A) frequency of marks from distal end (left) to proximal (right) end; (B) frequency of marks from lateral (left) to medial (right), and, (C) frequency of marks on caudal (left) to cranial (right).

Full size image

As was the case of the upper limb bones, radii-ulnae also exhibit a localized tooth mark pattern. The 95% confidence PCA ellipses overlap for both sides is more intense even than with the stylopodials. The only points falling outside the confidence ellipse are those that appear in the form of single marks and are caused stochastically. The wavelet coherence analysis indicates a strong pattern between both sides, with marks clustering in the proximal epiphysis and strong correlation in the exhibition of low-impact modifications (i.e., few isolates marks) in most of the shaft (high scale = low frequency). There is a high frequency of modifications on the proximal end (see black line sloping upwards in Fig. 7), which decreases as we go down the shaft. The low frequency is maintained throughout the length of the shaft. Only because a few more marks have been documented on the distal and proximal ends, do we see a lower scale (i.e., higher frequency) at the beginning and end of the plot. The high correlation spread along the element shaft indicates that both the right and left radii-ulnae display virtually the same modification pattern.
Figure 7

Bivariate wavelet coherence plot showing the correlation of tooth mark damage on the proximal and distal sections of left and right radius-ulna in moderate to high frequencies. Arrows indicate that in these two high-correlation areas, both femoral sides are in phase (i.e., the covary together in the same direction). The right radius-ulna is leading (arrows pointing to the right-down or left-up). Binning of histograms is described in Table 3. (A) frequency of marks from distal end (left) to proximal (right) end; (B) frequency of marks from lateral (left) to medial (right), and, (C) frequency of marks on caudal(left) to cranial (right).

Full size image

Tibiae also show similar tooth-marking patterns when comparing right and left sides of the skeleton. A PCA shows that a 95% confidence ellipse of samples from both sides overlap in most of their areas (Fig. 4). However, it should be remarked that there is more coordinate variation (i.e., variation in distribution) of tooth marks in tibiae compared to the other long bones. The reason may be double. On the one hand, the tibia exhibits the longest length dimensions of the appendicular skeleton. On the other side, the occurrence of tooth marks outside the area surrounding the tibial crest is commonly in the form of isolated marks that are more prone to occur randomly during defleshing because no muscle insertions occur on the cranial aspect of the element. Only in the proximal caudal side are tooth marks more prone to cluster because of the muscle insertions on that side. A wavelet coherence analysis shows that tibiae show a low density of modifications, similar to radii-ulnae but over a more widespread area. This creates a situation of high correlation between the left and right sides in the location of the few scattered marks (Fig. 8). The correlation is also similar in the proximal and distal ends when modifications are more clustered. Overall, the lack of intensive (i.e., abundant clustering) modifications on the shaft, makes both tibial sides to lack a pattern, with the exception of the lateral and caudal proximal shafts. This moderate clustering there creates the small peninsula between bins 5 and 11 of Fig. 8.
Figure 8

Bivariate wavelet coherence plot showing the correlation of tooth mark damage on the proximal and distal sections of left and right tibiae in moderate to high frequencies. Notice different location of proximal and distal ends compared to the other elements. Arrows indicate that in these two high-correlation areas, both tibial sides are in phase (i.e., the covary together in the same direction). Binning of histograms is described in Table 3. (A) frequency of marks from distal end (right) to proximal (left) end; (B) frequency of marks from lateral (right) to medial (left), and, (C) frequency of marks on caudal(left) and cranial (right).

Full size image

In summary, the humeri, femora and radii-ulnae exhibit strong patterning on how lions modify them after consumption, as reflected in tooth mark distribution on both sides of the same elements. The tibiae display a more variable pattern, which overall is reflected on fewer modifications, especially along the shaft. Given the commonly isolated nature of most marks created along the shaft, these respond more to stochastic processes and reflect higher variability than in the other elements. Exceptions to this observation are found in BSM observed on the tibial crest and proximal caudal-lateral portions of the shaft.
Multi-element comparison
The information contained in the three-dimensional coordinates of the toothmark pattern documented on each of the elements, when approached through the holistic consideration of the mean values of their global interrelation (as documented through the second-order functions), provides identity information (i.e., element-specific identification) for each of the bones analyzed. On a different scale, this could be applied to individual assemblages instead of individual elements as done here. In the comparison among the different elements and their sides, the way marks were distributed in each of their respective point processes (considering their intensity and distances per element) contained sufficient information to differentiate four different clusters corresponding to the four different elements (Table 2, Fig. 9). Within each element set, both sides were contained within the same node. This is of utmost interest, because in the variables used for this analysis, it is the patterns and not the raw coordinates of marks on each element that were used. This enabled the relativization of the actual location of marks on the different long bone elements and only the emergent properties of the mark assemblage in each of them (understood as individual point process) was considered. Thus, multi-element comparison was possible and different bones were successfully differentiated (Fig. 9).
Figure 9

Hierarchical clustering of the selected variables from the second-order functions, intensity, and nearest-neighbour distance. A phylogenetic dendrogram was used. Four groups were identified (different colors) corresponding to each of the four elements analyzed. Key: lHum (left humerus); rHum (right humerus); lFem (left femur); rFem (right femur); lRad (left radius-ulna); rRad (right radius-ulna); lTib (left tibia); rTib (right tibia).

Full size image More

Alcala N, Goudet J, Vuilleumier S (2014) On the transition of genetic differentiation from isolation to panmixia: what we can learn from GST and D. Theor Popul Biol 93:75–84
PubMed  Article  Google Scholar 

Barton NH (1996) Natural selection and random genetic drift as causes of evolution on islands. Philos Trans R Soc Lond B Biol Sci 351:785–795
CAS  PubMed  Article  Google Scholar 

Bellemain E, Ricklefs RE (2008) Are islands the end of the colonization road? Trends Ecol Evol 23:461–468
PubMed  Article  Google Scholar 

Benjamini Y, Hochberg Y (1995) Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc Ser B-Methodol 57:289–300
Google Scholar 

Bolger A, Lohse M, Usadel B (2014) Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics 30:2114–2120
CAS  PubMed  PubMed Central  Article  Google Scholar 

Bouckaert R (2010) DensiTree: making sense of sets of phylogenetic trees. Bioinformatics 26:1372–1373
CAS  PubMed  Article  Google Scholar 

Bouckaert R, Heled J, Kühnert D, Vaughan T, Wu C-H, Xie D et al. (2014) BEAST 2: a software platform for Bayesian evolutionary analysis. PLoS Comput Biol 10:e1003537
PubMed  PubMed Central  Article  CAS  Google Scholar 

Bryant D, Moulton V (2004) Neighbor-net: an agglomerative method for the construction of phylogenetic networks. Mol Biol Evol 21:255–265
CAS  PubMed  Article  Google Scholar 

Bryant D, Bouckaert R, Felsenstein J, Rosenberg NA, RoyChoudhury A (2012) Inferring species trees directly from biallelic genetic markers: bypassing gene trees in a full coalescent analysis. Mol Biol Evol 29:1917–1932
CAS  PubMed  PubMed Central  Article  Google Scholar 

Catchen J, Hohenlohe PA, Bassham S, Amores A, Cresko WA (2013) Stacks: an analysis tool set for population genomics. Mol Ecol 22:3124–3140
PubMed  PubMed Central  Article  Google Scholar 

Chapin III FS, Matson PA, Vitousek PM (2002) Principles of terrestrial ecosystem ecology. Springer, New York
Google Scholar 

Clark PU, Dyke AS, Shakun JD, Carlson AE, Clark J, Wohlfarth B et al. (2009) The last glacial maximum. Science 325:710–714
CAS  PubMed  Article  Google Scholar 

DeChaine EG, Martin AP (2005) Marked genetic divergence among sky island populations of Sedum lanceolatum (Crassulaceae) in the Rocky Mountains. Am J Bot 92:477–486
CAS  PubMed  Article  Google Scholar 

Drummond AJ, Rambaut A (2007) BEAST: Bayesian evolutionary analysis by sampling trees. BMC Evol Biol 7:214
PubMed  PubMed Central  Article  CAS  Google Scholar 

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

Foote AD, Morin PA (2016) Genome-wide SNP data suggest complex ancestry of sympatric North Pacific killer whale ecotypes. Heredity 117:316–325
CAS  PubMed  PubMed Central  Article  Google Scholar 

Frankham R (1997) Do island populations have less genetic variation than mainland populations? Heredity 78:311–327
PubMed  Article  Google Scholar 

Frichot E, François O (2015) LEA: an R package for landscape and ecological association studies. Methods Ecol Evol 6:925–929
Article  Google Scholar 

Frichot E, Mathieu F, Trouillon T, Bouchard G, François O (2014) Fast and efficient estimation of individual ancestry coefficients. Genetics 196:973–983
PubMed  PubMed Central  Article  Google Scholar 

Frichot E, Schoville SD, Bouchard G, François O (2013) Testing for associations between loci and environmental gradients using latent factor mixed models. Mol Biol Evol 30:1687–1699
CAS  PubMed  PubMed Central  Article  Google Scholar 

Frichot E, Schoville SD, De Villemereuil P, Gaggiotti OE, François O (2015) Detecting adaptive evolution based on association with ecological gradients: orientation matters! Heredity 115:22–28
CAS  PubMed  PubMed Central  Article  Google Scholar 

Funk WC, Lovich RE, Hohenlohe PA, Hofman CA, Morrison SA, Sillett TS et al. (2016) Adaptive divergence despite strong genetic drift: genomic analysis of the evolutionary mechanisms causing genetic differentiation in the island fox (Urocyon littoralis). Mol Ecol 25:2176–2194
CAS  PubMed  PubMed Central  Article  Google Scholar 

Garot E, Joët T, Combes M, Severac D, Lashermes P (2019) Plant population dynamics on oceanic islands during the Late Quaternary climate changes: genetic evidence from a tree species (Coffea mauritiana) in Reunion Island. N Phytol 224:974–986
CAS  Article  Google Scholar 

Gillespie R (2004) Community assembly through adaptive radiation in Hawaiian spiders. Science 303:356–359
CAS  PubMed  Article  Google Scholar 

Hamabata T, Kinoshita G, Kurita K, Cao PL, Ito M, Murata J et al. (2019) Endangered island endemic plants have vulnerable genomes. Commun Biol 2:244
PubMed  PubMed Central  Article  Google Scholar 

Hedrick PW (2005) A standardized genetic differentiation measure. Evolution 59:1633–1638
CAS  PubMed  Article  Google Scholar 

Huson DH, Bryant D (2006) Application of phylogenetic networks in evolutionary studies. Mol Biol Evol 23:254–267
CAS  PubMed  Article  Google Scholar 

Izuno A, Kitayama K, Onoda Y, Tsujii Y, Hatakeyama M, Nagano AJ et al. (2017) The population genomic signature of environmental association and gene flow in an ecologically divergent tree species Metrosideros polymorpha (Myrtaceae). Mol Ecol 26:1515–1532
CAS  PubMed  Article  Google Scholar 

James JE, Lanfear R, Eyre-Walker A (2016) Molecular evolutionary consequences of island colonization. Genome Biol Evol 8:1876–1888
CAS  PubMed  PubMed Central  Article  Google Scholar 

Juan C, Emerson BC, Oromı́ P, Hewitt GM (2000) Colonization and diversification: towards a phylogeographic synthesis for the Canary Islands. Trends Ecol Evol 15:104–109
CAS  PubMed  Article  Google Scholar 

Kamijo T, Hashiba K (2003) Island ecosystem and vegetation dynamics before and after the 2000-year eruption on Miyake-jima Island, Japan, with implications for conservation of the island’s ecosystem. Glob Environ Res 7:69–78
Google Scholar 

Kier G, Kreft H, Lee TM, Jetz W, Ibisch PL, Nowicki C et al. (2009) A global assessment of endemism and species richness across island and mainland regions. Proc Natl Acad Sci USA 106:9322–9327
CAS  PubMed  Article  Google Scholar 

Koch MA, Haubold B, Mitchell-Olds T (2000) Comparative evolutionary analysis of chalcone synthase and alcohol dehydrogenase loci in Arabidopsis, Arabis, and related genera (Brassicaceae). Mol Biol Evol 17:1483–1498
CAS  Article  Google Scholar 

Liu X, Fu YX (2015) Exploring population size changes using SNP frequency spectra. Nat Genet 47:555–559
CAS  PubMed  PubMed Central  Article  Google Scholar 

MacArthur RH, Wilson EO (1963) An equilibrium theory of insular zoogeography. Evolution 17:373–387
Article  Google Scholar 

Maekawa F (1949) Makinoesia and its bearing to Oriental Asiatic flora. J Jpn Bot 24:91–96. in Japanese
Google Scholar 

McGlaughlin ME, Wallace LE, Wheeler GL, Bresowar G, Riley L, Britten NR et al. (2014) Do the island biogeography predictions of MacArthur and Wilson hold when examining genetic diversity on the near mainland California Channel Islands? Examples from endemic Acmispon (Fabaceae). Bot J Linn Soc 174:289–304
Article  Google Scholar 

Ministry of the Environment of Japan (2019) Threatened wildlife of Japan – red list 2019. http://www.env.go.jp/press/106383.html

Murray M, Thompson W (1980) Rapid isolation of high molecular weight plant DNA. Nucleic Acids Res 8:4321–4326
CAS  PubMed  PubMed Central  Article  Google Scholar 

Nakamura K, Denda T, Kokubugata G, Suwa R, Yang TYA, Peng C-I et al. (2010) Phylogeography of Ophiorrhiza japonica (Rubiaceae) in continental islands, the Ryukyu Archipelago, Japan. J Biogeogr 37:1907–1918
Google Scholar 

Nei M, Tajima F, Tateno Y (1983) Accuracy of estimated phylogenetic trees from molecular data. J Mol Evol 19:153–170
CAS  PubMed  Article  Google Scholar 

Ossowski S, Schneeberger K, Lucas-Lledó JI, Warthmann N, Clark RM, Shaw RG et al. (2010) The rate and molecular spectrum of spontaneous mutations in Arabidopsis thaliana. Science 327:92–94
CAS  PubMed  Article  Google Scholar 

Peterson BK, Weber JN, Kay EH, Fisher HS, Hoekstra HE (2012) Double digest RADseq: an inexpensive method for de novo SNP discovery and genotyping in model and non-model species. PLoS ONE 7:e37135
CAS  PubMed  PubMed Central  Article  Google Scholar 

Petit JR, Jouzel J, Raynaud D, Barkov NI, Barnola J-M, Basile I et al. (1999) Climate and atmospheric history of the past 420,000 years from the Vostok ice core, Antarctica. Nature 399:429–436
CAS  Article  Google Scholar 

Purcell S, Neale B, Todd-Brown K, Thomas L, Ferreira MAR, Bender D et al. (2007) PLINK: a tool set for whole-genome association and population-based linkage analyses. Am J Hum Genet 81:559–575
CAS  PubMed  PubMed Central  Article  Google Scholar 

R Development Core Team (2019) R: a language and environment for statistical computing. https://www.R-project.org/

Rellstab C, Gugerli F, Eckert AJ, Hancock AM, Holderegger R (2015) A practical guide to environmental association analysis in landscape genomics. Mol Ecol 24:4348–4370
PubMed  Article  Google Scholar 

Stervander M, Illera JC, Kvist L, Barbosa P, Keehnen NP, Pruisscher P et al. (2015) Disentangling the complex evolutionary history of the Western Palearctic blue tits (Cyanistes spp.) – Phylogenomic analyses suggest radiation by multiple colonization events and subsequent isolation. Mol Ecol 24:2477–2494
CAS  PubMed  Article  Google Scholar 

Stuessy TF, Takayama K, López-Sepúlveda P, Crawford DJ (2014) Interpretation of patterns of genetic variation in endemic plant species of oceanic islands. Bot J Linn Soc 174:276–288
PubMed  Article  Google Scholar 

Sundqvist L, Keenan K, Zackrisson M, Prodöhl P, Kleinhans D (2016) Directional genetic differentiation and relative migration. Ecol Evol 6:3461–3475
PubMed  PubMed Central  Article  Google Scholar 

Taira A, Saito S, Aoike K, Morita S, Tokuyama H, Suyehiro K et al. (1998) Nature and growth rate of the northern Izu-Bonin (Ogasawara) arc crust and their implications for continental crust formation. Isl Arc 7:395–407
CAS  Article  Google Scholar 

Takahashi H (1971) Fossa Magna element plants. Res Reports Kanagawa Prefect Museum. Nat Hist 2:2–59. in Japanese
Google Scholar 

Takahashi H, Katsuyama T (1992) Natural hybrids between Rhododendron tsusiophyllum and R. kaempferi var. macrogemma (Ericaceae). Bull Kanagawa Prefect. Museum 21:59–71. in Japanese
Google Scholar 

Takezaki N, Nei M, Tamura K (2010) POPTREE2: software for constructing population trees from allele frequency data and computing other population statistics with Windows interface. Mol Biol Evol 27:747–752
CAS  PubMed  Article  Google Scholar 

Tanaka N (1999) Plant communities in Mt. Tenjo, Koudzu Island, Tokyo. Actinia 12:147–158. in Japanese
Google Scholar 

Vaxevanidou Z, González-Martínez SC, Climent J, Gil L (2006) Tree populations bordering on extinction: a case study in the endemic Canary Island pine. Biol Conserv 129:451–460
Article  Google Scholar 

Velo-Antón G, Zamudio K, Cordero-Rivera A (2012) Genetic drift and rapid evolution of viviparity in insular fire salamanders (Salamandra salamandra). Heredity 108:410–418
PubMed  Article  Google Scholar 

Wagner DB, Furnier GR, Saghai-Maroof MA, Williams SM, Dancik BP, Allard RW (1987) Chloroplast DNA polymorphisms in lodgepole and jack pines and their hybrids. Proc Natl Acad Sci USA 84:2097–2100
CAS  PubMed  Article  Google Scholar 

Warren BH, Simberloff D, Ricklefs RE, Aguilée R, Condamine FL, Gravel D et al. (2015) Islands as model systems in ecology and evolution: prospects fifty years after MacArthur-Wilson. Ecol Lett 18:200–217
PubMed  Article  Google Scholar 

Weigelt P, Jetz W, Kreft H (2013) Bioclimatic and physical characterization of the world’s islands. Proc Natl Acad Sci USA 110:15307–15312
CAS  PubMed  Article  Google Scholar 

Yamada T, Maki M (2012) Impact of geographical isolation on genetic differentiation in insular and mainland populations of Weigela coraeensis (Caprifoliaceae) on Honshu and the Izu Islands. J Biogeogr 39:901–917
Article  Google Scholar 

Yamamoto S, Kikuchi T, Yamagiwa Y, Handa T (2017) Genetic diversity of Lilium auratum var. platyphyllum endemic to the Izu archipelago and its relationship to a nearby population of L. auratum var. auratum by morphological and SSR analysis. Hortic J 86:379–388
Article  Google Scholar 

Yoichi W, Minamitani T, Oh S-H, Nagano AJ, Abe H, Yukawa T (2019) New taxa of Rhododendron tschonoskii alliance (Ericaceae) from East Asia. PhytoKeys 134:97–114
PubMed  PubMed Central  Article  Google Scholar  More

1.
Mulec, J. Phototrophs in caves. In Cave Ecology (eds Moldovan, O. T. et al.) 91–106 (Springer, Cham, 2018). https://doi.org/10.1007/978-3-319-98852-8_6
Google Scholar 
2.
Culver, D. C. & Pipan, T. The Biology of Caves and Other Subterranean Habitats (Oxford University Press Inc., New York, 2009).
Google Scholar 

3.
Engel, A. S. Chemoautotrophy. In Encyclopedia of caves 2nd edn (eds White, W. B. & Culver, D. C.) 125–134 (Elsevier, Amsterdam, 2012).
Google Scholar 

4.
Kinkle, B. K. & Kane, T. C. Chemolithotrophic microorganisms and their potential role in subsurface environments. In Ecosystems of the World 30 Subterranean Ecosystems (eds Wilkens, H. et al.) 309–319 (Elsevier, Amsterdam, 2000).
Google Scholar 

5.
Sarbu, S. M. Movile cave: A chemoautotrophically based groundwater ecosystem. In Ecosystems of the World 30 Subterranean Ecosystems (eds Wilkens, H. et al.) 319–343 (Elsevier, Amsterdam, 2001).
Google Scholar 

6.
Simon, K. S., Pipan, T. & Culver, D. C. A conceptual model of the flow and distribution of organic carbon in caves. J. Cave Karst Stud. 69, 279–284 (2007).
CAS  Google Scholar 

7.
Camassa, M. M. Food resources. In Encyclopaedia of Caves and Karst Science (ed. Gunn, J.) 755–760 (Fitzroy Dearborn, London, 2004).
Google Scholar 

8.
Poulson, T. L. & Lavoie, K. H. (The trophic basis of subsurface ecosystems. In Ecosystems of the World 30 Subterranean Ecosystems (eds Wilkens, H. et al.) 323–334 (Elsevier, Amsterdam, 2000).
Google Scholar 

9.
Gibert, J. & Deharveng, L. Subterranean ecosystems: A truncated functional biodiversity. Bioscience 52(6), 473–481. https://doi.org/10.1641/0006-3568(2002)052[0473:SEATFB]2.0 (2002).
Article  Google Scholar 

10.
Chen, B. & Wise, D. H. Bottom-up limitation of predaceous arthropods in a detritus-based terrestrial food web. Ecology 80(3), 761–772. https://doi.org/10.2307/177015 (1999).
Article  Google Scholar 

11.
Venarsky, M. P. & Huntsman, B. M. Food webs in caves. In Cave Ecology (eds Moldovan, O. T. et al.) 309–331 (Springer, Cham, 2018). https://doi.org/10.1007/978-3-319-98852-8_14
Google Scholar 

12.
Gnaspini, P. Guano communities. In Encyclopedia of caves 2nd edn (eds White, W. B. & Culver, D. C.) 357–364 (Elsevier, Amsterdam, 2012).
Google Scholar 

13.
Ipsen, A. The Segeberger Höhle—A phylogenetically young cave ecosystem in northern Germany. In Ecosystems of the World 30. Subterranean Ecosystems (eds Wilkens, H. et al.) 569–579 (Elsevier, Amsterdam, 2000).
Google Scholar 

14.
Stone, F. D., Howarth, F. G., Hoch, H. & Asche, M. Root communities in lava tubes. In Encyclopedia of Caves 2nd edn (eds White, W. B. & Culver, D. C.) 658–664 (Elsevier, Amsterdam, 2012).
Google Scholar 

15.
Mammola, S., Piano, E. & Isaia, M. Step back! Niche dynamics in cave-dwelling predators. Acta Oecol. 75, 35–42. https://doi.org/10.1016/j.actao.2016.06.011 (2016).
ADS  Article  Google Scholar 

16.
Mammola, S. & Isaia, M. Cave communities and species interactions. In Cave Ecology (eds Moldovan, O. T. et al.) 255–269 (Springer, Cham, 2018). https://doi.org/10.1007/978-3-319-98852-8_11
Google Scholar 

17.
Scheu, S. & Setälä, H. Multitrophic interactions in decomposer food webs. In Multitrophic Interactions in Terrestrial Systems (eds Tscharntke, T. & Hawkins, B. A.) 223–264 (Cambridge, Cambridge University Press, 2001).
Google Scholar 

18.
Wood, P. J. Subterranean ecology. In Encyclopaedia of Caves and Karst Science (ed. Gunn, J.) 1514–1519 (Fitzroy Dearborn, London, 2004).
Google Scholar 

19.
Pekár, S., García, L. F. & Viera, C. Trophic niche and trophic adaptations of prey specialised spiders of the Neotropics: A guide. In Behavioural Ecology of Neotropical Spiders (eds Viera, C. & Gonzaga, M. O.) 247–274 (Springer, Cham, 2017).
Google Scholar 

20.
Pohlman, J. W., Iliffe, T. M. & Cifuentes, L. A. A stable isotope study of organic cycling and the ecology of an anchialine cave ecosystem. Mar. Ecol. Prog. Ser. 155, 17–27 (1997).
ADS  CAS  Article  Google Scholar 

21.
Pohlman, J. W., Cifuentes, L. A. & Iliffe, T. M. Food web dynamics and biogeochemistry of anchialine caves: A stable isotope approach. In Ecosystems of the World 30 Subterranean Ecosystems (eds Wilkens, H. et al.) 345–357 (Elsevier, Amsterdam, 2000).
Google Scholar 

22.
Sarbu, S. M., Galdenzi, S., Menichetti, M. & Gentile, G. Geology and biology of the Frasassi caves in Central Italy: An ecological multi-disciplinary study of a hypogenic underground karst system. In Ecosystems of the World 30 Subterranean Ecosystems (eds Wilkens, H. et al.) 359–378 (Elsevier, Amsterdam, 2000).
Google Scholar 

23.
Eitzinger, B., Micic, A., Körner, M., Traugott, M. & Scheu, S. Unveiling soil food web links: New PCR assays for detection of prey DNA in the gut of soil arthropod predators. Soil Biol. Biochem. 57, 943–945. https://doi.org/10.1016/j.soilbio.2012.09.001 (2013).
CAS  Article  Google Scholar 

24.
Juen, A. & Traugott, M. Revealing species-specific trophic links in soil food webs: Molecular identification of scarab predators. Mol. Ecol. 16, 1545–1557. https://doi.org/10.1111/j.1365-294X.2007.03238.x (2007).
CAS  Article  PubMed  Google Scholar 

25.
King, R. A., Read, D. S., Traugott, M. & Symondson, W. O. C. Molecular analysis of predation: A review of best practice for DNA-based approaches. Mol. Ecol. 17, 947–963. https://doi.org/10.1111/j.1365-294X.2007.03613.x (2008).
CAS  Article  PubMed  Google Scholar 

26.
Symondson, W. O. C. Molecular identification of prey in predator diets. Mol. Ecol. 11(4), 627–641. https://doi.org/10.1046/j.1365-294x.2002.01471.x (2002).
CAS  Article  PubMed  Google Scholar 

27.
Traugott, M., Kamenova, S., Ruess, L., Seeber, J. & Plantegenest, M. Empirically characterising trophic networks: What emerging DNA-based methods, stable isotope and fatty acid analyses can offer. Adv. Ecol. Res. 49, 177–224. https://doi.org/10.1016/B978-0-12-420002-9.00003-2 (2013).
Article  Google Scholar 

28.
Kováč, Ľ. et al. Terrestrial arthropods of the Domica Cave system and the Ardovská Cave (Slovak Karst): Principal microhabitats and diversity. In Contributions to Soil Zoology in Central Europe I (eds Tajovský, K. et al.) 61–70 (ISB AS CR, České Budějovice, 2005).
Google Scholar 

29.
Kováč, Ľ. et al. The cave biota of Slovakia. Speleologia Slovaca 5. (Liptovský Mikuláš, State Nature Conservancy SR, Slovak Caves Administration, 2014). https://doi.org/10.13140/2.1.3473.0569

30.
Kováč, Ľ, Parimuchová, A. & Miklisová, D. Distributional patterns of cave Collembola (Hexapoda) in association with habitat conditions, geography and subterranean refugia in the Western Carpathians. Biol. J. Linn. Soc. Lond. 119(3), 571–592. https://doi.org/10.1111/bij.12555 (2016).
Article  Google Scholar 

31.
Smrž, J., Kováč, Ľ, Mikeš, J. & Lukešová, A. Microwhip scorpions (Palpigradi) feed on heterotrophic Cyanobacteria in Slovak caves: A curiosity among Arachnida. PLoS ONE 8(10), e75989. https://doi.org/10.1371/journal.pone.0075989 (2013).
ADS  CAS  Article  PubMed  PubMed Central  Google Scholar 

32.
Pekár, S., Coddington, J. A. & Blackledge, T. Evolution of stenophagy in spiders (Araneae): Evidence based on the comparative analysis of spider diets. Evolution 66(3), 776–806. https://doi.org/10.1111/j.1558-5646.2011.01471.x (2012).
Article  Google Scholar 

33.
Alderweireldt, M. Prey selection and prey capture strategies of linyphiid spiders in highinput agricultural fields. Bull. Br. Arachnol. Soc. 9, 300–308 (1994).
Google Scholar 

34.
Lukić, M., Collembola in caves. Croatian Biospeleological Society, DVD, 10.25 min (2012).

35.
Roewer, C. F. Palpigradi. In Klassen und Ordnungen des Tierreichs 5: Arthropoda IV: Arachnoidea (ed. Bronns, H. G.) 640–707 (Akademische Verlagsgesellschaft MBH, Leipzig, 1932).
Google Scholar 

36.
van der Hammen, L. Comparative studies in Chelicerata II. Epimerata (Palpigradi and Actinotrichida). Zool. Verh. 196, 3–70 (1982).
Google Scholar 

37.
Wheeler, W. M. A singular arachnid Koenenia mirabilis (Grassi) occurring in Texas. Am. Nat. 34, 837–850 (1900).
Article  Google Scholar 

38.
Harwood, J. D., Phillips, S. W., Sunderland, K. D. & Symondson, W. O. C. Secondary predation: quantification of food chain errors in an aphid–spider–carabid system using monoclonal antibodies. Mol. Ecol. 10(8), 2049–2057. https://doi.org/10.1046/j.0962-1083.2001.01349.x (2001).
CAS  Article  PubMed  Google Scholar 

39.
Szafranek, P., Lewandowski, M. & Kozak, M. Prey preference and life tables of the predatory mite Parasitus bituberosus (Acari: Parasitidae) when offered various prey combinations. Exp. Appl. Acarol. 61(1), 53–67. https://doi.org/10.1007/s10493-013-9701-y (2013).
Article  PubMed  PubMed Central  Google Scholar 

40.
Al-Amidi, A. H. K. & Downes, M. J. Parasitus bituberosus (Acari: Parasitidae), a possible agent for biological control of Heteropeza pygmaea (Diptera: Cecidomyiidae) in mushroom compost. Exp. Appl. Acarol. 8(1–2), 13–25 (1990).
Article  Google Scholar 

41.
Adams, B. J. & Nguyen, K. B. Nematode parasites of insects. In Encyclopedia of Entomology (ed. Capinera, J. L.) 2577–2584 (Springer, Cham, 2008).
Google Scholar 

42.
Cokendolpher, J. C. Pathogens and parasites of opiliones (arthropoda: arachnida). J. Arachnol. 21(2), 120–146 (1993).
Google Scholar 

43.
Kruse, P. D., Toft, S. & Sunderland, K. D. Temperature and prey capture: Opposite relationships in two predator taxa. Ecol. Entomol. 33(2), 305–312. https://doi.org/10.1111/j.1365-2311.2007.00978.x (2008).
Article  Google Scholar 

44.
Krooss, S. & Schaefer, M. How predacious are predators? A study on Ocypus similis, a rove beetle of cereal fields. Ann. Appl. Biol. 133(1), 1–16. https://doi.org/10.1111/j.1744-7348.1998.tb05797.x (1998).
Article  Google Scholar 

45.
Waldbauer, G. P. & Friedman, S. Self-selection of optimal diets by insects. Annu. Rev. Entomol. 36(1), 43–63. https://doi.org/10.1146/annurev.en.36.010191.000355 (1991).
Article  Google Scholar 

46.
Mayntz, D. & Toft, S. Nutrient composition of the prey’s diet affects growth and survivorship of a generalist predator. Oecologia 127, 207–213. https://doi.org/10.1007/s004420000591 (2001).
ADS  Article  PubMed  Google Scholar 

47.
Finke, D. L. & Denno, R. F. Intraguild predation diminished in complex-structured vegetation: implications for prey suppression. Ecology 83, 643–652. https://doi.org/10.2307/3071870 (2002).
Article  Google Scholar 

48.
Staudacher, K. et al. Habitat heterogeneity induces rapid changes in the feeding behaviour of generalist arthropod predators. Funct. Ecol. 32(3), 809–819. https://doi.org/10.1111/1365-2435.13028 (2018).
Article  PubMed  PubMed Central  Google Scholar 

49.
Finke, D. L. & Denno, R. F. Predator diversity and the functioning of ecosystems: the role of intraguild predation in dampening trophic cascades. Ecol. Lett. 8, 1299–1306. https://doi.org/10.1111/j.1461-0248.2005.00832.x (2005).
Article  Google Scholar 

50.
Schausberger, P. & Croft, B. A. Nutritional benefits of intraguild predation and cannibalism among generalist and specialist phytoseiid mites. Ecol. Entomol. 25(4), 473–480. https://doi.org/10.1046/j.1365-2311.2000.00284.x (2000).
Article  Google Scholar 

51.
Schausberger, P. Cannibalism among phytoseiid mites: a review. Exp. Appl. Acarol. 29(3/4), 173–191. https://doi.org/10.1023/a:1025839206394 (2003).
Article  PubMed  Google Scholar 

52.
Elgar, M. A. & Crespi, B. J. Cannibalism: Ecology and Evolution Among Diverse Taxa (Oxford University Press, Oxford, 1992).
Google Scholar 

53.
Polis, G. A. The evolution and dynamics of intraspecific predation. Annu. Rev. Ecol. Syst. 12(1), 225–251. https://doi.org/10.1146/annurev.es.12.110181.001301 (1981).
Article  Google Scholar 

54.
Fagan, W. F. et al. Nitrogen in insects: Implications for trophic complexity and species diversification. Am. Nat. 160(6), 784–802. https://doi.org/10.1086/343879 (2002).
Article  Google Scholar 

55.
Fagan, W. F. & Denno, R. F. Stoichiometry of actual vs. potential predator–prey interactions: Insights into nitrogen limitation for arthropod predators. Ecol. Lett. 7(9), 876–883. https://doi.org/10.1111/j.1461-0248.2004.00641.x (2004).
Article  Google Scholar 

56.
Denno, R. F. & Fagan, W. F. Might nitrogen limitation promote omnivory among carnivorous arthropods?. Ecology 84(10), 2522–2531. https://doi.org/10.1890/02-0370 (2003).
Article  Google Scholar 

57.
Snyder, W. E., Joseph, S. B., Preziosi, R. F. & Moore, A. J. Nutritional benefits of cannibalism for the lady beetle Harmonia axyridis (Coleoptera: Coccinellidae) when prey quality is poor. Environ. Entomol. 29(6), 1173–1179. https://doi.org/10.1603/0046-225x-29.6.1173 (2000).
Article  Google Scholar 

58.
Nováková, A. et al. Feeding sources of invertebrates in the Ardovská Cave and Domica Cave systems: preliminary results. In Contributions to Soil Zoology in Central Europe I (eds Tajovský, K. et al.) 107–112 (ISB AS CR, České Budějovice, 2005).
Google Scholar 

59.
Crossley, D. & Blair, J. M. A high efficiency, “low-technology” Tullgren-type extractor for soil microarthropods. Agric. Ecosyst. Environ. 34, 187–192 (1991).
Article  Google Scholar 

60.
Folmer, O., Black, M., Hoeh, W., Lutz, R. & Vrijenhoek, R. DNA primers for amplification of mitochondrial cytochrome c oxidase subunit I from diverse metazoan invertebrates. Mol. Mar. Biol. Biotechnol. 3(5), 294–299 (1994).
CAS  PubMed  Google Scholar 

61.
de Groot, A. G., Laros, I. & Geisen, S. Molecular identification of soil eukaryotes and focused approaches targeting protist and faunal groups using high-throughput meta-barcoding methods in molecular biology. Methods Mol. Biol. 1399, 125–140. https://doi.org/10.1007/978-1-4939-3369-3_7 (2016).
CAS  Article  Google Scholar 

62.
Aronesty, E. Comparison of sequencing utility programs. Open Bioinform. J. 7(1), 1–8. https://doi.org/10.2174/1875036201307010001 (2013).
MathSciNet  Article  Google Scholar 

63.
Martin, M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet J. 17, 10–12. https://doi.org/10.14806/ej.17.1.200 (2011).
Article  Google Scholar 

64.
Mahé, F., Rognes, T., Quince, C., de Vargas, C. & Dunthorn, M. Swarm: robust and fast clustering method for amplicon-based studies. PeerJ. 2, e593 (2014).
Article  Google Scholar 

65.
Belshaw, R., Lopez-Vaamonde, C., Degerli, N. & Quicke, D. L. J. Paraphyletic taxa and taxonomic chaining: Evaluation the classification of braconine wasps (Hymenoptera: Braconidae) using 28S D2–3 rDNA sequences and morphological characters. Biol. J. Linn. Soc. Lond. 73(4), 411–424. https://doi.org/10.1111/j.1095-8312.2001.tb01370.x (2001).
Article  Google Scholar 

66.
Hurlbert, S. H. The measurement of niche overlap and some relatives. Ecology 59(1), 67–77. https://doi.org/10.2307/1936632 (1978).
Article  Google Scholar 

67.
Novakowski, G. C., Hahn, N. S. & Fugi, R. Diet seasonality and food overlap of the fish assemblage in a pantanal pond. Neotrop. Ichthyol. 6(4), 567–576. https://doi.org/10.1590/S1679-62252008000400004 (2008).
Article  Google Scholar 

68.
Pianka, E. R. The structure of lizard communities. Annu. Rev. Ecol. Syst. 4(1), 53–74. https://doi.org/10.1146/annurev.es.04.110173.000413 (1973).
Article  Google Scholar 

69.
Pekár, S. & Brabec, M. Modern Analysis of Biological Data. Generalized Linear Models in R (MUNI Press, Brno, 2016).
Google Scholar 

70.
R Core Team. R: A Language and Environment for Statistical Computing. (R Foundation for Statistical Computing, Vienna). https://www.R-project.org/ (2017).

71.
Breheny, P. & Burchett, W. Visualization of regression models using visreg. R J. 9, 56–71 (2017).
Article  Google Scholar 

72.
Kučera, B. Krasová morfologie a vývoj Ardovské jeskyně v Jihoslovenském krasu. Československý Kras. 16, 41–56 (1964) ([in Czech]).
Google Scholar  More

Author notes
These authors jointly supervised this work: Elizabeth A. Archie and José M. Montoya.

Affiliations

Department of Biological Sciences, University of Notre Dame, Notre Dame, IN, USA
Johannes R. Björk & Elizabeth A. Archie

Theoretical and Experimental Ecology Station, CNRS-University Paul Sabatier, Moulis, France
Johannes R. Björk & José M. Montoya

Natural History Museum, London, UK
Cristina Díez-Vives

School of Biological Sciences, University of Auckland, Auckland, New Zealand
Carmen Astudillo-García

Authors
Johannes R. Björk

Cristina Díez-Vives

Carmen Astudillo-García

Elizabeth A. Archie

José M. Montoya

Corresponding authors
Correspondence to Johannes R. Björk or Elizabeth A. Archie or José M. Montoya. More