The Barton Springs segment provides water for ~60,000 people in the ~400 km² region that extends from the city of Kyle, Hays County, Texas, northward to the Colorado River in the city of Austin, Travis County, Texas (Krejca and Reddell 2019) (Fig. 1). The FWSWI defines the eastern boundary of the segment. To date, the phreatic and vadose zones of the Barton Springs segment contain at least 16 described species of stygobiontic invertebrates, as well as two endemic salamanders, the Barton Springs salamander, Eurycea sosorum Chippindale, Price & Hillis, 1993, and the Austin blind salamander, Eurycea waterlooensis Hillis, Chamberlain, Wilcox & Chippindale, 2001 (Krejca and Reddell 2019), both listed as endangered species under the U.S. Endangered Species Act.

Groundwater in the Barton Springs segment discharges at several discrete spring openings that form the Barton Springs Complex at the northern boundary of the aquifer segment. These springs are the terminus for regional flow paths, and each spring has different saline water contributions based on the flow path. Old Mill Spring has the highest saline contribution, estimated to be from 3 to 9% (Garner 2005; Mahler et al. 2006; Hunt et al. 2019). Eliza and Main Barton springs receive a mixture of three flow paths (i.e., Manchaca, Sunset Valley, and Cold Springs) depending on hydrologic conditions, with a ~0.5% saline contribution for Eliza Spring and from 1 to 6% for Main Barton Spring. There is no known saline contribution for Upper Barton Spring. There are five surface drainage sub-basins in the Barton Springs Watershed: Barton, Williamson, Slaughter, Bear, and Onion Creeks. Surface flow in these creeks each contribute up to 36% of the total aquifer recharge, as well as surficial OM input to the aquifer (Slade Jr 2014; Hauwert 2016), with the remaining recharge occurring diffusely across upland regions between the creeks. A series of wells in the recharge, confined, and saline zones, including several wells along the FWSWI, allow direct access for sampling the deep aquifer along and across the regional flow paths.

Water and stygofauna were collected throughout the Barton Springs segment from October 2019 to May 2021 (Table 2). Surface waters in Barton, Bear, and Onion Creeks were sampled from along the western Barton Springs segment boundary near the upstream margin of the recharge zone. Slaughter and Williamson Creeks were dry during the study period and were not sampled. Four springs in the Barton Springs Complex were sampled, including Eliza, Old Mill, Main Barton, and Upper Barton. A fifth spring in the Complex, Cold Spring, was not sampled because it is in a different sub-basin (Hauwert et al. 2004). Twenty-two wells spanning a range of zones and hydrologic conditions (i.e., recharge, unconfined, shallow-confined, deep-confined, and saline), and at varying distances from the FWSWI were sampled at least once over the study period. At least two borehole volumes were purged and field physicochemistry (e.g., dissolved oxygen, electrical conductance, temperature, and pH) was monitored with a Eureka Manta 30+ water quality sonde or YSI ProDSS sonde for stability in a flow-through cell prior to well sampling (USGS 2006). Depending on the well, stabilization took from 0.5 to 5 hours.

At each sampling location, field parameters were collected with a Eureka Manta 30+ water quality multiprobe sonde placed in flowing stream water or in the spring orifice. Samples for fine particulate OM (FPOM) were collected at each of the recharge stream sites and at spring orifices by filling a pre-rinsed (with deionized water, followed by site water) 20 L carboy. At each well, bulk water samples were collected in 1 L bottles that were acid-washed, triple-rinsed with deionized (DI) water, and then rinsed twice with sample water. At all sites, dissolved sulfide, sulfate, ammonium, and dissolved oxygen (DO) concentrations were measured directly in the field using CHEMetrics® kits and a V-2000 multi-analyte photometer. Water samples were transported on ice and stored at < 4 °C until analysis within 48 hours of collection.

To sample the microbial community at each stream and spring location, up to 2 L of water was manually filtered through at least two sterile 0.22 μm Sterivex Millipore polyethersulfone (PES) filter cartridges using a syringe. Filters were preserved with 2 mL of 2% cetyltrimethyl ammonium bromide (CTAB) (Engel et al. 2020), then kept on ice during transport and stored at -20 °C until DNA extractions.

Metazoans were collected from the springs using 100 μm mesh drift nets set at spring openings. Some nets had sample cups attached to the ends to reduce turbulence and provide refuge for salamanders and fragile metazoans. Drift nets were deployed and typically checked every 12 hours, when samples were examined to remove and release live salamanders but retain and preserve invertebrates. Before release, salamanders were identified to species level and tail clips were taken from each individual using a clean razor blade. Additional salamander tail clips were collected by City of Austin field crews. Tail clips were preserved in 95% ethanol (EtOH).

Except for two wells with anoxic conditions (Kai Frech and United Gas; Table 2), invertebrate collection was attempted at all wells using plastic bottle traps modified from Fenolio et al. (2017) and phreatobiologic nets modified from Cvetkov (1968). Traps contained poly-cotton mophead material and bait (e.g., frankfurters, peanuts, raw potatoes, and/or cheese), but invertebrates could not directly access nor consume the bait. Traps were recovered after 48 to 72 hours and contents were preserved on-site in 95% EtOH. Phreatobiologic nets (100 μm mesh-size) were lowered and raised five times through the borehole water column after trap removal. Contents were preserved in 95% EtOH after each tow.

Metazonas were sorted and identified to the lowest practical taxonomic level using a Nikon dissection microscope, relevant keys (e.g., Holsinger and Longley 1980; Merritt et al. 2019), and various reference images of stygobionts (courtesy of Pete Diaz, United States Fish and Wildlife Service, San Marcos Aquatic Resources Center, San Marcos, Texas). Chironomidae were identified to subfamily (i.e., Chironominae, Orthocladiinae, and Tanypodinae).

To determine if cyclopoid copepods were of epigean or hypogean origin, we examined 100 individuals from Eliza Spring to determine if they possessed an eye spot, a lack of which suggests hypogean origin (Hose et al. 2014).

