Stygobiont Diversity in the San Marcos Artesian Well and Edwards Aquifer Groundwater Ecosystem, Texas, USA

Abstract

The Edwards Aquifer and related Edwards-Trinity Aquifer of Central Texas, USA, is a global hotspot of stygobiont biodiversity. We summarize 125 years of biological investigation at the San Marcos Artesian Well (SMAW), the best studied and most biodiverse groundwater site (55 stygobiont taxa: 39 described and 16 undescribed) within the Edwards Aquifer Groundwater Ecosystem. Cluster analysis and redundancy analysis (RDA) incorporating temporally derived, distance-based Moran’s Eigenvector Mapping (dbMem) illustrate temporal dynamics in community composition in 85 high-frequency samples from the SMAW. Although hydraulic variability related to precipitation and discharge partially explained changes in community composition at the SMAW, a large amount of temporal autocorrelation between samples remains unexplained. We summarize potential mechanisms by which hydraulic changes can affect community structure in deep, phreatic karst aquifers. We also compile information on 12 other Edwards and Edwards-Trinity Aquifer sites with 10 or more documented stygobionts and used distance-based RDA to assess the relative influences of distance and site type on three measures of β-diversity. Distance between sites was the most important predictor of total dissimilarity and replacement, although site type was also important. Species richness difference was not predicted by either distance or site type.

1. Introduction

The karstic Edwards Aquifer of Central Texas (USA) supplies water for more than 2 million people [1] and is recognized for high stygobiont biodiversity [2]. The 10,500 km² Edwards Aquifer occurs in a broad arc of Cretaceous limestones that stretch approximately 400 km across Central Texas, USA, and is hydrologically connected to the 91,744 km² Edwards-Trinity Aquifer (Figure 1). Edwards-equivalent limestones also extend into northern Coahuila, Mexico [3]. Late Cretaceous through early Miocene uplift of the Edwards Plateau exposed Edwards limestones along the Balcones Fault Zone: a series of en-echelon, high-angle faults downthrown to the southeast [4,5]. Increased permeability along faults in exposed Edwards limestones allowed meteoric recharge and dissolution, forming complex west–east and southwest–northeast flowpaths within hydrologically connected segments of the aquifer [6]. Present-day flowpaths are overprinted on hypogenically-derived permeability [7]. To the south and east, freshwater in the aquifer is confined below non-karstic units and juxtaposed against a lower-permeability zone of sulfide-rich, saline water along a steep freshwater–saline water interface (FWSWI) [8]. Shallower flowpaths and less faulting dominate in the Edwards-Trinity system to the north and west of the Edwards Aquifer. This region also contains many active stream caves.

The San Marcos Artesian Well (SMAW, state well number 6701828), on the Texas State University, San Marcos campus, is a flowing freshwater artesian well completed in late 1895 or early 1896 in the confined portion of the Edwards aquifer. The well intersects a 1.5 m tall phreatic conduit at a depth of −59.5 m [9]. Nearby saline wells illustrate proximity (i.e., <100 m) to the FWSWI. Dye tracing showed hydraulic connectivity between nearby Ezell’s Cave (2.9 km to the southwest), the SMAW, and Deep Hole Spring (part of the San Marcos Springs Complex, 500 m northeast) [10]. From 2015 to 2020, discharge from the well averaged 16 L/s. During November 2013, 15-min continuous data documented average water properties of: temperature 22.3 °C (± 0.007), dissolved oxygen 5.3 mg/L (±0.01), and electrical conductivity 608 μS/cm (±0.5).

The SMAW and, to a lesser extent, springs, caves, and other wells in the Edwards and Edwards-Trinity Aquifers have been the focus of numerous studies ranging from taxonomic to ecological investigations. In the literature, the SMAW has been referred to as SWTSU Well, Texas State Artesian Well, and artesian well at San Marcos, and it is probably the source for most of the data for ‘San Marcos Springs’ in the first global list of subterranean biodiversity hotspots [11]. San Marcos Springs are hydrologically connected to the SMAW but inundated by a shallow reservoir, making them more difficult to sample. Consequently, fewer stygobionts are documented from San Marcos Springs relative to the SMAW.

