Alveolates (dinoflagellates, ciliates and apicomplexans) and Rhizarians are the most common microbial eukaryotes in temperate Appalachian karst caves.

The purpose of this study was to survey the eukaryotic microbiome of two karst caves in the Valley and Ridge physiographic region of the Appalachian Mountains. Caves are known to harbour eukaryotic microbes but their very low densities and small cell size make them difficult to collect and identify. Microeukaryotes were surveyed using two methodologies, filtering water and submerging glass microscope slides mounted in periphytometers in cave pools. The periphyton sampling yielded 13.5 times more unique amplicon sequence variants (ASVs) than filtered water. The most abundant protist supergroup was Alveolata with large proportions of the ASVs belonging to dinoflagellate, ciliate and apicomplexan clades. The next most abundant were Rhizarians followed by Stramenopiles (diatoms and chrysophytes) and Ameobozoans. Very few of the ASVs, 1.5%, matched curated protist sequences with greater than 99% identity and only 2.5% could be identified from surface plankton samples collected in the same region. The overall composition of the eukaryotic microbiome appears to be a combination of bacterial grazers and parasitic species that could possibly survive underground as well as cells, cysts and spores probably transported from the surface.

Introduction

Caves are extremely oligotrophic environments yet harbour a multitude of cave‐adapted organisms, from charismatic macroscopic fauna to extremophilic prokaryotes (Culver and Pipan, 2009; Engel, 2010; Tomczyk‐Żak and Zielenkiewicz, 2016). These survive without the benefit of sunlight or abundant primary producers so biological energy is limited to dissolved organic material and microbes transported from the surface by water movement (Laiz et al., 1999, Simon et al., 2003, Engel and Northup, 2008) and chemolithoautotrophic cave‐adapted microbes (Northup and Lavoie, 2001; Chen et al., 2009; Ortiz et al., 2014).

Cave‐dwelling prokaryotic microbes have been the focus of numerous surveys, beginning with culturable species and most recently molecular metagenomic barcoding techniques. These have revealed a wide array of prokaryotes including many potential endemic species, suggesting caves may act as refugia for various prokaryotic microbes, some of which are pathogenic or have unusual physiologies (Engel, 2010; Igreja, 2011). Recent metabarcoding surveys focused on the microbial ecology of caves has revealed that only 11.2%–21.4% of the prokaryotic barcodes found on speleothems (cave formations) and 53.8% of the microbes in cave sediments were found in nearby surface soils (Ortiz et al., 2013; Lavoie et al., 2017; Thompson et al., 2019). The use of barcoding has also revealed the movement of microbes from the surface and epikarst into karst caves (Morse et al., 2021). These results suggest that, despite the constant transfer of microbes from the surface, the cave ecosystem, especially speleothems, are inhabited by many cave‐endemic microbes.

The eukaryotic microbiome of cave systems has received less attention than macroscopic fauna and prokaryotic microbiomes. Surveys of these organisms are challenging as they are difficult to find due to their low population density and they are often much smaller than their surface‐dwelling counterparts, making them difficult to identify. When systematic microscopic screens of protists have been conducted, numerous uncharacterized species were found (Coppellotti and Guidolin, 2003; Bastian et al., 2009; Sigala‐Regalado et al., 2011; Baković et al., 2019). Metabarcoding has recently been used in some eukaryotic surveys, including cave walls in a freshwater algae biodiversity survey of the Hawaiian archipelago (Sherwood et al., 2014), surveys of green algae and diatoms in show cave Lampenflora (Pfendler et al., 2018; Burgoyne et al., 2021), soils from the Manao‐Pee Cave in Thailand (Wisechart et al., 2019), and a survey of microbial mats growing in low oxygen airbells in Movile Cave Romania (Reboul et al., 2019). Microeukaryotes are essential components of every food web (Corliss, 2002) as they link bacteria and small metazoans so describing them is important for a full understanding of the underground food web. The purpose of this study was to attempt a survey of the eukaryotic microbiome of two relatively shallow temperate caves in the Appalachian mountain chain using environmental DNA metabarcoding. Our goals were to refine sampling techniques to enhance the recovery of these organisms and attempt to identify as many barcodes as possible. We hope this study can serve as a productive step in the exploration of cave eukaryotic microbial ecology.

Results and discussion

Sampling strategy