In the laboratory, dissolved ions and nutrients were analyzed from samples filtered through pre-ashed Pall A/E filters into clean HDPE bottles. Alkalinity, as total titratable bases dominated by bicarbonate at the range of measured pH values, was measured via end-point titration with 1.6 N sulfuric acid (Rounds and Wilde 2012). Dissolved major cations and anions were measured from filtered and acidified water samples using Dionex ICS-1600 ion chromatographs (Bannockburn, IL). Nitrate was also quantified on a Varian Cary 50 UV-Vis spectrophotometer by using second-derivative spectroscopy (Crumpton et al. 1992). Charge balance was calculated for each analysis from the alkalinity, cation, and anion concentrations, with +/- 10% considered a complete analysis. Analyses exceeding this range are reported for comparison. Unfiltered water was placed into clean, acid-washed HDPE bottles to determine total nutrient concentrations following Wetzel and Likens (2000) to measure ammonium by using the sodium hypochlorite method, total nitrogen (TN) as nitrate by using the Crumpton et al. (1992) method after digestion with alkaline potassium persulfate, and soluble reactive phosphorus (SRP) as phosphate and total phosphorus (TP) by using the ascorbic acid method, with TP being measured after digestion with potassium persulfate and the molybdenum blue method. Non-Purgeable Organic Carbon (NPOC) was quantified by using a Shimadzu TOC-V CSH and Total Organic Carbon Analyzer (Columbia, Maryland), following a modified NPOC method to yield the dissolved organic carbon fraction (USGS 2015).

Fluorescence spectroscopy of chromophoric dissolved OM (CDOM) was analyzed from water filtered using 0.45-μm-micron Fisher brand nylon syringe filters with both excitation and emission slits set to 5 nm and 0.1 sec integration time in a quartz 1 cm cuvette on a Horiba Jobin Fluoromax-4 (Birdwell and Engel 2009, 2010). Excitation-emission matrix spectra were collected as the signal of the sample normalized to the reference detector ratio across a series of 43 emission scans from 250 nm to 550 nm in 2.5-nm increments at 5 nm excitation wavelength intervals from 240 nm to 450 nm in 5 nm increments (Birdwell and Engel 2009, 2010; Cory et al. 2010). Absorbance was measured on a Thermo Scientific Evolution 201 UV-Visible spectrophotometer from 240 nm to 820 nm in 0.5-nm increments with ultrapure DI water as a blank (>18.1 mΩ). Inner-filter effects from the absorbance of excitation light or emitted light before being detected by the fluorometer were corrected by using the absorbance data of each sample (Ohno 2002). Raman and Rayleigh scattering were removed following McKnight et al. (2001) and Cory et al. (2010) blanked with ultrapure Millipore DI water. The Biological Index (BIX), which is a measure of microbially-produced OM contributions, was determined from the ratio of fluorescence intensity at 380 nm to that at 430 nm and at the excitation wavelength of 310 nm from corrected EEM spectra. Fresh, unhumified and proteinaceous microbial DOM will have values between 0.8 and 1.0, whereas terrestrial, photosynthetic, and humified or humic-like DOM has values < 0.6 (Birdwell and Engel 2010). Due to CDOM photodegradation in streams and springs, surface water fluorescence indices were treated as estimates (Birdwell and Engel 2010). Indices were calculated using the StaRdom package (v1.1.21) (Pucher et al. 2019) in R (R Core Team 2021). Student’s t-tests were done to determine if there were significant differences between BIX values and site (α = 0.05).

The stable isotope ratios for oxygen (¹⁸O/¹⁶O) and hydrogen (²H or ²D/¹H) in water are expressed in delta (δ) notation as δ¹⁸O and δ²D, respectively, and were measured from filtered, unacidified water on a Picarro Cavity Ring-Down Spectrometer L2130-I Analyzer with Vaporization Module A0211 (Santa Clara, CA). Analytical uncertainty for δ¹⁸O is ± 0.02‰ and ± 0.05‰ for δ²D.

The stable isotope ratios of carbon (¹³/¹²C) and nitrogen (¹⁵N/¹⁴N), expressed as δ¹³C and δ¹⁵N, respectively, were measured from fine particulate OM (FPOM) collected by filtering water (up to 6 L or until filter clogging) onto pre-weighed and pre-ashed 47-mm diameter Whatman GF/F filters (nominal pore size = 0.7 µm). FPOM can include microbial cells, such as filamentous cells that would be floating in the water column, as described by Engel and Randall (2011), as well as microeukaryotes. Filters were dried at 60 °C for at least 48 hours, after which time the filters were cut into halves. One half of a filter was designated for δ¹³C_FPOM and placed in an HCl fumigation chamber for at least 8 hours to remove inorganic C and then dried at 60 °C for 24 hours. Halves designated for δ¹⁵N_FPOM were not fumigated. After treatments, filter halves were packed into tin capsules for isotope analysis.

Because stygofauna richness and abundance were very low at well sites, we focused on exploring patterns in isotopic values at the three terminal spring sites that had the greatest number of individuals collected and, therefore, sample masses sufficient for isotopic analyses: Eliza, Main Barton, and Old Mill springs. The number of taxa and number of individuals analyzed per site were determined by rarity and sample mass needed to obtain isotope values. Single individuals were used for larger taxa (e.g., large amphipods), but small taxa (e.g., copepods, isopods, and some amphipods) required composite samples of n = 2 to 4,050 individuals, depending on taxa. Salamander tail clips were also used for δ¹³C and δ¹⁵N analysis. Prior to analyses, tissue and tail clip samples were dried at 60 °C for 48 hours, separately ground to a fine powder, weighed, and packed into tin capsules.

The δ¹³C and δ¹⁵N for filters, invertebrates, and salamanders were measured at the University of California–Davis Stable Isotope Facility using a PDZ Europa ANCA-GSL elemental analyzer interfaced with a PDZ Europa 20-20 isotope ratio mass spectrometer (Sercon Ltd., Cheshire, U.K.). Vienna PeeDee Belemnite (VPDB) and atmospheric N₂ were used for the δ¹³C and δ¹⁵N international standards, respectively. Results are expressed in δ notation in parts per thousand (‰), using δX = [(R_sample/R_standard) – 1] * 1000, where R_sample and R_standard are the isotopic ratio (R) for the sample and standard, respectively; X represents the heavy isotope, whereby δX is the isotopic value, either enriched or depleted relative to the standard. Analytical uncertainty for δ¹³C ranged from ± 0.03‰ to ± 0.04‰ and from ± 0.03‰ to ± 0.06‰ for δ¹⁵N.