The SMAW was the site of the first biospeleological investigations in Texas with the descriptions of the salamander Eurycea rathbuni (Stejneger, 1896), the shrimp Palaemon antrorum (Benedict, 1896), the isopod Cirolanides texensis (Benedict, 1896), and the amphipod Stygobromus flagellatus (Benedict, 1896). Soon after, Ulrich [12] described the isopod Lirceolus smithii (Ulrich, 1902) and two species of cyclopoid copepods. After those initial descriptions, taxonomic work in the Edwards Aquifer slowed until the late 1970s when Glenn Longley began a second phase of investigation with systematic sampling of the well via drift nets. Longley initiated important collaborations with taxonomists including John Holsinger (Amphipoda) and Robert Hershler (Gastropoda) and supported graduate students such as Henry Karnei that investigated Edwards Aquifer biodiversity. Through his collaboration and directorship of the Edwards Aquifer Research and Data Center, San Marcos, TX, USA (EARDC), created in 1979, the number of stygobionts recorded from the SMAW increased from eight to 25 between 1976 and 2000. Holsinger and Longley [9] and Longley [2] illustrated that the SMAW had a globally diverse stygobiont fauna, although previous studies had illustrated that faunal composition at the well was apparently distinct compared to other USA stygofauna sites [12,13,14,15,16,17]. Holsinger [15] and Holsinger and Longley [9] emphasized the presence of both marine and freshwater derived species at the site.

A third phase of biologic investigation began in the mid-2010s, when Benjamin Schwartz became director of the EARDC upon Longley’s retirement. Schwartz also facilitated systematic sampling and taxonomic collaborations, most importantly with Okan Külköylüoğlu (Ostracoda) [18,19,20,21,22,23]. He also initiated morphometric and molecular analyses [23,24] to identify new records of previously described species. Since 2015, published stygobiont richness at the SMAW has increased from 26 to 39 (

Table 1. Stygobionts from the San Marcos Artesian Well and basis of record. § = single site endemic; † = type locality.

Phylum	Class	Order	Family	Species	Reference
Nematoda	Enoplea	Trichocephalida	Capillariidae	Amphibiocapillaria texensis Moravec and Huffman, 2000†	[52]
Acanthocephala	Eoacanthocephala	Neoechinorhynchida	Dendronucleatidae	Dendronucleata americana Moravec and Huffman, 2000†§	[52]
Platyhelminthes	Rhabditophora	Tricladida	Kenkiidae	Sphalloplana mohri Hyman, 1938	[53]
Annelida	Clitellata	Arhynchobdellida	Erpobdellidae	Erpobdella sp.	[2]
		Lumbriculida	Lumbriculidae	gen. sp. undet.	Steve Fend, pers. comm.
Mollusca	Gastropoda	Littorinimorpha	Cochliopidae	Phreatodrobia micra (Pilsbry and Ferriss, 1906)	[54]
				Phreatodrobia nugax (Pilsbry and Ferriss, 1906)	[55]
				Phreatodrobia plana Hershler and Longley, 1986	[55]
				Phreatodrobia rotunda Hershler and Longley, 1986	[55]
Arthropoda	Ostracoda	Podocopida	Candonidae	Cabralcandona mixoni Külköylüoğlu, Yavuzatmaca, Akdemir, Schwartz and Hutchins, 2019†§	[18]
				Candonopsis sp.	Okan Külköylüoğlu pers. comm.
				Comalcandona sp.	Okan Külköylüoğlu pers. comm.
				Cypria lacrima Külköylüoğlu, Akdemir, Yavuzatmaca, Schwartz and Hutchins, 2017†§	[19]
				Lacrimacandona wisei Külköylüoğlu, Yavuzatmaca, Akdemir, Schwartz and Hutchins, 2017†§	[20]
				Namiotkocypria haysensis Külköylüoğlu, 2018†§	[56]
				Rugosuscandona scharfi Külköylüoğlu, Akdemir, Yavuzatmaca, Schwartz and Hutchins, 2017†	[21]
				Ufocandona hannaleeae Külköylüoğlu, Yavuzatmaca, Akdemir, Schwartz and Hutchins, 2017†§	[22]
				Schornikovcandona bellensis Külköylüoğlu, Yavuzatmaca, Akdemir, Diaz and Gibson, 2017	Okan Külköylüoğlu pers. comm.
			Darwinulidae	Darwinula sp.	Okan Külköylüoğlu pers. comm.
				Vestalenula sp.	Okan Külköylüoğlu pers. comm.
	Hexanauplia	Cyclopoida	Cyclopidae	Cyclops cavernarum Ulrich, 1902†§	[12]
				Cyclops learii Ulrich, 1902†§	[12]
		Harpacticoida		gen. sp. undet.	Diego Figueroa, pers. comm.
	Malacostraca	Bathynellacea	Bathynellidae	Hobbsinella edwardensis Camacho, Hutchins, Schwartz, Dorda, Casado and Rey, 2017†	[57]
		Thermosbaenacea	Monodellidae	Tethysbaena texana (Maguire, 1965)	[17]
		Isopoda	Microcerberidae	Microcerberidae sp. 1	[58]
				Microcerberidae sp. 2	pers. obs.
			Asellidae	Lirceolus hardeni Lewis and Bowman, 1996	[24]
				Lirceolus pilus (Steeves, 1968)	[24]
				Lirceolus smithii (Ulrich, 1902)†	[12]
				Lirceolus sp.	William Coleman pers. comm.
			Cirolanidae	Cirolanides texensis Benedict, 1896†	[13]
				Cirolanides wassenichae Schwartz, Hutchins, Schwartz, Hess and Bonett, 2019†	[59]
				Cirolanides sp.	pers. obs.
		Amphipoda	Artesiidae	Artesia subterranea Holsinger, 1980†	[9]
			Crangonyctidae	Stygobromus bifurcatus (Holsinger, 1967)	[23]
				Stygobromus flagellatus (Benedict, 1896)†	[13]
				Stygobromus longipes (Holsinger, 1966)	[23]
				Stygobromus russelli (Holsinger, 1967)	[9]
			Hadziidae	Allotexiweckelia hirsuta Holsinger, 1980†	[9]
				Holsingerius samacos (Holsinger, 1980)†	[9]
				Texiweckelia texensis (Holsinger, 1973)†	[15]
				Texiweckeliopsis insolita (Holsinger, 1980)†	[9]
			Parabogidiellidae	Parabogidiella americana Holsinger, 1980†	[9]
				Parabogidiella sp. Holsinger, 1980	[9]
			Seborgiidae	Seborgia relicta Holsinger, 1980†	[9]
		Decapoda	Palaemonidae	Calathaemon holthuisi (Strenth, 1976)	[60]
				Palaemon antrorum (Benedict, 1896)†	[13]
	Arachnida	Trombidiformes	Chappuisididae	Chappuisides sp.	Ian Smith, in prep.
				Uchidastygacarus sp.	Ian Smith, in prep.
			Halacaridae	gen. sp. undet.	Ian Smith, in prep.
			Mideopsidae	gen. sp. undet.	Ian Smith, in prep.
			Nudomideopsidae	Allomideopsis wichitaensis (Smith, 1990)	Ian Smith, in prep.
	Insecta	Coleoptera	Dytiscidae	Haideoporus texanus Young and Longley, 1976†	[16]
Chordata	Amphibia	Urodela	Plethodontidae	Eurycea rathbuni (Stejneger, 1896)†	[14]

).