Our first attempts to detect microeukaryotes using metabarcoding utilized DNA extracted from 49 L of filtered water collected from Panel and Bolling Caves from 2016 to 2019. These were screened for microeukaryotes using the ‘universal’ 18S rDNA primer set and collectively, 54 unique amplicon sequence variants (ASVs) were detected from all samples from both caves. Hoping to increase yields, we attempted a different year‐long assay using periphytometers and DNA extracted from glass slides submerged in cave water, which yielded 729 unique ASVs. For Bolling Cave 261 ASVs were detected using the periphyton sampler and 16 ASVs from scooping water. For Panel Cave 468 ASVs were detected using the periphyton sampler and 38 from scooping water. There was no overlap in ASVs detected from the two sampling methods so the two were combined for a total of 783 that were used for all subsequent analyses.

Our yields using periphytometers were comparable to another 18S rDNA metabarcoding survey from Movile cave, Romania (Reboul et al., 2019) which reported microeukaryotes from hypoxic air pockets and noted that they had higher yields from ‘fresh’ samples, meaning microbial mats whose DNA was extracted very soon after collection as opposed to mats that had been cultured for 2 months. In our study, the fact that the glass slides were immersed in DNA extraction buffer immediately after removal from the periphytometer in the cave may have also contributed to the higher ASV yield than those from filtered water.

ASV identification

GenBank records were used to create a reference tree from a combination of taxa predicted to be the most closely related to the cave ASVs as well as other taxa from microeukaryotic groups (Supplemental File 1). When a phylogenetic tree was produced using the reference sequences and the cave ASVs, 502 of the ASVs were broadly identifiable (summarized in Fig. 1A, all listed in Supplemental File 2). This approach provided taxonomic identities for 394 of the ASVs (Supplemental File 3). Only 1.5% of the ASVs matched a named GenBank record with greater than 99% identity, others matched records named ‘environmental sample’ or ‘uncultured’. We also used the PR2 reference database to make taxonomic assignments and 13.2% could be identified to the genus level with 99% certainty, 17.4% with 95% certainty, or 20.4% at 90% certainty. This is similar to Reboul et al. (2019) who found that few of their OTUs matched curated/cultured 18S microeukaryotic rDNA sequences.

Fig. 1

Distribution of 18S rDNA ASVs found in Panel and Bolling caves.

A. Proportion of ASVs from each taxonomic group found in each cave from filtered water and periphytometers combined. Am. = Amoebozoa, Api. = Apicomplexa, Archae. = Archaeplastida, Cil. = Ciliophora, Dino. = Dinoflagellata, M. = Metazoa, Rhi. = Rhizaria, S. = Stramenopiles, Ukn. = Unknowns.

B. Distribution of the number of periphytometer collected ASVs and the relative proportion of each group, based on read counts, collected at each time point from each cave.

This lab has been collecting plankton from surface water samples from various sites in and around Southwestern Virginia (the location of Panel and Bolling Caves) since 2015 and has an archive of freshwater 18S rDNA sequences which, to date, has 3283 unique ASVs (Cahoon et al., 2018, Fawley et al., 2021, and unpublished data). Only 2.5% of the cave samples matched one of these sequences with greater than 99% identity.

Temporal differences

The microscope slides were installed/removed at three time points over the course of a year so we were able to roughly estimate temporal changes in the most common groups in each cave (Fig. 1B). Notably, Rhizarians were the only group detected in both caves at all three time points. Dinoflagellates were detected in both caves but only from slides left submerged from May to September or September to January. A larger diversity of microeukaryotes for the three time points was found in Panel Cave. The temporal differences in microeukaryote variety and abundance are intriguing but it is not possible to associate them with abiotic factors since none were measured during his project.

Alveolates

The microeukaryotic supergroup with the greatest ASV representation were the alveolates. Dinoflagellates were well represented in the cave samples with 72 ASVs forming a clade with Dinoflagellate reference sequences (Fig. 2 and Supplemental File 4). When compared to the PR2 protist database, the same 72 ASVs were classified as Dinoflagellates with >70% certainty and 51 of the ASVs with >90% certainty (Supplemental File 3). Dinoflagellates have only been recorded in one other cave survey and it was a minor component (Reboul et al., 2019). About 88% of described freshwater dinoflagellates are photosynthetic with very few being benthic or heterotrophic (Tang, 1996; Stoecker, 1999; Gómez, 2012), attributes that would be useful for cave‐dwelling species. The majority of the ASVs, 68, formed a clade with the genus Peridinium, which was confirmed by the PR2 taxonomic predictions which identified some ASVs as Peridinium volzii and Peridinium willei. Peridinium is a common freshwater genus (Carty and Parrow, 2015) but it is unknown whether they can heterotrophically maintain homeostasis as would be required in a cave. These data do suggest that conditions either in the caves or in the surface hydrological environment have been favourable for the radiation of this genus. Eleven of the ASVs matched Asulcocephalium miricentonis which was first characterized in a temperate freshwater pond in Japan (Takahashi et al., 2015).