Total nucleic acids were extracted from the Sterivex filters obtained from each site using a modified protocol for Qiagen DNeasy Blood and Tissue Kits (Spens et al. 2016). The quality of genomic DNA yielded from each extraction was evaluated by using agarose gel electrophoresis, and the quantity and purity of recovered DNA was estimated by using a NanoDrop 2000 spectrophotometer (ThermoFisher Scientific). Paired-end reads (2 x300 bp) were sequenced on the Illumina MiSeq sequencing platform at the Molecular Research (MrDNA) Laboratory (Shallowater, TX) with the primers 515F (Apprill et al. 2015) and 806R (Parada et al. 2016) to target the V4 region of the 16S rRNA gene. Raw 16S rRNA gene sequences generated for this study were submitted to the NCBI Sequence Read Archive (http://www.ncbi.nlm.nih.gov) under the BioProject PRJNA1445879.

16S rRNA gene sequences were analyzed using MOTHUR (v1.44.3), following the MiSeq SOP (Kozich et al. 2013). Briefly, paired-end reads were trimmed to overlapping read length and used to generate contiguous reads using a high minimum quality score threshold (insert ≥ 24, deltaq = 6). Sequences were further screened for appropriate length, and sequences with ambiguous base calls and/or statistically erroneous homopolymers (i.e., 8+ sequential calls for the V4 region) were removed. Remaining 16S rRNA reads were aligned to the SILVA v138 database (https://www.arb-silva.de/). To remove reads that could have resulted from sequencing error, only unique and aligned reads were processed using the pre.cluster simple algorithm based on Huse et al. (2010), and chimeric sequences were identified and removed using VSEARCH (Rognes et al. 2016).

After preprocessing, bacterial and archaeal sequences were taxonomically classified using a Bayesian classifier with the SILVA database trimmed to the V4 region with the MOTHUR default cutoff (default = 80). Reads identified as eukaryote, chloroplast, mitochondria, or unknown were removed. Operational Taxonomic Units (OTUs) were clustered using the OptiClust algorithm in MOTHUR at a similarity level of 99% (0.01 distance cut-off). OTUs with less than five reads were removed, and total OTU reads were rarefied to 49,135 sequences (based on the library with the fewest number of reads) using the phyloseq package (McMurdie and Holmes 2013) in R (R Core Team 2021). Shannon diversity and Chao1 diversity indices were analyzed using the vegan package (v2.6) in R (Oksanen et al. 2025). Higher Shannon index values indicate higher diversity, whereas Chao1 estimates total richness after correcting for rare or low abundance OTUs. Alpha-diversity indices were analyzed for statistical significance using a Kruskall-Wallis nonparametric test and differences between groups were analyzed post-hoc using the Dunn test. Beta-diversity was visualized using a Bray-Curtis distance matrix and nonmetric multidimensional scaling (NMDS) in phyloseq. Environmental variables were fit to NMDS space using the env.fit function and significance was analyzed using adonis2 in vegan.

FAPROTAX (Louca et al. 2016), with manual curation of the database based on literature, was used to explore putative microbial metabolic functions of major taxonomic groups, specifically related to chemolithoautotrophic metabolic processes that have known carbon isotope fractionation factors. FAPROTAX compares taxonomic information from 16S rRNA genes in a sample to the metabolisms of cultured representatives within an assigned taxonomic group. One caveat is that poorly characterized metabolisms for some microbial groups, as well as groups with unique taxonomic identities, may be underrepresented in the metabolic profiles for some samples (Sun et al 2020; Djemiel et al. 2022). The likelihood of each metabolism identified for a sample and location was assessed based on the evidence of electron donors from water samples analyzed in our study or reported in the literature. For example, methane has been detected in saline and sulfidic wells in the past (Zhang et al. 1998), thus the presence of methylotrophs and methanogens would be likely. For reductive acetogenesis, the higher relative abundances of taxonomic groups associated with fermentation, which generates H₂ as a byproduct, were used as an indicator of electron donor variability.

We assessed the abundances of microbial 16S rRNA genes using quantitative PCR (qPCR) with iTAQ Universal SYBR Green Supermix on a CFX96 Touch Real-Time PCR Detection System (Bio-Rad, Hercules, CA). DNA was diluted to 3 to 5 ng/μL before amplification. Bacterial 16S rRNA genes were amplified using the sequencing primers and archaeal biomass was amplified using the 915F/1059R primer set (Takai and Horikoshi 2000). PCR conditions for bacteria and archaea included an initial denaturation step of 95 °C for 3 min, followed by 40 cycles of denaturation at 95 °C for 30 sec, annealing for 30 sec at 50 °C for bacteria or 53 °C for archaea, and extension at 72 °C for 30 sec. The qPCR reaction mixtures (20 μL) contained 10 μL of supermix, 200 nM of forward and reverse bacterial primers or 500 nM of forward and reverse archaeal primers, and 2 μL of DNA template. All qPCR reactions were performed in triplicate.

The qPCR data were analyzed using Bio-Rad’s CFX Manager Software (v3.1). The 10-log fold standard curves used to convert the threshold cycle (C_t) value into copy numbers ranged from 10² to 10⁷ copies. The standards were produced from Escherichia coli DNA for bacterial biomass quantification (Dhanasekaran et al. 2010) and gBlocks Gene Fragments (Integrated DNA Technologies) derived from Nitrosopumilus maritimus and prepared following manufacturing instructions for archaeal biomass quantification. Amounts of standard DNA were checked by using a NanoDrop2000 and copy numbers in each standard were calculated based on Whelan et al. (2003). Percent amplification efficiencies were calculated by the instrument manufacturer’s instructions (Bio-Rad), which ranged from 89% to 92% for bacteria and 91% to 101% for archaea. The average R² value was 0.99 for both primers. The average slopes and y-intercepts for bacterial biomass were -3.61 and 38.36, and -3.46 and 37.04 for archaeal biomass. Extraction and no-template controls (three per plate) were used to check for cross-contamination. PCR products were verified by using gel electrophoresis.

Because differences in cell size may affect DNA extraction and the potential presence of extracellular polymeric substances containing extracellular DNA could affect biomass calculations by using qPCR (e.g., Cole et al. 1993; Nisha et al. 2007; Lima et al. 2022), we conservatively estimated the total number of cells/mL in water based on the relative abundance of each genera, as described in Epp Schmidt et al. (2023). To account for differences in copy numbers for 16S rRNA genes, we used the closest relative based on rrnDB (v5.8; Stoddard et al. 2015). If the number of copies of 16S rRNA gene were unknown, then the average for the entire domain was used (i.e., 5.3 copies for bacteria and 1.7 copies for archaea).

Based on the relative abundances of bacterial and archaeal OTUs and the total biomass based on the qPCR results, in cells/mL, we estimated the proportion of each putative metabolic function identified from FAPROTAX. The cells/mL of each putative metabolic function were log transformed and analyzed using the pheatmap package in R (v1.012).

We first explored spatial patterns in geochemistry, nutrient concentrations, and FPOM isotopic composition by using principal component analysis (PCA) in the vegan package (Oksanen et al. 2025) in R (R Core Team 2021). The initial data set contained 21 nutrient and geochemical variables from four spring sites, 19 wells, and three recharge creeks (Fig. 2A, Table 2). A Pearson correlation matrix was used to exclude highly correlated variables (|r| ≥ 0.70) prior to analysis to reduce multicollinearity. Data were z-score transformed prior to analyses. The final PCA dataset included 10 variables for water quality, fluorescence indices, and ionic and stable isotope values (water temperature, DO, conductivity, NPOC, δ¹⁸O of water, TP, SO₄^2-, BIX, NO₃^-, and NH₄⁺). We then ran a second PCA restricted to the 19 locations where FPOM was collected to examine how δ¹³C_FPOM and δ¹⁵N_FPOM values co-varied with geochemistry and nutrients.

Trends in community composition and geochemistry were assessed with a redundancy analysis (RDA) performed using vegan in R. For each site, data from a single sampling event were chosen for analysis, and in most cases, sites were sampled on the same day for both geochemistry, microbial composition, and metazoan taxa. In a few instances where water and animal samples were not collected on the same day, metazoan data were selected from the closest sampling date. PCA was used to ordinate scaled, centered microbial richness data (i.e., summed reads of filtered OTUs at the order level) for each sample. Site scores for the first three PCA axes were used as environmental variables in RDA. Geochemical variables with high missingness (>10%) or no variability were removed prior to analysis, as were variables with high correlation (Pearson R > 0.75). Missing values representing less than 10% of observations were assigned mean values. Environmental variables were z-score transformed and centered prior to analysis. To minimize the outsized influence of rare taxa (Zhao et al. 2017), singletons representing less than 1% of total abundance were removed prior to analysis. In four instances where a taxon was not identified to the same taxonomic level as in other samples, those animals were assigned to lower taxonomic levels proportionally, based on their abundance in other samples from the same site. Animal count data were Hellinger-transformed prior to analysis (Legendre and Gallagher 2001). Forward selection using the ordistep was used to parameterize RDA. Variables were added until variance inflation factors reached 20. Global RDA significance and significance of individual axes was assessed using a permutation test (999 permutations) (Borcard et al. 2018).

Isotope mixing models were run using the three most abundant metabolisms identified using FAPROTAX for all sites, with the metabolic endmembers being reductive acetogenesis via the WLP, oxidation of sulfur compounds using the CBB cycle, and aerobic ammonia oxidation using the rTCA or 3-HP/4-HB cycle. All models were run with MixSIAR (v3.1.12) using a run set to Very Long, which has three Markov chains, set to a chain length of 1,000,000, with a 500,000 burn in and a 500 thinning (Stock et al. 2018). Process error was set to true and residual error was set to false. Source δ¹³C was either that of surface OM for heterotrophs or inorganic carbon for autotrophs. Source δ¹⁵N was either that of surface OM for heterotrophs or dissolved nitrate (δ¹⁵N-NO₃) from the literature or from the City of Austin records (Musgrove et al. 2016) for denitrification. δ¹⁵N-NO₃ values were used to back-calculate the δ¹⁵N of ammonia (NH₃) substrate for nitrifiers, which required assuming that groundwater nitrate was predominantly produced via microbial ammonia oxidation. This assumption was based on δ¹⁵N-NO₃ values from the Barton Springs segment (range: 4.5 to 7.5‰) being above the typical range for inorganic fertilizer-derived nitrate (-4‰ to 4‰; Kendall and Aravena 2000). As such, fertilizer input does not appear to be the dominant source of nitrate in the aquifer. Nitrification is also mechanistically plausible, as organic nitrogen can be oxidized to nitrate subaqueously when entering the saturated zone during rapid recharge (DeSimone and Howes 1998). Multiple priors were compared and are provided in Suppl. material 1: table SS1. Generalized priors assume all sources contribute equally to the FPOM. High priors were set by using the relative percent of the top three or four microbial metabolic processes that may contribute to FPOM, based on FAPROTAX identification (Djemiel et al. 2022). Moderate priors were calculated by normalizing the relative percentages by the number of microbial metabolic sources used in each model (Stock et al. 2018). Trophic fractionations were set based on highest relative abundance of genera associated with each metabolism using values reported in the literature.

The geochemistry of surface streams, springs, and wells was evaluated to understand energy gradients and conditions across the Barton Springs segment that could shape metazoan and microbial communities and food webs (Tables 2, 3). In general, most samples were calcium-bicarbonate type waters, with four saline wells being sodium-potassium-chloride type that is more typical for groundwater wells east of the FWSWI (Fig. 1C). Surface stream water had higher δ¹⁸O values that would be expected for water subjected to evaporation, in comparison with generally lower δ¹⁸O values from springs and groundwater wells that represent the long-term mean composition of water that was rapidly recharged during larger precipitation events (Table 3). As noted in previous studies, Old Mill Spring had higher specific conductance values and sulfate concentrations compared to the other springs in the Barton Springs Complex due to having more saline water contribution (Garner 2005; Mahler et al. 2006; Hunt et al. 2019). Among the wells examined, five (Dowell, Kai Frech, Sweeney, Old San Antonio 1, and United Gas) were within 1 km of the FWSWI, intersected the FWSWI, or were in the saline zone (Suppl. material 2). Of all the wells studied, those nearest the FWSWI were characterized by higher specific conductance, chloride, and sulfate concentrations, and had lower DO and higher sulfide concentrations, as were expected based on prior studies (e.g., Lambert et al. 2010). In contrast, wells influenced by the recharge zone (i.e., near the western aquifer boundary) were generally characterized by higher DO, nitrate, and alkalinity concentrations (Fig. 1C).

Based on NPOC concentrations, the amount of OM was relatively low (i.e., oligotrophic) throughout the Barton Springs segment, at less than 1.0 mg C/L for most samples except for a couple of confined freshwater and saline wells (Table 1). Creek and spring BIX values were significantly lower than BIX values from wells (p-value = 0.006) and indicated humified CDOM, which could have been produced through photodegradation. However, there were no significant differences among well types for BIX values and distance from the FWSWI. Higher BIX values > 0.8 correspond to fresh proteinaceous CDOM attributed to microbial activity, and most likely from chemolithoautotrophy (Birdwell and Engel 2009, 2010; Hutchins et al. 2016). All but one well (Gregg) had BIX values > 0.8 and six other wells had values > 1.0.

Principal components (PC) 1 and 2 cumulatively explained 54.9% of environmental variation among sites. Creeks and freshwater wells were separated from saline wells along PC1 (32.9%), which described gradients in ammonium, sulfate, total phosphorus, and DO concentrations (Fig. 2A). Kai Frech and United Gas wells exemplified solute-rich end-members, versus the oxygenated spring end-members. PC2 (22.0%) corresponded to an isotopic and nutrient gradient for creeks, springs, and some unconfined freshwater wells, with negative loadings for nitrate and BIX and positive loadings for δ¹⁸O_water, pH, and temperature. There were no significant differences among well types for δ¹³C_FPOM and distance from the FWSWI, but δ¹³C_FPOM and δ¹⁵N_FPOM loadings explained 21.6% of variation on the first two PCs (Fig. 2B). Generally, δ¹³C_FPOM values varied by 10‰ across all sample types, from -25.34‰ at the Sweeny well to -35.23‰ at the Dowell well, and δ¹⁵N_FPOM varied by about 14‰, from -8.66‰ at the Holiday Inn well to 6.91‰ at both the Old San Antonio #1 and Sweeny wells (Table 4).

The first three PC axes of a PCA on scaled microbial taxonomic orders (results described below) described 89% of the variance in the dataset and site scores were used as environmental variables in RDA. The first axis of a global RDA was significant and explained 24.4% of variance in community structure. Only two significant variables, temperature (p-value = 0.005) and microbial taxonomic order PCA axis 2 (p-value = 0.02), were retained in the final RDA model. The RDA described weak gradients in community structure, specifically higher abundances of cyclopoid copepods at sites with cooler water temperatures and microbial orders Thermoanaerobaculales and Geobacterales, which contrasts to warmer sites with aquatic mites and the Firmicutes orders Oscillospirales, Erysipelotrichales, and Salinisphaerales.

A total of 125 different samples representing separate collections at four springs and 18 wells yielded 50,460 metazoans (Table 5). Taxa included a mixture of spring benthic and stygobiontic taxa (e.g., Isopoda, Amphipoda, Gastropoda, Ostracoda) or benthic-only taxa commonly found in spring habitats (e.g., Cnidaria, Coleoptera, Diptera, Odonata, Trichoptera, Ephemeroptera, Lepidoptera, Hemiptera, Hyallela, and Platyhelminthes) (Table 5). Eliza Spring was visited more than any location (n = 36) and had the highest number of taxa recovered (45 taxa, of which 16 are stygobionts). In general, more species were retrieved from locations having < 700 µS/cm specific conductance, higher DO, and lower to no dissolved sulfide (Tables 3, 5). Macrofauna were recovered from all six sulfidic wells, with Gregg and Old San Antonio #2 wells each having seven different taxa and despite the Gregg well having DO that was below detection (Table 2). The freshwater Target and Bliss Spillar wells each had nine different taxa recovered, and 15 other wells had at least one species recovered. Despite oxygenated conditions, no metazoans were collected from the Republic of Texas well.

The most common taxa recovered from all springs and 14 wells across geochemical conditions were undetermined species of Arachnida (mites), ostracods, and copepods. Mite occurrences and abundances were also associated with higher TN, TP, and nitrate concentrations. A deeper examination of the mite collections for this study may provide more information about possible differences in mite assemblages from freshwater versus saline zone wells. Nearly 3,800 ostracods were collected from all four springs and five wells. The majority were troglomorphic but unidentified. There were tens of thousands of copepods collected, and although calanoids are rare among groundwater fauna, they were recovered from three of the springs at the Barton Springs Complex, although they could represent contamination from the surface aquatic environment. Cyclopoid copepods are grasping and raptorial feeders that consume FPOM, algae, and smaller crustaceans (Williamson and Reid 2001). Ninety-eight percent of the cyclopoids examined did not have eye spots, suggesting they are stygobionts having a hypogean origin.

Isotopic composition (δ¹³C and δ¹⁵N) of terminal spring communities at the three sites varied in the community composition and taxa available or detectable for isotopic analysis. Eliza Spring had the most taxa-rich food web (Fig. 3A) and we were able to collect isotopic data for 12 different stygobiontic and surficial taxa. Copepods, which were almost exclusively stygobiontic, were the consumer group with the lowest δ¹³C and δ¹⁵N, and surficial algivorous invertebrates (e.g., chironomid and psephenid larvae) had similarly low δ¹³C and δ¹⁵N. Stygobiontic snails (Phreatodrobia spp.) had some of the highest δ¹³C values (-29‰) and comparable δ¹⁵N values to the flatworm Dugesia sp. Salamanders had the largest range in isotopic values, but means and ranges differed between the two species; E. sosorum δ¹³C ranged from -38 to -29‰ and δ¹⁵N from 9 to 12‰, whereas E. waterlooensis had a similar δ¹³C range but δ¹⁵N was from 4 to 10‰. The δ¹⁵N values of the stygobiontic amphipod Stygobromus spp. were similar to E. sosorum. Other consumers in the Eliza assemblage (i.e., mites, annelids, isopods) had δ¹³C and δ¹⁵N values between the end-member groups. The Old Mill and Main Barton assemblages had less taxonomic richness but exhibited similar overall ranges in isotopic values (Fig. 3B, C). At Old Mill, copepods and the surficial amphipod Hyalella sp. had the lowest δ¹⁵N values, whereas annelids had the highest δ¹³C values (-25‰). At the Main Barton site, copepods had the lowest δ¹³C and δ¹⁵N values among the metazoans. However, tail clips from two E. sosorum individuals had low δ¹³C values (-38‰) similar to individuals from Eliza Spring.

We generated 4,386,276 16S rRNA gene sequence reads that represented 77 bacterial and archaeal phyla. Bacteria dominated all sample types and locations, with archaeal reads having relative abundances from 0.4% to 5% and being mostly from the phylum Thermoproteota. In general, community differences were significant when comparing sample types across the recharge to discharge regions, with freshwater and saline water wells having higher OTU-level diversity than the spring and creeks based on Shannon (Fig. 4A) and Chao1 (Fig. 4B) diversity indices. Two wells had communities with low diversity overall, with 50% of the reads from the freshwater Wyldwood well community being affiliated with the family Williamwitmanaceae (phylum Bacteriodota) and 48% of the saline Lovelady well community being comprised of members of the Nitrospirota.

An NMDS and multilevel pairwise comparison of microbial community compositions supported that the same sample types were more similar to each other than to other sample types (adonis Pr(>F) = 0.0001, R² = 0.325) (Fig. 4C). Spring and creek communities were also significantly distinct from well communities (adjusted p-value = 0.006, Suppl. material 1: table SS2), which slightly overlapped with each other in NMDS ordination space. Communities from freshwater wells corresponded to higher nitrate concentrations, creek and spring communities with higher δ^¹⁸O values, and saline and sulfidic wells with higher sulfate and ammonium but lower DO concentrations. Community composition similarities and differences were attributed to specific taxonomic groups (Fig. 5, Suppl. material 1: table SS2). For instance, despite members of the class Gammaproteobacteria (phylum Pseudomonadota) representing from 18% to 65% of sequence reads in the entire dataset, their relative abundances were significantly higher in creeks and springs than in all wells (Kruskal-Wallis p-value = 0.016), with relative abundances decreasing from freshwater to sulfidic and saline wells. Members of the order Burkholderiales were more abundant in creeks, springs, and freshwater wells, and Frankiales (phylum Actinomycetota) were more abundant in creeks. Members of the class Campylobacteria were more prevalent in springs, by contributing up to 20.5% of all sequence reads in some samples. Members of the phylum Bacteroidota (class Bacteroidia) had the highest relative abundances from creeks and freshwater wells. Members of the class Alphaproteobacteria (phylum Pseudomonadota), class Bacilli (phylum Bacillota), the class Nitrososphaeria (phylum Thermoproteota), and the class Nitrospiria (phylum Nitrospirota) also had high relative abundances from freshwater wells. In contrast, sulfidic and saline wells had distinctly different communities, with notably higher relative abundances of the order Thermodesulfovibrionales (phylum Thermodesulfobacteriota), the class Omnitrophia (phylum Omnitrophota), the class Parcubacteria, and the class Anaerolineae (phylum Chloroflexota).

We conservatively calculated the total microbial biomass that could be available as a food source for metazoans at each location. Biomass of bacterial and archaeal cells varied by up to four orders of magnitude across wells (Fig. 6A, B), ranging from 2.9 × 10⁴ to 3.4 × 10⁷ bacterial cells/mL and 9.2 × 10³ to 3.4 × 10⁷ archaeal cells/mL. The highest biomass was from the sulfidic Gregg well (3.4 × 10⁷ bacterial cells/mL and 1.9 × 10⁷ archaeal cells/mL), whereas the lowest was from two freshwater wells, Holiday Inn (9.2 × 10² archaeal cells/mL) and Circle C (2.94 × 10⁴ bacterial cells/mL). Creeks and springs had lower spatial variability, ranging from 2.4 × 10⁶ to 2.8 × 10⁷ bacterial cells/mL to 5.5 × 10³ and 1.9 × 10⁵ archaeal cells/mL (Fig. 6C, D).

According to FAPROTAX, 18 metabolic functions were identified based on the taxonomic profiles from 16S rRNA genes (Fig. 7). Using conservative interpretations (e.g., Djemiel et al. 2022), predicted metabolic functions for aerobic heterotrophy, fermentation, and anaerobic heterotrophy reflected the geochemistry across the aquifer. Based on the relative abundances of taxonomic groups and their assigned metabolisms, aerobic heterotrophy was most prevalent in the creeks (median = 81.51%), springs (median = 68.39%), and freshwater wells (median = 43.82%), but contributed less to the metabolic potential for saline and sulfidic wells (median = 29.63%). This pattern matches decreasing DO with increasing proximity to the FWSWI and anoxic conditions in the saline and sulfidic wells (Suppl. material 2). A shift away from heterotrophic metabolism for the saline and sulfidic wells coincided with greater evidence for putative chemolithoautotrophic pathways based on the relative abundances of taxonomic groups capable of acetogenesis, reduced sulfur compound oxidation (herein, simplified to sulfide oxidation), aerobic ammonia oxidation, and nitrate reduction, among other metabolisms (Fig. 7). Despite the variety of pathways, collectively, chemolithoautotrophy only explained a median of 9.51% of the metabolic potential for all freshwater wells and a median of 11.32% for saline and sulfidic wells, although there was variability. For instance, the sulfide oxidation explained as much as 52.35% of the metabolic potential in the Lovelady well, compared to 1.7% to 22.57% for the rest of the saline and sulfidic wells. There was also evidence for ammonia oxidation from anoxic wells (DO ≤ 0.2 mg/L) (Sweeney, Kai Frech, and United Gas), which represented up to 8.24% of the communities. Evidence for metabolic contributions from methanogenesis and methylotrophy was low for all sites, except at the United Gas well where methanogenesis could contribute 1.09% of the metabolic potential, and methylotrophy could contribute 1.39%. In contrast, sulfide oxidation represented a much larger fraction (22.57%) of the microbial community biomass in this well.

Although isotope mixing models suggested that at least a portion of the FPOM was derived from chemolithoautotrophy in all sites tested, most of the FPOM at springs (73.3% to 87.4%) and wells (55.3% to 88.5%) was derived from heterotrophy not chemolithoautotrophy (Table 6). These findings were supported across all three models with different priors. Sulfide oxidation coupled with nitrate reduction using CBB could have contributed about 20% or more to the δ¹³C_FPOM at the Main Barton and Eliza springs, but about 10% to δ¹³C_FPOM at Old Mill Spring. Photoautotrophy could also explain the δ¹³C_FPOM from Upper Barton Spring. For the wells, combinations of the most prevalent autotrophic pathways (WLP, CBB, and rTCA with 4-HP/4-HB cycles) could have contributed up to 44.8% of the δ¹³C_FPOM signature, which was similar to isotope mixing results from the San Antonio segment (Hutchins et al. 2016) but without assigning metabolic potential. The suboxic Dowell well had the lowest heterotrophic contribution for δ¹³C_FPOM and the highest contribution from reductive acetogenesis using WLP, likely attributed to the relatively high abundance of anaerobic taxa such as Candidatus Omnitrophus, and members of Thermoproteota and Thaumarchaeota (Table 1, Fig. 5). Several metagenome-assembled genomes have shown that members of the order Omnitrophacaea, including Ca. Omnitrophus, have genes coding for WLP (Perez‐Molphe‐Montoya et al. 2022; Seymour et al. 2023). Although aerobic ammonia oxidation and heterotrophy could explain the δ¹³C_FPOM from the Kai Frech and Bliss Spillar wells, the lowest contribution to bulk δ¹³C_FPOM for all springs and all other wells was from aerobic ammonia oxidation using the 3-HP/4-HB pathway, which would only impart a δ¹³C fractionation of -0.2 to -3.8‰ from dissolved inorganic carbon (Table 1).

We evaluated the diversity of stygofauna and microbes in the Barton Springs segment of the Edwards Aquifer, focused on characterizing microbial communities associated with chemolithoautotrophic OM production, and assessed the potential contributions of that OM in the diets of stygobionts and salamanders in the Barton Springs segment of the Edwards Aquifer. Despite a long history of biological investigations throughout the Edwards Aquifer, including its major terminal springs (e.g., Holsinger and Longley 1980; Gibson et al. 2008; Nissen et al. 2018; Krejca and Reddell 2019), several new invertebrate taxa and new taxonomic distributions were recorded from this study (Table 5). At the Barton Springs Complex, new detections were made for the thermosbaenacean Tethysbaena texana (Maguire, 1965) (which was also a new county record), the isopod Texicerberus castilloi Hutchins, Schwartz & Coleman, 2021 (Hutchins et al. 2021), the amphipod Parabogidiella americana Holsinger, 1980, the snail Phreatodrobia cf. micra, undetermined bathynellids, and an undescribed ingolfiellid. The Barton Springs groundwater snail fauna consisted of the widespread species Phreatodrobia nugax (Pilsbry & Ferriss, 1906), which has at least two distinct shell morphologies, including a spectacularly ribbed and partially uncoiled shell shape (Hershler and Longley 1986). The ribbed form, along with the rarer Phreatodrobia punctata Hershler & Longley, 1986, were only collected as empty shells during the study. At the Bliss Spillar well, shells of P. nugax (also a new record) were more depressed and similar to the snail populations common to San Marcos Springs, Texas. Five stygobiont species were recovered from the Target well, including a new county record for the amphipods Texiweckeliopsis insolita (Holsinger, 1980) and Artesia subterranea Holsinger, 1992, a northern range extension for T. insolita, as well as a new location for the recently described isopod T. castilloi (Hutchins et al. 2021). Three stygobiont taxa were recovered from the Bliss Spillar well, P. nugax, an undetermined bathynellid, and the isopod Cirolanides nr. texensis that had been previously noted from this location (Nissen et al. 2018). Lastly, like the bathynellids recovered in this study, there are many other taxa that were unidentifiable to species. There is potential for additional novel diversity and species distributions in this segment, which could have important implications for conservation and aquifer resource management.

In the San Antonio segment of the Edwards Aquifer, minimal mixing between surface recharge and confined groundwater establishes a steep hydrogeochemical gradient between oxygenated freshwater and confined, anoxic and sulfidic saline water along the FWSWI. Geochemical energy gradients along the FWSWI are exploited by chemolithoautotrophic microbial communities (Engel and Randall 2011; Gray and Engel 2013) and chemolithoautotrophic primary productivity partially supports stygobiont food webs with higher species richness along the FWSWI (Hutchins et al. 2016). For the Barton Springs segment, the overall stygobiontic distribution pattern mirrored that observed from the San Antonio segment. Stygofaunal species richness was highest at the Barton Spring Complex, followed by the Target well at the FWSWI. Twenty-two different taxa were also recovered from saline and sulfidic wells, despite the lower DO concentrations (Table 5). However, we also recovered relatively high diversity from the Bliss Spillar well from the unconfined zone near Bear Creek, and diverse taxa (e.g., unclassified aquatic mites, copepods, ostracods, and annelids) across a range of hydrogeologic conditions, including from sulfidic wells and wells far from the FWSWI (Fig. 1, Table 5). Fauna recovered from near recharge zone features could be due to being more frequently flushed into the subsurface during precipitation and recharge events. Stygofauna near recharge features also have greater access to allochthonous, photosynthetically-derived OM (e.g., Humphreys 2009; Venarsky et al. 2014; Shapouri et al. 2016; Saccò et al. 2022).

Previously in the San Antonio segment, borehole camera surveys documented microbial biofilms, with taxa linked to chemolithoautotrophic metabolisms, floating in wells and coating exposed carbonate rocks (Engel and Randall 2011). Without conducting visual surveys of wells in the Barton Springs segment, there are several lines of evidence that point to chemolithoautotrophic microbial activity in the Barton Springs segment, including BIX measurements of CDOM, well and spring microbial diversity, and δ¹³C_FPOM patterns. First, despite low NPOC concentrations, CDOM from many of the groundwater wells had high BIX values (>1.0) that indicate freshly produced OM from microbial processes (Birdwell and Engel 2010). Second, taxonomic groups with high relative abundances in the Barton Springs segment were linked to chemolithoautotrophic carbon fixation pathways, such as sulfide oxidation and ammonia oxidation, which was similar to findings from the San Antonio segment (Engel and Randall 2011; Gray and Engel 2013) and from other groundwater systems (e.g., Anantharaman et al. 2016; Probst et al. 2017; Kadnikov et al. 2020; Atencio et al. 2024). Based on the relative abundances of these and other presumed chemolithoautotrophic groups, plus biomass estimates from qPCR, chemolithoautotrophs could account for up to 24.7% of the total cellular biomass across the aquifer. This is slightly lower than the roughly 36% chemolithoautotrophic biomass calculated for a mixed carbonate/siliciclastic bedrock hillslope aquifer in central Germany (Herrmann et al. 2020) but still suggests that chemolithoautotrophic processes have the potential to contribute OM in the Barton Springs segment. Third, the δ¹³C_FPOM from wells ranged from -25.79‰ ± 0.45‰ to -34.50‰ ± 0.38‰. Chemolithoautotrophic enzymatic discrimination against ¹³C generally produces biomass with low δ¹³C values (e.g., Preuß et al. 1989; Whitaker et al. 2019; Garcia et al. 2021). However, δ¹³C_FPOM throughout the Barton Springs segment was not as depleted as was measured from the San Antonio segment, where δ¹³C_FPOM was as low as -58.2‰ (Hutchins et al. 2013; Hutchins et al. 2016). The differences between the two segments could be due to differences in the chemolithoautotrophic pathways used to fix carbon, or even due to the microbial communities having a mixture of metabolisms (Table 1).

Although a steady stream of research from groundwater and specifically karst systems is adding to our knowledge about subsurface microbiology and autotrophic metabolisms (e.g., Herrmann et al. 2020; Atencio et al. 2024), one of the challenges we faced was how to assign metabolic potential to taxa with novel identities or with poorly understood metabolisms. Even used conservatively, FAPROTAX is limited if taxonomic identity is unclassified beyond family-level. This may have affected the distribution of some metabolic processes among the taxonomic groups and affected the isotope modeling interpretations (Djemiel et al. 2022). For instance, one of the prevalent taxonomic groups from sulfidic and saline wells belonged to the genus Sulfuritalea (class Gammaproteobacteria). Some species can fix carbon via the CBB pathway, but others are mixotrophs and able to supplement carbon requirements with organic carbon when available (Kojima and Fukui 2011). Simply, recovering 16S rRNA genes for taxonomic groups capable of chemolithoautotrophic metabolism does not mean that those chemolithoautotrophic pathways are active in the aquifer. As such, and because of uncertainty, biomass estimates based on the relative abundances of Sulfuritalea were divided into CBB and heterotrophic metabolisms, which would decrease the calculated amount of isotopic fractionation attributed to CBB for this taxonomic group in the δ¹³C_FPOM mixing models. Similarly, members of the bacterial phylum Nitrospirota and archaeal phylum Thermoproteota were present in freshwater wells and sulfidic/saline wells, but their taxonomic assignments to the family-level were unclassified. Members of both phyla are capable of nitrite oxidation, complete ammonia oxidation to nitrate, as well as sulfur and hydrogen oxidation, with carbon fixation pathways being rTCA, 3-HP/4-HB, or DC/4-HB cycles (e.g., Könneke et al. 2014; Alfreider et al. 2018). Taken collectively, nitrite and/or ammonia oxidation had the potential to comprise up to 17.03% of the total cellular biomass for the well communities in the Barton Spring segment. But the δ¹³C_FPOM mixing models did not indicate high contributions of these pathways, which is likely due to the smaller isotopic fractionation factors for the carbon fixation pathways used by the taxa (Table 1). Future research that targets known genes for specific metabolic steps, accompanied by metagenomics and transcriptomics studies, can directly probe pathways, assess completeness of each pathway for different metabolic guilds, and improve estimates of specific metabolic contributions to OM production and impact to carbon isotope systematics.

The δ¹³C values of invertebrates and salamanders at the three terminal springs in the Barton Springs Complex (Eliza, Old Mill, Main Barton) were not as depleted in ¹³C compared to FPOM or stygofauna from the San Antonio segment (Hutchins et al. 2016) and reveal similar isotopically-defined food webs and basal diet subsidies (Fig. 3). Our isotope mixing models based on δ¹³C_FPOM, using established fractionation factors for microbial metabolic pathways, indicate that chemolithoautotrophy is a small contributor to δ¹³C_FPOM at all three springs. These findings differ from previous studies that suggest chemolithoautotrophically-sourced OM could be a more reliable basal energy source for stygofauna (e.g., Opsahl and Chanton 2006; Herrmann et al. 2020). Given the complexities interpreting microbial metabolisms and contributions to OM production, as well as overlapping bulk isotope compositions among consumers, other approaches used in aquatic food web studies, such as compound specific isotope analyses (Nair et al. 2021) or lipid/fatty acid profiles (Kainz et al. 2023), may provide a clearer picture of OM origin, use, and relative importance for consumers in the Barton Springs segment.

Our study followed historically low aquifer levels for the Barton Spring segment (Hunt et al. 2019) and meteorological conditions that shifted from “near normal” to “severe drought” over time (https://waterdatafortexas.org/drought). Although we did not quantify the potential OM flux into the aquifer throughout this study, decreased recharge over time during drying periods would have resulted in progressively less allochthonous OM entering the aquifer and should have pushed stygofauna to seek out other OM sources. For the San Antonio segment, Hutchins et al. (2013) link drying to shifts in δ¹³C_FPOM through time because of decreased contributions of allochthonous OM in recharge and being transported in the aquifer. Specifically, stream δ¹³C_FPOM becomes more negative with drought severity, likely reflecting reduced overland flow and a corresponding decline in C4-plant-derived OM reaching the C3-dominated streams and rivers (Bird et al. 1998; Marwick et al. 2014). In contrast, δ¹³C_FPOM in groundwater wells near the FWSWI becomes more positive with input of dissolved inorganic carbon caused by water-rock interactions, specifically carbonate dissolution (e.g., Li and Li 2018), and limited input of modern, isotopically negative δ¹³C from OM decomposition (Hunt et al. 2019). Therefore, as recharge declines, the chemolithoautotrophic OM becomes progressively more ¹³C enriched. Over time, these changes would eventually impact the isotopic compositions of stygofaunal diets as they metabolize and incorporate shifting OM sources. However, the isotopic signal recorded in the diets could also be blurred depending on how much the δ¹³C_FPOM changes as drought conditions progress and the isotopic fractionation values for the different chemolithoautotrophic pathways responsible for OM production. As such, isotopic signatures used to determine potential relative contributions of different OM sources for the groundwater ecosystem during periods of drying or drought could be misinterpreted.

In addition to the possibility that increased drying and drought affect the flux of OM into the aquifer, which could shift δ¹³C_FPOM compositions, hydrological differences between the Barton Spring and San Antonio segments of the Edwards Aquifer may also explain the weaker isotopic signal for chemolithoautotrophically-sourced OM in the Barton Springs segment. Flow velocities in the Barton Springs segment are three times higher than those in the San Antonio segment (Smith et al. 2005). Faster groundwater travel times would limit hydrologic exchange along the FWSWI and could increase the availability of surface-derived OM farther from the recharge zone. Under the low-flow conditions of this study, however, longer travel times could have increased exchange with the FWSWI, and thereby incorporate more isotopically negative FPOM produced from chemolithoautotrophy (e.g., Hutchins et al. 2016). This was not evidenced from our data likely because our study did not take place after a prolonged drought.

Hunt et al. (2019) suggest the rapid recharge rates provide long-term stability for aquifer levels and conditions during drying periods but not during extended periods of drought. Drought caused by long-term climate change that shifts precipitation patterns, elevates temperatures, and increases evapotranspiration rates (Wootten et al. 2024) will affect the quantity of water entering the aquifer and will impact the quantity and quality of allochthonous OM due to vegetation responses (Loáiciga and Schofield 2019). Overpumping the aquifer during drought also affects the quantity and quality of groundwater (Smith et al. 2005). For stygofauna in the San Antonio segment, chemolithoautotrophy provides a source of OM that is relatively well-buffered from surface conditions, overpumping freshwater, and the short-term impact of drying or the long-term impact of drought (Hutchins et al. 2016). The possibility that chemolithoautotrophic OM may play a bigger role in sustaining food webs in the Barton Springs segment during extended drought is speculative at best. Nevertheless, these findings highlight the importance of hydrological context in determining energy pathways in subterranean food webs and when considering how ecosystems may change over time due to external environmental and climatic pressures.

In summary, the results from our study suggest that the diets of stygofauna in the Barton Springs segment are primarily supported by allochthonous OM instead of chemolithoautotrophic OM. The rapid flow rates for this segment of the aquifer are likely the reason why there is limited access to and use of chemolithoautotrophically-sourced OM during normal recharge periods or even during minor drying periods. However, during extended drought, food webs may be sustained by chemolithoautotrophic OM because allochthonous OM input is reduced. At present, reliance on allochthonous OM presents potential complications for the future of the groundwater ecosystem. Karst aquifers are susceptible to contamination due to rapid recharge and flow rates (Guerra and Debbage 2021). Proposed urban development and construction projects within the City of Austin and neighboring cities, such as the expansion of transit corridors, will add impervious surfaces that reduce infiltration, speed up the rate of recharge, and add new avenues for contamination to enter the recharge and even confining zones. These immediate changes in land use adversely impact stygofauna (TxDOT 2026). Long-term studies characterizing the ecology and metabolic preferences of stygofauna, and monitoring species distribution and population sizes over time, especially during construction projects and land-use alteration, are needed to understand habitat stability and conservation needs for the stygofauna in the Barton Springs segment. Lastly, the diet resources that individual taxa rely on may be more complex than currently understood. With a lack of taxonomic resolution for many groups, combined with their small sizes, we were forced in this study to lump individuals for both taxonomic and isotopic analyses. Future work that better delineates species or closely-related species groups may provide additional insight into food web OM sources and dynamics in the Barton Springs segment.