Longley [2] hypothesized that the Edwards Aquifer foodweb is not supported by allochthonous organic matter from the surface, but rather by ‘fossil’ organic matter originating at depth in the saline portion of the aquifer. Longley’s hypothesis foreshadowed later studies that revealed the role of chemolithoautotrophy within the aquifer. Birdwell and Engel [25] characterized microbially derived dissolved organic matter along the FWSWI with a chromophoric signature distinct from terriginous surface and soil porewaters. Gray and Engel [26] identified microbial communities along the FWSWI with taxonomic composition similar to other chemolithoautotrophic systems. Finally, Hutchins et al. [27] reported isotopic signatures of carbon that indicated that chemolithoautotrophic production, in addition to photosynthetic organic matter, supports the metazoan community at the SMAW. They also suggested that chemolithoautotrophy might facilitate reduced extinction rates during climatically unfavorable periods. Hutchins et al. [28] suggested that, as a spatially and temporally stable food source, chemolithoautotrophy might support high biological diversity by increasing resource exploitation and reducing competition. By combining stable isotope and mouthpart morphologic data, the authors illustrated trophic niche partitioning among amphipod species, making the SMAW one of a small but growing number of sites with evidence of niche partitioning among stygobionts [29,30,31].

Given the hydraulic and geologic complexity of the Edwards Aquifer, the SMAW likely integrates water and stygobionts from multiple flowpaths, discrete locations, and microhabitats within the aquifer. Because hydraulic conditions along any discrete flowpath vary in response to precipitation and antecedent conditions, groundwater assemblage composition at the well may vary temporally as well. Temporal dynamics in groundwater community structure have been investigated in alluvial aquifers [32], epikarst [33], and karst aquifers [34,35,36], although to our knowledge, studies in the latter have been limited to vadose and shallow saturated systems rather than deep, phreatic sites. In high gradient, hydraulically ‘flashy’ karst aquifers, studies have emphasized the role of flood pulses in the ‘spatial redistribution’ of species [34,35] and how species-specific responses depend on hydraulic differences in microhabitats (e.g., transmissive conduits versus peripheral fracture networks). However, as Gibert et al. [34] noted, flow conditions may be more stable in the phreatic zone. Therefore, it is unclear whether stygobionts in more hydraulically stable, phreatic aquifers exhibit similar temporal variability. If so, then from a species-accumulation perspective, sample events might capture different components of a dynamic groundwater meta-community and species richness more likely describes the meta-community at a local, rather than site scale.