Fig. 2

Cave 18S rDNA ASVs classified as dinoflagellates. A maximum likelihood phylogeny was generated using IQTREE and the ModelFinder option. Reference sequences are black, cave ASVs purple, ASVs collected from scooped water are marked with asterisks and an outgroup sequence is green. ASVs found in surface samples are enclosed by rectangles. Bootstrap values are listed at each node. The scale represents number of substitutions.

Reference sequences from Ciliophora formed a clade with 65 cave ASVs (Supplemental Fig. 1 and Supplemental File 5) all of which were classified as Ciliophora with >70% certainty (Supplemental File 3). The majority of those were within the Orders Hymenostomatida and Peritrichia. Comparisons with the PR2 database suggested 50 of the ASVs matched a record with greater than 90% certainty. In both microscopic and metabarcoding cave studies, ciliates were the most common microeukaryote (Sigala‐Regalado et al., 2011; Reboul et al., 2019). The most common ciliate in Panel and Bolling caves appear to be within the genus Tetrahymena. Three ASVs from Bolling and Panel caves had greater than 99% identity to Tetrahymena (T. eeiotti, T. bergeri and T. farleyi) and as many as 15 others formed a clade that may represent unnamed members of this genus. Tetrahymena sp. are ubiquitously present in aquatic environments where they graze on bacterial and viral particles (Verni and Gualtieri, 1997; Pinheiro et al., 2007). They are also a valuable environmental bioindicator representative of healthy aquatic environments (Maurya and Pandey, 2020). Their presence in our survey could mean this species may be useful in evaluating the conservation status of a cave ecosystem. Four of these Tetrahymena ASVs were also found in surface samples. Two ASVs closely matched the genus Choreotrichia and four other ASVs formed a clade consistent with the class Spirotrichea which are found in freshwater environments but are best known in the marine benthosphere where they consume bacteria and microalgae and are prey for small metazoans (Pierce and Turner, 1992; Calbet and Saiz, 2005; Santoferrara et al., 2017). One other ASV closely matched a Peritrichia sp. record which along with 19 other ASVs grouped within the class Oligohymenophorea, a species‐rich but undersequenced group of ciliates (Sun et al., 2021). To the best of our knowledge, there are no records of cave endemic Tetrahymena, Choreotrichia, or Peritrichia but their lifestyles and physiology (Lynn, 2017) would be suitable to life in a lightless environment as long as there are sufficient numbers of prokaryotic and microeukaryotic prey available

Fifty nine ASVs formed a clade with Apicomplexan reference sequences (Fig. 3 and Supplemental File 6). Most were in the order Eugregarinorida. Two ASVs, Boll_ASV_155 and Boll_ASV_197 formed a monophyletic branch within the Apicomplexan phylogeny without a closely associated GenBank record but PR2 based taxonomy suggested they were in the Order Colpodellida, a free‐living non‐parasitic branch of Apicomplexans. No individual ASVs closely matched reference sequences and only 29 were classified within the PR2 database (Supplemental File 3). The taxonomic predictions based on PR2 did suggest that several of the Eugregarinorida ASVs were assigned to the genus level with greater than 90% certainty and included Monocystis sp., Syncystis mirabilis and Paraschneideria metamorphosa. The ASV predicted to be closely related to Syncystis mirabilis was also found in surface samples. Apicomplexans are a widespread (Cavalier‐Smith, 2014) but undercharacterized group of microeukaryotes best‐known as parasites (Simdyanov et al., 2017; Janouškovec et al., 2019). A recent metabarcoding screen of surface‐dwelling apicomplexans found very low diversity in freshwater samples (del Campo et al., 2019) but they were the most common microeukaryote found in the soils of the Manao‐Pee Cave in Thailand (Wisechart et al., 2019) and three were found in anoxic microbial mats (Reboul et al., 2019). Our data along with other cave studies suggest these environments could be a source of apicomplexan diversity. We attempted to produce apicomplexan specific 18S rDNA barcodes using published primer sequences (Huggins et al., 2019) but failed to produce amplicons. We believe this failure was due to these primers having been designed to recognize Apicomplexan parasites using available GenBank records where there is a paucity of sequences similar to those found in these cave environments.