At a regional scale, the SMAW and surrounding groundwater environments are part of the larger Edwards Aquifer Groundwater Ecosystem. We use the term ‘Edwards Aquifer Groundwater Ecosystem’ to highlight the aquifer as a spatially extensive, discrete ecosystem with unique hydrogeologic, biological, and socio-economic elements. In contrast, conceptualizations of the aquifer have tended to emphasize isolated components of the system (i.e., groundwater as a resource, single sites as critical habitat for subsets of species). Although other sites in the aquifer have not been investigated as thoroughly as the SMAW, previous work has illustrated that diverse assemblages of stygobionts occur throughout the aquifer system [17,37]. Multiple diverse groundwater sites within a contiguous but heterogeneous aquifer system afford opportunity to investigate patterns of beta-diversity at the aquifer scale. Previous work has demonstrated that globally, most stygobionts are small-range endemics [38], and that overall diversity from local to continental scales is explained in part by the species-turnover component of beta (β)-diversity [39]. Although distance between sites, coupled with small ranges and limited dispersal potential [40], is a parsimonious explanation for species turnover, environmental differences among sites have also been proposed to explain differences in species richness and composition in epikarst copepod communities [41].

Here, we present a species list for the SMAW, using published and unpublished data. We also investigate temporal dynamics in community composition at the site via high-frequency sampling across 3 years and 85 samples. We hypothesized that samples from the well do not suggest a temporally stable community, but rather, show temporal variability, as in other groundwater systems. Specifically, we hypothesized that assemblage structure at the SMAW would vary seasonally and in response to precipitation-driven hydraulic changes. To our knowledge, this is the first assessment of temporal community dynamics within a deep phreatic karst system. We also investigate species richness at the SMAW within the context of the greater Edwards Aquifer Groundwater Ecosystem, highlighting other diverse sites. We assess whether (1) regional-scale patterns of β-diversity are explained by species turnover (replacement) or regional/ site-based differences in species richness and (2) whether those dissimilarities are explained by distance or site-type (i.e., springs versus wells and caves). We predicted that species turnover rather than differences in species richness drive β-diversity patterns, and that species turnover would be affected by distance rather than site type.

2. Materials and Methods

2.1. SMAW Diversity and Temporal Dynamics

We conducted a literature review to compile a stygobiont species list for the SMAW. In a few instances, taxonomic experts were consulted to determine whether species should be considered stygobionts (as defined by Trajano and de Carvalho [42]). Unpublished, undescribed taxa were also included, based on communication with collaborators and personal observations in the authors’ taxonomic area of expertise. For temporal analysis, 85 samples were collected via 60 µm drift nets attached to the outflow of the well for between 24 and 72 h between 13 February 2013 and 20 November 2015 (Table S1). Samples were preserved in 95% EtOH and sorted at 10× magnification. Species were identified to the lowest taxonomic level by the authors or taxonomic experts (see acknowledgements) although some undescribed taxa were lumped as a single taxon (e.g., Microcerberidae, Trombidiformes). Voucher specimens for most taxa are retained in the Aquifer Biodiversity Collection of the EARDC, Texas State University, San Marcos, TX, USA.

Statistical tools in R v4.0.3 were applied to explore stygobiont time-series data from the SMAW. A Ward’s minimum variance dendrogram based on species abundances was created using the vegan package following the method of Borcard et al. [43]. A graph of silhouette widths was visually examined to estimate the optimal number of clusters, which was statistically assessed via analysis of similarity using a Bray–Curtis similarity matrix and 9999 permutations. Species associated with each group at p < 0.05 were identified using Indicator Species Analysis run with 999 permutations.

To disentangle the potential effects of discharge, season, and unexplained temporal influences on community composition, redundancy analysis was performed using the vegan package. Independent variables included season, discharge, and distance-based Moran’s eigenvector maps (dbMEM) derived from time (day) since the first sampling event (T₀) (Table S1). Season was coded as a categorical variable according to Kollaus and Bonner [44]: winter = December–February, spring = March–May, summer = June–August, fall = September–November. Discharge data were derived from mean daily discharge recorded by the U.S. Geological Survey at the San Marcos River gaging station USGS 08170500 [45], converted to liters per second. The gaging station is approximately 0.5 km downstream from the San Marcos Springs, and discharge at the station is a surrogate for local aquifer levels, which are correlated with average flow velocities in the Edwards Aquifer [46]. Distance-based Moran’s eigenvector maps (dbMEMs) were created using the method of Legendre and Gauthier [47] as implemented using the dbmem function in the adespatial package. dbMEMs were derived from a distance matrix of the number of days between sampling events and represent a spectral decomposition of the temporal relationships among samples [47]. Significance of dbMEMs was assessed using the moran.randtest function in the adespatial package, and only dbMEMs significant at p < 0.05 were used in the RDA. Periodicity of significant dbMEMs was not calculated because of many missing values over the sampling period. RDAs on discharge and season only, and on the dbMEMs only, were performed prior to a global RDA with both sets of independent variables and variance partitioning. The species data matrix was Hellinger transformed and singletons and doubletons were removed, but data were not detrended prior to analysis. Significance of RDAs, independent variables, and canonical axes was assessed using the anova.cca function in vegan, with 1000 permutations. Variance explained by discharge and season versus dbMEMs was assessed using the varpart function in vegan.