Fig. 3

Cave 18S rDNA ASVs classified as apicomplexans. A maximum likelihood phylogeny was generated using IQTREE and the ModelFinder option. Reference sequences are black and cave ASVs red. ASVs found in surface samples are enclosed by rectangles. Dashed lines represent unknown taxonomy within Apicomplexa. Bootstrap values are listed at each node. The scale represents number of substitutions.

Rhizarians

Within Rhizaria, 78 ASVs were identified using the phylogenetic approach (Supplemental Fig. 2 and Supplemental File 7) and 71 using the PR2‐based taxonomy approach (Supplemental File 3). Two of the barcode sequences had greater than 99% similarity to GenBank archived sequences suggesting they may be the species Orciraptor agilis – KF207875 and Neocercomonas sp. – AY884313, which was confirmed by PR2‐based taxonomy (Supplemental File 3). PCR amplifications using primers targeting foraminifera 18S rDNA were attempted but no amplicons were produced. The majority of the Rhizaria ASVs were Cercozoans which are most likely bacterial grazers (Burki and Keeling, 2014) which could survive in an oligotrophic environment. There were also ASVs classified as Vampyrellida, also known as predatory amoebae that parasitize other microeukaryotes (Hess, 2017; More et al., 2019) and Plasmodiophorida best known as plant parasites (Hwang et al., 2012). Rhizaria were reported as a significant component of microbial mats collected in the Movile cave but they were not the predominant group, except in a cultured sample (Reboul et al., 2019). Freshwater foraminifera have been microscopically identified from European karst caves (Mazei et al., 2012; Baković et al., 2019) but were not detected in our survey using either universal 18S rDNA or foraminifera targeted primers.

Stramenopiles

Fifty two ASVs were predicted to be Stramenopiles according to PR2‐based taxonomy (Supplemental File 3) while 30 formed a clade with stramenopile sequences using the phylogenetic approach (Supplemental Fig. 3 and Supplemental File 8) with all grouping in either Bacillariophyta (diatoms) or Chrysophyceae. Six of the 16 Bacillariophyta ASVs were found in surface samples while six of the 14 Chrysophyceae were also found on the surface. Two of the diatom ASVs had greater than 99% identity with Gomphonema micropus – JN790282 and Achnanthidium pyrenaicum – KY863466, while one of the chrysophytes matched the record for Uroglenopsis americana – MK153242. PR2‐based taxonomic predictions agreed with these identities at the genus level. The additional Stramenopiles identified by PR2‐based taxonomy were predicted to be in the Classes MAST‐12, Oomycota and Labyrinthulomycetes (Supplemental File 3). Since these species have chloroplasts, a primer pair targeting plastid 23S rRNA was used to confirm their presence in the cave samples. Ninety eukaryotic 23S sequences were identified (Supplemental Fig. 4 and Supplemental File 9). The vast majority of these appeared to be embryophytes and green algae. Nine were grouped within the diatom clade but none matched Gomphonema or Achnanthidium similar to the 18S rDNA barcodes, but four did group with reference sequences from the genus Synedra. The remaining five were not identifiable with the 23S barcode. Fifteen were chrysophytes. Diatoms are well documented in caves as they are a component of show cave lampenflora and a cave conservation concern (reviewed in Falasco et al., 2014; Piano et al., 2015; Pfendler et al., 2018; Burgoyne et al., 2021). Chrysophytes are not a common component of lampenflora or of cave surveys but are mixotrophs (Nicholls and Wujek, 2015) so it is conceivable they could survive in a lightless environment. The predominating hypothesis regarding diatoms is that they are regularly introduced into caves as living cells or cysts which can opportunistically colonize under artificial light and that there are no endemic cave species. We think it likely that the diatom sequences detected in our screen were from spores or cysts that had been transported from the surface. We believe this since we found no evidence of non‐photosynthetic diatoms (Kamikawa et al., 2018) which would be the only ones capable of colonizing a cave. Also, a high proportion of them, 40%, were also found in planktonic surface samples.