2.2. Edwards Aquifer Groundwater Ecosystem β-Diversity

To place the SMAW community within the broader context of the Edwards Aquifer Groundwater Ecosystem, distance-based RDA (dbRDA) was performed to assess the influence of distance, aquifer pool, and site type on β-diversity in R v.4.0.3. A literature review identified additional Edwards Aquifer and Edwards-Trinity Aquifer sites with high diversity (with an arbitrary cut-off of 10 or more stygobiont species). Primary and grey literature resources and previously unreported records of described species represented by specimens in the collections of the author and colleagues were included (Table S2). Unpublished, undetermined taxa were not included. β-diversity (Jaccard index total dissimilarity) was estimated and partitioned into two components (replacement, richness) using the beta function in the package BAT [48]. Linear trends in β-diversity measures were assessed using linear regression against log-transformed Euclidean geographic distance between sites, with Bonferroni correction applied for multiple comparisons. Dissimilarity matrices were used as the response variable in dbRDA. Predictor variables included site type (spring or well/cave), aquifer pool, and a distance-based Moran’s eigenvector map (dbMem) derived from site coordinates (UTM). Because of high multicolinearity (VIF > 10) between aquifer pool and dbMem variables, aquifer pool was removed prior to analysis. The single dbMem derived from the coordinates of the 13 assessed sites was calculated using the dbmem function in the adespatial package. Because sites are unevenly spaced and often widely distributed, the truncation threshold was set at 217039.8 m, limiting assessment of spatial structure to broad spatial scales. dbRDAs were conducted using the capscale function in the vegan package. Significance of the dbRDA and predictor terms was assessed via the anova.cca function in vegan, with 1000 replications.

3. Results

3.1. SMAW Diversity and Temporal Dynamics

Currently, 55 taxa have been documented from the SMAW, including 39 species recorded in literature (Table 1), two of which remain undescribed (Parabogidiella sp. Holsinger, 1980 and Erpobdella sp.). An additional 16 stygobiont taxa (2 described species and 14 undetermined taxa) are reported here for the first time (Table 1). Crustaceans dominate the SMAW fauna, comprising 75% of documented species. These include the only thermosbaenacean in the United States [49] and a globally significant amphipod fauna [9] of 12 species in five families. Additionally, 11 ostracod and nine isopod species occur at the site. Unique soft-bodied taxa include the only North American stygobitic leech [50], two vertebrate parasites (see below), and five species in the gastropod genus Phreatodrobia. The beetle, Haideoporus texanus (Young and Longley, 1976) is one of five species of stygobitic dytiscid beetles in the United States (four of which are associated with the Edward Aquifer [51]). A single vertebrate, the Texas blind salamander, E. rathbuni, occurs at the well. The SMAW is the type locality for 25 of the published taxa (64%), and eight of these (21%) are single-site endemics (Table 1). New SMAW site records, (including species descriptions and new records for described species) have accumulated unevenly over time (Figure 2).

The stygobiont community at the SMAW is characterized by some taxonomic and ecological uncertainty. In addition to the two undescribed species already mentioned, the two copepod species reported from the well are poorly described and have been designated as nomen inquirendum or nomen dubium by some [61]. Two species, the nematode Amphibiocapillaria texensis Moravec and Huffman, 2000 and the acanthocephalan Dendronucleata americana Moravec and Huffman, 2000 are parasites of the salamander E. rathbuni. Both parasites complete their life cycles using stygobiont invertebrates and should be considered stygobionts (David Huffman, pers. comm.). Neither parasite species have been conclusively collected from the epigean salamander Eurycea nana Bishop, 1941, which occurs in the San Marcos Springs and San Marcos River.

Eighty-five samples collected between 2013 and 2015 contained 42,814 individuals with an average of 275 individuals and 15 taxa per 24 h period (Table S1). Incidental epigean taxa were removed prior to temporal analysis (n = 61). Undetermined individuals (n = 87), mostly ostracods, were also removed because they could not be confidently assigned a taxonomic identity. An additional 415 juvenile and damaged individuals were assigned to taxa present in the sample.