Amoebozoans

Twenty five ASVs grouped with amoebozoan reference sequences (Supplemental Fig. 5 and Supplemental File 10). All appeared to be within the phylum Discosea but none had greater than 99% identity with a reference sequence and none were found in surface samples. We only found non‐testate Discosea amoebae in our screen and none matched a GenBank record with greater than 99% similarity and none were found among surface ASVs. Non‐testate amoebae have been reported from Lascaux Cave in France but the largest number of amoebae identified from caves is testate amoebae (Mazei et al., 2012; Baković et al., 2019).

Are there endemic cave microeukaryotes?

The majority of the ASVs detected in our cave survey did not occur in surface water samples, except for photosynthetic stramenopiles. Based on this alone it is tempting to conclude that most of these ASVs were cave endemic species but there are limitations in our comparisons. For example, the library of local surface planktonic microeukaryotes may not represent a true comparison to the benthic cave samples collected on the glass slides. Therefore we cannot estimate the proportion of the microeukaryotes found in our survey that are truly cave endemics, and it may not be possible until a surface periphyton survey is completed. Another limitation is the lack of archived microeukaryotic freshwater barcode sequences, which reflects the lack of sequence information from many branches of microeukaryotes (del Campo et al., 2014) and is even more problematic in under‐sampled environments. Metabarcoding technology allows us to detect the diversity of microeukaryotes from caves and make relative comparisons but a major challenge will be to isolate and/or culture these organisms to expand our sequence databases.

Baković et al. (2019) argue that the hydrological connections between the surface and karst cave systems prevent spatial isolation which would be essential for speciation and endemism of water‐dwelling microeukaryotes. The caves we sampled could align with this hypothesis since these are relatively shallow caves and it was previously established that 86.9% of the planktonic prokaryotic microbes overlap with those in surface water, which along with dye trace experiments of Panel Cave shows that surface water is continuously entering and cycling through what is presumed to be an epikarst layer (Morse et al., 2021). We also detected photosynthetic eukaryotes such as Stramenopiles, green algae and embryophytes, which would only be present if introduced from the surface by way of percolating water, air movement, or animal vectors. Although both caves are wild and have no permanent artificial lighting, they are in a state park and guided tours of the caves are offered to visitors during summer months who could transport spores, cysts and pollen into these caves.

Conclusions

We used two sampling methods to survey microeukaryotes in two karst caves and found periphytometers with glass slides provided a greater number and diversity of ASVs than scooping and filtering water when sampling this low‐density fraction of the microbiome. Using both methods we identified 784 unique microeukaryotic ASVs from two caves representing a wide array of bacterial grazers and parasites that could survive in the cave environment as well as other transient cells and/or spores and cysts transported from the surface. This study provides a baseline survey of microeukaryotes for the relatively common karst caves found in temperate regions.

Experimental procedures

Samples were collected from Panel and Bolling Caves in the Natural Tunnel State Park in Scott County, VA, USA using two different methodologies, scooping water and submerging periphytometers in cave streams and pools. The attributes of each cave, hydrogeology and sample sites were described in Thompson et al. (2019) and Morse et al. (2021). Briefly, they are gated wild caves located approximately 1.2 km apart on opposite ridges along Stock Creek. To our knowledge, the caves are hydrologically separated and developed in lower to middle Ordovician‐aged carbonate rock (Miller and Brosge, 1954; Brent, 1963). Permission to perform this research and collect samples was granted by the Commonwealth of Virginia Department of Conservation and Recreation (Research and Collecting Permit NT‐RCP‐121819).

The scooping method consisted of collecting three litres of water in three sterile 1‐L screw‐cap bottles at each collection site at each time. Water was collected from still pools approximately 25 cm deep and running streams 5–10 cm deep. Each litre of water was filtered through a disposable microfunnel with a 0.45 μm mesh filter (Daigger & Co., Vernon Hills, IL, USA) within 24 h of collection in a lab at UVA‐Wise. Filters were stored at −20°C until DNA extractions could be completed. Filters were thawed and suspended in the extraction buffer and bead‐beater tubes provided with Qiagen's DNeasy PowerWater Kit (Germantown, MD, USA). Using this method, Panel Cave was sampled once in April 2016, once in April 2017, and monthly for 14 months from April 2018 to May 2019 for a total of 16 sampling events. Bolling Cave was sampled once in April 2016 and once in April 2017. The 2017 samples from Panel and Bolling Caves and the monthly 2018–2019 Panel Cave samples were used in previous studies exploring the prokaryotic microbial ecology of these caves (Thompson et al., 2019; Morse et al., 2021).