Log-transformed total abundances show a nearly normal distribution (Figure 3), and only two species in the 85 samples (Lirceolus pilus (Steeves, 1968) and an unidentified nematode) were found in only one or two samples, suggesting that rare species are mostly accounted for by the sampling effort. Three species (Phreatodrobia micra (Pilsbry and Ferriss, 1906), Lirceolus hardeni Lewis and Bowman, 1996, and Stygobromus bifurcatus (Holsinger, 1967)) are known from the site by single specimens but were not present in any of the 85 analyzed samples. The two parasitic taxa present in the aquifer were not detected in our sampling strategy, which did not involve dissection. Additionally, copepods and mites were not identified to species, so all copepods and mites were each lumped into a single category. A histogram of the number of samples in which taxa occur shows a bimodal distribution, suggesting that most species are either common or rare (Figure 3). The shrimp, P. antrorum, makes up 44% of individuals, and just four taxa: P. antrorum, Copepoda, Cypria lacrima Külköylüoğlu, Akdemir, Yavuzatmaca, Schwartz and Hutchins, 2017, and Texiweckeliopsis insolita (Holsinger, 1980) make up over 90% of individuals. The 23 most infrequent taxa (61% of all taxa) collectively make up less than 1% of the total number of individuals.

Samples form significant clusters based on species abundances (Figure 4). Two clusters produced an optimal silhouette width, but between two and five clusters had similarly high silhouette widths. ANOSIM confirmed that groupings of samples into five and two clusters were both significant (R: 0.45 and p < 0.001; R = 0.26 and p < 0.001, respectively). Several taxa were significantly associated (p < 0.05) with one or two clusters although some species were associated with two clusters that were not nearest neighbors (i.e., species could be associated with two clusters that were otherwise compositionally dissimilar to one another).

Although there was a significant temporal trend in the community data (p < 0.05), only a small proportion of variance was constrained by sampling date alone. Separate redundancy analyses (RDA) using only discharge plus season and only distance-based Moran’s eigenvector maps (dbMEMs) were both significant (F = 8.589, p = 0.001; F = 3.930, p = 0.001, respectively). A global RDA incorporating discharge, season, and dbMEMs was also significant (F = 9.274, p = 0.001). Variance partitioning between discharge plus season versus dbMEMs showed that shared explained variance was not significant. Discharge and season explained 32% of variance in community structure, and both variables were significant at p = 0.001 (F = 21.084 and 8.013, respectively). dbMEMs explained 18% of variance in community structure, and three of four dbMEMs were significant at p = 0.001 (F = 6.635–10.878). The first three axes of the global RDA were significant at p = 0.001 and cumulatively accounted for 45% of the total explained variance. The first axis illustrated a gradient between low-flow samples (primarily during spring and summer) with larger numbers of the shrimp P. antrorum and the amphipod T. insolita and high-flow samples (primarily during fall and winter) with higher numbers of copepods, the snail Phreatodrobia plana Hershler and Longley, 1986, and the amphipod Texiweckelia texensis (Holsinger, 1980) (Figure 5). The second axis illustrated unexplained temporal gradients (dbMEM1 and dbMEM3) between samples (primarily in the summer) with higher numbers of the ostracod C. lacrima and samples with higher numbers of the amphipods S. flagellatus and Seborgia relicta Holsinger, 1980 (Figure 5). Clustering of sites in RDA space (influenced by community composition and environmental variables) reflects clustering on the Ward’s minimum variance dendrogram (based only on community composition, Figure 4).

3.2. Edwards Aquifer Groundwater Ecosystem β-Diversity

Thirteen sites with 10 or more stygobionts were identified across the Edwards and Edwards-Trinity Aquifers (Figure 1, Table S2). Sites included three flowing artesian wells, one cave, and nine springs varying from 1st magnitude springs (e.g., Comal Springs) to a small, intermittent spring (Sessom Creek Spring). All sites are hydraulically connected to the phreatic zone of the Edwards or Edwards-Trinity Aquifers. Total dissimilarity and replacement increased with increasing distance (R² = 0.57 and 0.32, respectively) at p < 0.001. Differences in species richness did not exhibit a spatial trend. dbRDA revealed significant site-type and distance-based effects on total dissimilarity (F = 2.49, p = 0.001, R² = 0.20) and replacement (F = 2.95, p = 0.001, R² = 0.41), but not on differences in species richness. Thirty-three percent of variance in total dissimilarity was constrained by the two canonical axes (axis 1: 19%, axis 2: 14%), and both site type (F = 2.14, p = 0.008, loadings: axis 1 = 0.15, axis 2 = −0.99) and dbMem (F = 2.84, p = 0.001, loadings: axis 1 = 0.89, axis 2 = 0.46) were significant terms. Thirty-seven percent of variance in replacement was constrained by the two canonical axes (axis 1: 22%, axis 2: 15%), and both site type (F = 2.47, p = 0.003, loadings: axis 1 = 0.09, axis 2 = −1.00) and dbMem (F = 3.43, p = 0.001, loadings: axis 1 = 0.91, axis 2 = 0.41) were significant terms. For each RDA, the first axis describes differences in dissimilarity explained by distance between sites and the second axis describes differences in dissimilarity explained by site type (Figure 6).

4. Discussion

Despite over one century of intensive sampling and study, knowledge of diversity at the SMAW remains incomplete. ‘Orphan’ taxa (e.g., Cyclopoida and Harpacticoida) are completely or largely unassessed. Even in better-studied groups (e.g., Isopoda), undescribed taxa have been identified. We conclude that reported species richness is underestimated. The continuing increase in species recorded at the site over time (Figure 2) and the near absence of taxa that are important elements of many groundwater communities (e.g., Copepoda) supports this assertion. The semi-normal distribution of species abundances in high-frequency samples (Figure 3A) and species accumulation curves (data not shown) suggest that sampling has been adequate for known species. Additional species discovery at the well will mostly likely result from additional taxonomic assessment of the orphan taxa discussed above, and cryptic species [23,24]. Nevertheless, with 55 groundwater-obligate taxa (Table 1), the SMAW is the most diverse groundwater site in North America and among the most diverse sites globally [63]. Proposed explanations for high-biodiversity within the Edwards Aquifer Groundwater Ecosystem include the role of marine embayments producing relic taxa (i.e., a richer colonist pool sensu Cardinale et al. [64]) [9] and high rates of primary productivity [27] supported, in part, by chemolithoautotrophy. Long-term productive energy, which is linked to climate and climatic-variability, has emerged as an important driver of groundwater diversity patterns in Europe and North America [65,66]. South of Pleistocene ice sheets and permafrost, Central Texas climate was cooler and wetter during glacial periods [67], potentially having a positive effect on productivity in the region. The region has become warmer and dryer since the last glacial maximum and although pronounced aridity during the mid-Holocene Altithermal Period would have reduced surface productivity, Hutchins et al. [27] hypothesized that chemolithoautotrophic production may have mitigated effects of aridity-related changes in surface productivity on stygobionts. In Europe, habitat heterogeneity is also an important predictor of stygobiont richness at mid and southern latitudes [65] and may be important in the hydrogeologically complex Edwards Aquifer.

To our knowledge, our analysis of a high-frequency dataset spanning multiple seasons and flow regimes provides the first illustration of temporal variability in a deep phreatic aquifer community. Community composition-based clustering of samples suggests that some species abundances vary in synchrony over time, which would not be expected in random samples from a temporally stable community. The bimodal distribution of taxa frequencies in repeated samples (Figure 3B) is interesting and may reflect the presence of two or more metacommunities: one comprised of species that are ubiquitous in flowpaths intersected by the well (e.g., Lacrimacandona wisei Külköylüoğlu, Yavuzatmaca, Akdemir, Schwartz and Hutchins, 2017, P. antrorum, S. flagellatus, T. insolita), and one or more communities comprised of species in more remote habitats that are only infrequently washed out, typically during high flows. Alternatively, the distribution may not reflect spatial or temporal heterogeneity in community structure, but rather, biological differences (e.g., benthic versus pelagic habitat, or swimming ability) in species’ propensity to be expelled from the well, which would also produce a relationship between flow and sample composition.

Redundancy analysis illustrates that, like other groundwater systems, hydrologic regime plays an important role in structuring the SMAW stygobiont assemblage over time. Discharge was the most important predictor in the global RDA. The influence of discharge is also apparent when viewed alongside the community-based dendrogram (Figure 4), although several low-flow samples fall within otherwise high-flow clusters and vice versa. Samples 71 and 72 were high-flow samples collected immediately after an extended period of low-flow, and so may represent a ‘piston-effect’ in which rapid recharge by meteoric water pushes resident groundwater (and associated fauna) through the system [38]. Conversely, the low-flow samples 38, 40, and 83 all contained uncommon species. Since species richness at the SMAW is positively correlated with discharge (data not shown), uncommon species more strongly effect overall community composition in otherwise less-diverse, low-flow samples, and probably drive the clustering of these samples with more-diverse, high-flow samples. Sample 84, a high-flow sample clustering with low-flow samples, had a near absence of copepods, potentially reflecting an undetected sample processing error.