The second method utilized periphytometers, devices that provide an artificial substrate for waterborne microbes to adhere to and colonize. These are commonly used to monitor microbial biomass and diversity for environmental bioassessment assays of surface water (Aloi, 1990; Barbour et al., 1999). We used plexiglass periphytometers designed to hold standard sterilized glass microscope slides. Beginning in January 2020 the samplers were submerged in water in the same two locations that had been sampled by scooping in Panel Cave and the two locations in Bolling Cave. Water depth ranged from 10 to 50 cm depending upon the collection site and water levels within the cave which varied throughout the year. The substrate supporting the periphytometers was solid rock, pebbles, or sand. The glass slides were collected and replaced approximately every 4 months for one calendar year, three collection events. At each retrieval time, the glass slides were removed from the samplers, placed directly into DNeasy PowerWater Kit bead‐beater tubes containing DNA extraction buffer at the sampling site, shaken vigorously by hand, and then transported out of the cave. Samples were stored in light tight bags for transport out of the caves. Tubes were taken to a laboratory at UVA‐Wise and DNA extractions were completed the day of each collection.

Samples of DNA from all the replicates collected from each cave at each time (e.g. six per cave for each periphyton collection event) were individually used as templates for three PCR reactions. These replicates were combined to minimize PCR bias. Phusion DNA polymerase (Thermo Fisher, Waltham, MA, USA) was used for all amplifications. The V4 region of the 18S rDNA gene was amplified using ‘universal’ microeukaryote primers, TAReuk454FWD1 (5′CCAGCASCYGCGGTAATTCC) and TAReukREV3 (5′ACTTTCGTTCTTGATYRA) (Stoeck et al., 2010), as well as the ‘universal’ 23S rDNA primers, p23SrV_f1 (5′ GGACAGAAAGACCCTATGAA) and p23SrV_r1 (5′ TCA GCCTGTTATCCCTAGAG) (Sherwood and Presting, 2007), apicomplexan‐targeted 18S rDNA primers, ApicomplexF: (5′‐CRAGGAAGTTTRAGGCAATAACAG) and ApicomplexR: (5′‐CTAGGCATTCCTCGTTHAHGATT) (Huggins et al., 2019), and foraminifera‐targeted 18S rDNA primers, S14F1 (5′‐CCATCTCATCCCTGCGTGTCTCCGAC) and S19F (5′‐GTACRAGGCATTCCTRGTT) (Morard et al., 2018). Primers were synthesized by Integrated DNA Technologies (Coralville, IA, USA). Amplicons were purified using a Select‐a‐Size DNA Clean and Concentrator kit (Zymo Research, Irvine, CA, USA). Amplicon mixtures were paired‐end sequenced using Genewiz's Amplicon EZ Illumina MiSeq service (South Plainfield, NJ, USA). Paired‐end reads were processed and ASVs were produced using DADA2 using default parameters except for truncLen = c(250,250) (Callahan et al., 2016). The ASV outputs were transferred to Geneious Prime (Biomatters, Auckland, NZ) where primer sequences and duplicate ASVs were identified and removed. Barcode data from this study have been archived in Mendeley (DOI: 10.17632/r5fv5txtfs.1).