Why certain species appear to be associated with low or high flows is unclear. Gibert et al. [34] suggested that stygobionts are heterogeneously distributed across distinct microhabitats in aquifers, and flood pulses initiate transport of organisms in transmissive zones through the karst system. However, other direct and indirect mechanisms may also facilitate assemblage structure changes through demographic shifts or passive or active movement of species (

Table 2. Potential mechanisms by which hydraulic changes affect stygobiont assemblage composition in deep phreatic aquifers.

	Positive Association			Negative Association
	Hydraulic/ environmental effect	→	Community response	Hydraulic/ environmental effect	→	Community response
Increased flow	Activates intermittent/ alternative flowpaths and hydraulic connections	→	Animals transported from typically less-connected aquifer areas to sampling location	Increased flow within conduit	→	Animals within conduit seek refuge in adjacent fracture network
	Increased hydraulic gradient from fracture network to conduit network	→	Animals transported to sampling location from adjacent fracture network
	Increased sediment transport and conduit boundary velocities	→	Dislodges benthic/ interstitial organisms already within sampling area	Changes in nutrient availability, community composition, and geochemistry	→	Direct and indirect demographic (e.g., reproduction) or species interaction (e.g., ecologic release) effects
	Changes in nutrient availability, community composition, and geochemistry	→	Direct and indirect demographic (e.g., reproduction) or species interaction (e.g., ecologic release) effects
Decreased flow	Decreased flow within conduit	→	Animals within fracture network move into conduit	Changes in nutrient availability, community composition, and geochemistry	→	Direct and indirect demographic (e.g., reproduction) or species interaction (e.g., ecologic release) effects
	Changes in nutrient availability, community composition, and geochemistry	→	Direct and indirect demographic (e.g., reproduction) or species interaction (e.g., ecologic release) effects

). Analysis of hydrographs, geochemical dynamics, morphologic/ trophic patterns, and population dynamics within groups of concordant taxa may provide additional insight into causal relationship between hydrologic variability and stygobiont community dynamics.

As Gibert et al. [34] also acknowledged, groundwater flow alone does not explain temporal dynamics of groundwater communities. The significant contribution of temporal dbMEMs in our RDA demonstrates unexplained temporal dynamics in community composition at multiple scales. We did not attempt to correlate dbMEMs with potential explanatory phenomenon because we did not have a priori predictions about potential mechanisms, and because our dataset did not span multi-year cyclical weather oscillations like the Atlantic Multi-decadal Oscillation or El Niño–Southern Oscillation.

The 55 taxa from the SMAW represent about half of the approximately 102 stygobionts recorded from the Edwards Aquifer, and most occur at other sites. Within an extensive aquifer, assessment of biodiversity at local to regional scales, rather than at single sites, makes more sense ecologically and for aquifer management. Cardoso et al. [48] discuss issues with analyses of β-diversity based on incomplete/uneven sampling, including underestimation of similarity. However, assessment of randomized accumulation curves to control for sample effort (sensu Cardoso et al. [68]) was not possible in this study because species lists for sites other than the SMAW were based on literature and not samples. Consequently, β-diversity data are interpreted with the acknowledgement that uneven sampling effort across sites (and across taxa) obscures patterns. Given that caveat, we did observe increasing dissimilarity and replacement with increasing distance between sites, as predicted, and total dissimilarity and replacement was greater between springs and wells/caves than within site types (although the number of diverse wells and caves was limited in number and spatial extent compared to springs). However, biplot scores showed that for both total dissimilarity and replacement, distance between sites was more important than site type. Dissimilarity among habitat types was also a relatively unimportant component of β-diversity in European stygobionts [69]. Importantly, the species richness difference component of dissimilarity did not vary by distance or site type, suggesting the lack of ‘hotspot’ regions within the Edwards Aquifer with respect to location or site type, despite uneven sampling effort across regions. Because of uneven site distribution across the aquifer, however, fine-scale patterns may be obscured. For example, Hutchins et al. [27] detected a positive relationship between species richness and proximity to the FWSWI, although biodiverse sites far from the FWSWI (e.g., Caroline and Finegan-Blue Springs) raise questions about the importance of that relationship. Certainly, species richness at fine scales varies in response to hydraulic and geochemical properties, as evidenced by low- and high-diversity sites in close proximity to one another (pers. obs.).

Relative to other groundwater habitats, knowledge of stygobiont diversity in deep phreatic karst aquifers is lacking. The SMAW and the Edwards Aquifer Groundwater Ecosystem illustrates that springs and wells can be particularly productive sites for sampling these habitats, and we suspect that biodiversity within the aquifer is not anomalous relative to other deep phreatic karst aquifers on a global scale. Although spatial biodiversity patterns have received a good deal of attention from groundwater ecologists, increasingly sophisticated analytical methods [70] afford more opportunity to assess spatial and temporal patterns in community structure. In the face of global climate change and increasing anthropogenic pressures on groundwater ecosystems [71], analysis of spatial and temporal trends in groundwater communities will be increasingly important.