For phylogenetic analyses a reference set of 18S rDNA sequences was constructed with GenBank records (https://www.ncbi.nlm.nih.gov/genbank/, listed in Supplemental File 1). First, accessions most similar to cave ASVs were identified using BLAST searches to roughly identify the taxonomic groups present in the samples. Second, once a taxonomic group was identified all GenBank records of 18S rDNA sequences from that taxonomic group (e.g. Rhizaria) with a defined species were downloaded and included in the reference set. The overall reference set and the reference sets used for each taxonomic group's phylogeny were refined over numerous cycles until a suitable reference set containing key named taxa, accessions with the most similarity to the cave ASVs, and only as much taxonomic redundancy as necessary was formed. ASV and the GenBank derived reference sequences were aligned using MUSCLE with default parameters (Edgar, 2004) and sequence ends trimmed within the Geneious environment. Maximum Likelihood analyses were completed using IQTREE 1.6.12 (Trifinopoulos et al., 2016) using either ModelFinder (Kalyaanamoorthy et al., 2017) or the HKY substitution model (Hasegawa et al., 1985). Ultrafast bootstrapping was used to support the ML topology (Hoang et al., 2018). These data are represented as cladograms in this report for the simplest expression of the hypothetical classifications of the unknowns. Identity and Similarity scores were calculated using LALIGN (Huang and Miller, 1991) using the web‐based version at (https://molbiol-tools.ca/Alignments.htm). ASVs were also classified based on the PR2 protist database, version 4.14.0 SSU (Guillou et al., 2013; del Campo et al., 2018) using QIIME2 (Bolyen et al., 2019) functions to trim the database to the V4 SSU barcode region [feature‐classifier extract‐reads], train the classifier [feature‐classifier fit‐classifier‐naïve‐bayes] and assign taxonomy to the ASVs [feature‐classifier classify‐sklearn] using default parameters.

Supplemental Fig. 1. Cave 18S rDNA ASVs Classified as Ciliophora. A Maximum Likelihood phylogeny was generated using IQTREE and the HKY base substitution model. All reference sequences were drawn from the phylum Ciliophora according to records archived in GenBank. Reference sequences are black, the cave ASVs red, ASVs collected from scooped water are marked with asterisks, and an outgroup is green. ASVs found in surface samples are enclosed by rectangles. Bootstrap values are listed at each node. The scale represents number of substitutions.

Click here for additional data file.

Supplemental Fig. 2. Cave 18S rDNA ASVs Classified as Rhizarians. Maximum Likelihood phylogeny generated using IQTREE and the HKY substitution model. Reference sequences are black, cave ASVs purple, and an outgroup sequence is green. Dashed lines represent unknown taxonomy within Rhizaria. Im. = phylum Imbricatae. Bootstrap values are listed at each node. The scale represents number of substitutions.

Click here for additional data file.

Supplemental Fig. 3. Stramenopiles identified by their similarities to archived 18S rDNA sequences. Maximum Likelihood phylogeny generated using IQTREE and HKY base substitution model. Bacillariophyta ASVs are tan and Chrysophytes red. ASVs collected from scooped water are marked with asterisks. ASVs found in surface samples are enclosed by rectangles. Bootstrap values are listed at each node. The scale represents number of substitutions.

Click here for additional data file.

Supplemental Fig. 4. Stramenopiles identified by their similarities to archived 23S rDNA sequences. This Maximum Likelihood phylogeny was generated using IQTREE and HKY base substitution model. ASVs found in surface samples are enclosed by rectangles. Bootstrap values are listed at each node. The scale represents number of substitutions.

Click here for additional data file.

Supplemental Fig. 5. Cave 18S rDNA ASVs Most Closely Related to Amoebozoan Reference Sequences. Maximum Likelihood phylogeny generated using IQTREE and the ModelFinder option. Reference sequences are black and cave ASVs orange. Bootstrap values are listed at each node. The scale represents number of substitutions.

Click here for additional data file.

Supplemental File 1. The list of GenBank files used to build the 18S rDNA reference trees.

Supplemental File 2. IQTREE output file for phylogenetic analysis of all nuclear 18S rDNA barcodes collected during this study.

Supplemental File 3. Eukaryotic taxonomic predictions for the cave nuclear 18S rDNA ASVs based on the PR2 database.

Supplemental File 4. IQTREE output file for Dinoflagellate nuclear 18S rDNA barcodes collected during this study.

Supplemental File 5. IQTREE output file for Ciliophora nuclear 18S rDNA barcodes collected during this study.

Supplemental File 6. IQTREE output file for Apicomplexan nuclear 18S rDNA barcodes collected during this study.

Supplemental File 7. IQTREE output file for Rhizaria nuclear 18S rDNA barcodes collected during this study.

Supplemental File 8. IQTREE output file for Stramenopile nuclear 18S rDNA barcodes collected during this study.

Supplemental File 9. IQTREE output file for plastid 23S rDNA barcodes collected during this study.

Supplemental File 10. IQTREE output file for Amoebozoan nuclear 18S rDNA barcodes collected during this study.

Click here for additional data file.