Synecological response of desert spring benthic prokaryotes and macroinvertebrates to Sierra Nevada roof pendant-derived calcium

Ariel D. Friel University of Nevada, Las Vegas Khaled Pordel University of Nevada Zach Meyers Purdue University Cale O. Seymour University of Nevada, Las Vegas Nicole J. Thomas University of Nevada, Las Vegas Fred M. Phillips New Mexico Tech Jeffrey R. Knott California State University Donald W. Sada Desert Research Institute Laura Rademacher University of the Paci c Marty Frisbee Purdue University Brian P. Hedlund (  brian.hedlund@unlv.edu ) University of Nevada


Introduction
Springs emanate where groundwater encounters geologic features such as faults or contacts that direct groundwater to the surface [1] . A combination of factors determines spring hydrogeochemistry, including the geology of the recharge zone/aquifer, source geomorphology, and climate [2] . Variable spring hydrogeochemistry, physicochemically distinct microenvironments, and the often-isolated nature of springs contribute to the uniqueness and complexity of these spring-supported ecosystems [2,3] . Despite their ecological signi cance, springs are historically understudied, and the connections between spring hydrogeochemistry and ecology are not well understood. Ecohydrogeology is a new, interdisciplinary approach to spring research and characterization that emphasizes integration of geology, chemistry, and biology [2,3,4] , with an ultimate goal of improving management strategies to protect springs from threats such as groundwater extraction and climate change [4,5] . These threats are particularly pronounced for desert springs that serve as critical habitat given they are the predominant freshwater features in most desert systems [6,7] where climate change is predicted to decrease the magnitude of mountain-system recharge that supports ow in these springs [8] . Although there is a dire need for integrated ecohydrogeological studies of springs, especially within desert spring systems, few studies have explored these ecosystems within that framework.
Bacteria and archaea are important in any ecosystem due to their abundance, diversity, and critical roles in biogeochemical cycles [9,10] , yet these communities remain understudied in desert springs. Several studies have focused on bacteria and archaea in desert springs within the Cuatro Cienegas Basin (Coahuila, Mexico) that are supported by oligotrophic geothermal groundwater [11,12,13] . These studies wide range of Cretaceous granitoid rocks with a patchy distribution of Paleozoic metasedimentary and metavolcanic roof pendants and septa also present [35,36] . These roof pendants and septa were originally deposited during the Paleozoic Era as marine sediments, beneath which granitic plutons intruded during the Cretaceous [34,37] . They formed the 'roof' of the plutons and were subsequently eroded during the late Cretaceous through the Eocene, leaving only patches of the now-metamorphosed roof pendants. The granitoid rocks have an average weight percent (wt.%) CaO and Na 2 O of 3.3 and 3.55, respectively [36] . In contrast, the CaO concentration in Paleozoic roof pendant rocks such as marble and calc-hornfels is >40 wt.% and 19.0 wt.%, respectively, with Na 2 O concentrations 0.2-0.9 wt.% in these same rocks [36] .
The objective of our study was to determine whether recharge and rock-water interaction through Sierra Nevada roof pendants affects spring geochemistry and aquatic communities across multiple trophic levels in springs that emerge thousands of meters downslope from the pendants. To address this objective, we collected ~60 physicochemical parameters, 843,111 microbial 16S rRNA genes, and 2,660 individual BMIs from ten mountain-front springs in the Owens Valley that were partially recharged through the roof pendants or only through granitoid rocks and compared their geochemistry and community ecology.

Study Area and Spring Sites
The Owens Valley (~1400 m above mean sea level [masl]) is an endorheic basin in the southwestern Great Basin between the Sierra Nevada (max. elev. 4421 masl) and . Average precipitation in the Sierra Nevada is >76 cm/yr whereas precipitation in Owens Valley (13 cm/yr) and  illustrate the effect of the Sierra Nevada rain shadow [34] . Precipitation and surface runoff from the surrounding mountains collect in the Owens River and ow to Owens Lake (dry since 1927) lakebed playa. Annually, Owens Valley average temperatures range from 4℃ to 35℃, which explains its cold-arid to cold desert climatic classi cation [38] .
Ten reference springs were sampled at the foot of the eastern Sierra Nevada in the Owens Valley. These springs are broadly distributed from the northern to southern end of the western side of the valley in order to capture geologically heterogeneous recharge zones. Geochemical and microbial samples were collected during March, 2016, May, 2016and October, 2017. BMI samples and associated aquatic habitat parameters were collected from all of the springs in a single sampling trip in summer 2017. The slope of the spring emergences was calculated using ArcGIS from 1/3 arc second digital elevation models. Field notes, sampling dates, and geographic coordinates are included in Supplementary Table 1 and  Supplementary Table 2. Field Measurements and Hydrogeochemistry Sampling Physicochemical characteristics were recorded for each spring source concurrently with the collection of microbial samples, as outlined previously [39] . A YSI Professional Plus Quatro multiparameter probe (YSI, Yellow Springs, OH) was used for measuring electrical conductance (EC; µS cm −1 ), dissolved oxygen (DO; % and mg L −1 ), pH, and temperature at each spring. Speci c conductance (SPC; µS cm −1 ) was calculated by the YSI from: SpC = EC * 1.91; where 1.91 is the temperature coe cient for waters at 25 C. Aquatic habitat parameters (current velocity, water depth, wetted width, bank cover, head cover, and substrate characteristics) were recorded for each spring while BMI samples were collected. Current velocity was measured using a Marsh-McBirney Flo-Mate 2000 in-stream ow meter (Field Environmental Instruments, INC., Pittsburgh, PA). Both current velocity and water depth represent averaged values of ve measurements collected at the center of quadrats where BMIs were sampled. Wetted width represents an averaged value of ve measurements of width collected along evenly spaced transects oriented perpendicular to the spring channel. Bank cover and head cover were qualitatively estimated. Benthic substrate was characterized in percentages as silt, sand, gravel, cobble, boulder, bedrock, and muck within the BMI sampling area using the Wolman pebble count method. These benthic substrate characteristics represent an averaged value of ve measurements of percent substrate present. Quaternary ammonium and 50% ethanol (95% ethanol diluted with NanoPure DI water) were utilized on sampling equipment between locations to prevent cross-contamination during sampling trips.
Spring water was collected with a Geopump peristaltic pump (Geotech, Denver, CO) directly from the spring emergence or as close to the source as possible using Master ex platinum-cured silicone tubing (Cole-Parmer, Vernon Hills, IL) and ltered using 0.2 µm polyethersulfone membrane Sterivex-GP pressure lters (Millipore Sigma, Burlington, MA) when relevant. Filtered water samples for measurement of major ions were collected in 250 mL high-density polyethylene (HDPE) bottles that were rinsed 3 times with NanoPure deionized water (DI) and soaked in DI for a minimum of 24 hours before sampling trips. The bottles were pre-rinsed in the eld with fresh spring water 3 times and samples were refrigerated upon collection until analysis. Filtered water samples for trace elements were collected in 60 mL HDPE bottles.
Bottles for trace elements analysis were soaked in 10% nitric acid (EMD, Omni Trace) for a minimum of 72 hours, rinsed 3 times with NanoPure DI, and dried in a laminar ow hood with a HEPA air lter. 400 µL of ultra-pure nitric acid (EMD, Omni Trace Ultra) was added prior to eld work to prevent precipitation of metals. Un ltered water samples were collected for tritium ( 3 H) and radiocarbon ( 14 C) analyses, whereas ltered water samples were collected for chlorine-36 ( 36 Cl) analysis. 3 H, 14 C, and 36 Cl samples were collected with minimal headspace in 1 L HDPE bottles and the caps were sealed with electrical tape to minimize atmospheric exchange.
The major ionic composition of each spring was assessed at the New Mexico Bureau of Geology and Mineral Resources Chemistry Lab. Cations were measured using a PerkinElmer Optima 5300 DV inductively coupled plasma optical emission spectrometry (ICP-OES), according to EPA 200.7. Anions were measured using ion chromatography (IC), according to EPA 300.0. Alkalinity and hardness were measured according to EPA 310.1 and SM 2340B, respectively. Duplicates were run on every tenth sample. Trace elements were measured for most samples at MM Lab (Las Vegas, NV), but the trace element sample for one spring (Lubkin Canyon Spring 2 -IES054) was measured at the New Mexico Bureau of Geology and Mineral Resources Chemistry Lab. Trace elements were measured using an Agilent 7500 inductively coupled plasma mass spectrometry (ICP-MS), according to EPA 200.8. Tungsten was measured using an Agilent 7900 ICP-MS system, according to EPA 200.8. 3 H samples were analyzed at the Miami Tritium Lab using low-level electrolytic enrichment techniques (https://tritium.rsmas.miami.edu/), 14 C analyses were completed by the University of Arizona AMS Laboratory (https://ams.arizona.edu/radiocarbon), and 36 Cl analyses were completed at the Purdue PRIME Lab (http://www.physics.purdue.edu/primelab).

Analysis of Spring Hydrogeochemistry Data
Prior to hierarchical clustering analysis (HCA) or plotting of Piper diagrams, milliequivalents/kilogram (meq/kg) major ion data (Ca 2+ , Mg 2+ , K + , Na + , Cl − , SO 4 2− , and HCO 3 − ) were normalized by dividing each cation by the sum of all the cations in the sample and each anion by the sum of all the anions in the sample to get the nal units of % meq/kg. Values for speci c conductance and Si included in the HCA were normalized prior to the analysis. Speci c conductance was log-transformed initially, then both the speci c conductance and Si values were normalized to a value between zero and one using the R package SuperPCA version 0.3.0 [40] . To identify the hydrogeochemical groups, an HCA was performed using a Euclidean-based dissimilarity matrix of the transformed data with the R package stats version 4.0.3. Ward's Linkage was used as the linkage type. To identify major spring water types and geochemical differences between the hydrogeochemical groups, Piper diagrams were generated using Geochemist's Workbench Community Edition 15.0 (https://community.gwb.com/). Additionally, a twosample student's t-test was conducted using the R package stats on the physicochemical data (units as shown in Supplementary Table 1) to determine which speci c parameters were signi cantly different between the hydrogeochemical groups. An expanded discussion of the hydrogeochemical framework and analysis is available in another study [41] . Major ion or trace element values observed to be below the detection limit (DL) were set to (DL/2) if they were included in an analysis. Analytes with >50% of the concentrations below the DL were excluded from analyses. DLs and analytical uncertainties are included in Supplementary Table 1.

Analysis of Spring Residence Times
Spring waters were dated using a sequential multi-tracer approach using 3 H, 14 C, and 36 Cl in that order, as outlined previously [42] . First, 3 H dating was used to identify the youngest spring waters (e.g., the shortest residence times). Suitable spring water samples with 3 H activities > 0.80 TU were dated using TracerLPM [43] . A time series of atmospheric 3 H from Modesto, CA was chosen as the input function to the convolution integral since it represents the closest and most complete times series for the study area. A dispersion model was selected as the response function because it more accurately represents mixing of different groundwater owpaths in spring discharge [43] . The lower and upper age constraints in the dispersion model were kept constant at 0 and 400 years, respectively, and the dispersion parameter (DP), equal to the inverse of the Peclet number, was initialized to 0.50. Because TracerLPM minimizes error as it seeks a solution within the speci ed boundaries, the analytical error reported for each water sample by the lab was used to assess numerical uncertainty in the model. The results from TracerLPM were compared to residence times calculated using a recharge-weighted, steady-state, backward-in time tritium mixing model [42] to identify discrepancies in residence times. The raw data, analytical uncertainty, estimated residence time, and numerical uncertainty associated with 3 H dating are reported in Supplementary Table 1 and Supplementary Table 3.
Second, 14 C dating was used to date spring water samples having 3 H activities < 0.80 TU and 14 C pMC (percent modern carbon) and d 13 C values suitable for 14 C corrections. The graphical method described in a previous study [44] was used to gain insight into the geochemical reactions affecting the δ 13 C and 14 C pMC measured in the spring water samples and to inform the selection of radiocarbon correction models. 14 C corrections were completed using Netpath XL [45] by 1) inputting the measured δ 13 C and 14 C pMC of each water sample, 2) specifying the δ 13 C and 14 C pMC of carbonate in soil and bedrock, and CO 2 in the soil, and 3) testing different 14 C correction models [46,47,48,49,50] . Hand calculations were performed by allowing the carbonate endmembers to vary substantially and then seeking suitable residence times within the potential solutions. For this step, δ 13 C rec (the δ 13 C of recharge) was allowed to vary from 0‰ to -22‰, δ 13 C carb (the δ 13 C of bedrock calcite) was allowed to vary from 0‰ to -10‰. However, the δ 13 C carb values of roof-pendant marble and disseminated calcite in the granitic plutonic bedrock likely fall within the range of 0‰ to -5‰, but disseminated calcite precipitated in fractures by groundwater could be relatively modern and hence, more negative. The raw data, analytical uncertainty, estimated residence time, and numerical uncertainty associated with 14 C dating are reported in Supplementary Table 1 and  Supplementary Table 3.
Finally, a chlorine-36 chronometer ( 36 Cl/Cl chronometer) was utilized to estimate the residence times of spring water samples that weren't suitable for either 3 H or 14 C dating [42] . The 36 Cl/Cl chronometer was produced by tting a trendline between the 36 Cl/Cl ratio and the contributing groundwater residence time of springs in the study that had young residence times (estimated using 3 H) and springs that had old residence times (estimated using 14 C). The power-law trendline (r 2 = 0.76) of the chronometer was found to be: R.T. = 700,670*( 36 Cl/Cl) −1.06 where; R.T. is residence time in years and 36 Cl/Cl is the chlorine-36 ratio of the spring water sample. An envelope of uncertainty surrounds the chronometer since there is analytical uncertainty in the measured 3 H and 36 Cl/Cl compositions of the spring water and numerical uncertainty in the estimation of 3 H and 14 C residence times. Uncertainty in the 36 Cl/Cl chronometer residence times was determined by linearizing the trendline [42] . The raw data, analytical uncertainty, estimated residence time, and numerical uncertainty associated with 36 Cl dating are reported in Supplementary Table 1 and Supplementary Table 3.

Organic and Inorganic Nutrient Analysis of Spring Sediment
One set of the duplicate sediment samples collected for microbial community analysis was utilized to measure the concentration of total organic carbon (TOC), total nitrogen (TN), and total sulfur (TS) in spring sediment. Frozen sediment samples were thawed on ice and once fully thawed were dried at 60 o C in an incubator until at a constant weight. Dried sediment samples were then ground with a steel mortar and pestle and sieved through an 88-micron sieve. All equipment was rinsed with Barnstead™ Nanopure™

Microbial Field Sampling
Water microbial communities were sampled by pumping spring water (≥ 2 L total) onto replicate 0.2 µm polyethersulfone membrane Sterivex-GP pressure lters (4 lters/spring), as described in detail in the geochemistry 'Field Measurements' section. Excess water was cleared from each lter after sampling with a sterile syringe before sample storage. Benthic microbial biomass was gathered using a sterilized shovel by collecting the top > 2 cm of spring sediment or benthic microbial mat. Four benthic samples were collected in duplicate at each spring. Sampling locations were selected based on the substrate diversity (e.g., ne-grain sediment, coarse-grain sediment, or microbial mat) of each spring to maximize sampling of unique habitats at every spring. All samples were frozen immediately after sampling and kept frozen on dry ice in the eld, and then stored in a -80 ℃ freezer until DNA extraction or organic/inorganic nutrient analysis. Field notes, including descriptions of benthic sample, are detailed in Supplementary Table 2.

Sample Preparation and Sequence Processing for 16S rRNA Gene Sequencing
Microbial DNA was extracted from the samples using the FastDNA™ SPIN Kit for Soil (MP Biomedicals, Santa Ana, CA) according to the manufacturer's instructions. The V4 region of the 16S rRNA gene was ampli ed and sequenced using the updated bacterial-and archaeal-speci c 515F/806R primer set [19,51,52] . Library prep and paired-end sequencing was conducted at Argonne National Laboratory (Lemont, IL) using an Illumina MiSeq platform (2x151 bp) per the Earth Microbiome Project pipeline (https://earthmicrobiome.org/protocols-and-standards/16s/). Analysis of the 16S rRNA gene sequences was carried out using QIIME2 version 2019.4 [53] . Paired-end reads were joined, quality ltered, and assigned to amplicon sequence variants (ASVs) using DADA2 [54] . Taxonomy was assigned using a naïve-Bayesian classi er trained on the V4 region of the Silva NR99 132 alignment [55] via the q2-featureclassi er plugin [56] . Sequence variants were aligned using Mafft via the q2-alignment plugin at the default settings. Phylogenetic trees were constructed from this alignment using FastTree [57] in QIIME2 version 2019.10. ASVs classi ed as chloroplast or mitochondria were removed. BMI Field Sampling and Sample Processing BMI communities were sampled no more than 15 m from each spring source because the source represents the most stable temperature and chemical environment for these communities [58,59] . A composite BMI community sample was collected from ve subsamples at each spring. Subsamples were collected by placing 120 cm 2 quadrats and 10 cm x 12 cm 500 µm mesh D-frame netting along the same sequence of ve equally spaced transects positioned 3 m apart: center, right bank, center, left bank, and center. Benthic substrate was collected and gently agitated using a quadrat to release BMIs from the substrate and allow them to oat downstream into the netting. Samples were preserved immediately after collection in the eld using 90% ethyl alcohol. Processing, identi cation, and enumeration via microscopic examination was conducted at Rhithron Analytical Laboratory (Missoula, MT). A random subsample of at least 300 BMI individuals were identi ed and enumerated to determine the structure of BMI communities, this number of individuals has been shown to provide adequate coverage of BMI community diversity [60] . BMI taxa were identi ed to varying levels of taxonomy. Crenobiotic taxa were identi ed to the species level, whereas other aquatic insects and mites were identi ed to the genus level. Triclads (order -Tricladida), oligochaetes (subclass -Oligochaeta), ostracods (class -Ostracoda), and nematodes (phylum -Nematoda) were not identi ed any further.

α-and β-Diversity Analysis of Benthic Microbial and Macroinvertebrate Communities
For microbial samples, rarefaction curves were constructed and plotted with the R package vegan version 2.5-6 [61] . Faith's phylogenetic diversity (Faith's PD) values were calculated for microbial samples using the R package picante version 1.8 [62] . Three other α-diversity metrics (Observed, Shannon Index, and Gini-Simpson Index) were calculated for both microbial and BMI communities using the R packages phyloseq version 1.30.0 and all metrics were plotted using ggplot2 version 3.3.0 [63,64] . All α-diversity data were tested for normality with the Shapiro-Wilks test using the R package stats. Data that passed the check for normality were then tested for homoscedasticity using the Levene's Test of Equality of Variance (implemented from R package car version 3.0-10) to ensure that the assumption of homogeneity of variance necessary to run ANOVA was not violated [65] . If variance was homogenous, data were analyzed with ANOVA and a post-hoc Tukey's Honest Signi cant Difference Test (where appropriate) using the R packages stats and agricolae version 1.3-3 [66] . Metrics that failed the Shapiro-Wilks test for normality or Levene's Test of Equality of Variance were instead analyzed with a Kruskal-Wallis Test and a post-hoc Dunn Test (where appropriate) using the R packages stats and dunn.test version 1.3.5.
Three approaches were used to explore and compare the patterns in benthic microbial and BMI Bray-Curtis community dissimilarity: (1) nonmetric multidimensional scaling (NMDS) analysis, (2) principal coordinate analysis (PCoA), and (3) analysis of similarity (ANOSIM). These analyses were carried out using the R packages phyloseq and vegan. Environmental vectors were generated from major ion and trace metals data (in mol/L) and physicochemical parameters via the env t function of R package vegan. The Mantel test was used to determine if the benthic microbial and BMI Bray-Curtis community dissimilarity was associated with geographic distance. The distance between springs was calculated using the R package geodist version 0.0.7 [67] . Spearman's rank-order correlation Mantel tests (9999 permutations) were conducted on: (1) unaveraged benthic microbial samples, (2) averaged benthic microbial samples, and (3) BMI samples using the package vegan. Prior to the Mantel test, benthic microbial community data were averaged together to control for the presence of zero-distance comparisons when analyzing multiple samples from the same spring.
The taxonomic composition of the water and benthic microbial communities in the two geochemical groups (granitic and roof pendant springs) were assessed independently of each other and independently of the BMI community. Microbial taxa were agglomerated at the phylum, order, and genus level and BMI taxa at the order and genus level. Microbial ASV and BMI count data were then converted to percent relative abundance and a cutoff was applied to merge low abundance taxa. Taxonomic agglomeration and data transformation for the biological communities was carried out using phyloseq and taxonomic barplots were plotted using ggplot2. To identify differentially enriched benthic microbial taxa in granitic and roof pendant springs, a linear discriminant analysis effect size (LEfSe) analysis was conducted using the R package microbiomeMarker version 0.0.1.9000 [68] . As suggested, benthic microbial relative abundance data was normalized using a counts per million (CPM) transformation prior to conducting the LEfSe analysis [68] .
Triplicate microbial water samples were originally sequenced for two springs: Grover Anton Spring (IES 028) and Elderberry Canyon Spring (IES 033), as part of an effort to constrain the variability in microbial community structure between replicate water samples. Given the similarity of the triplicate water samples to each other, the microbial ASV counts for those samples were averaged together for each spring individually. These averaged ASV counts were utilized for most analyses, except for the rarefaction curves and Faith's PD analysis, which were run prior to averaging the triplicate water samples. Several plots were generated for the water microbes and included in the supplement with supporting text, however, the water microbial community plots were not analyzed in detail in-text due to the weak relationship observed with the geochemical groups ( Supplementary Figures 1-6).

Co-Occurrence Analysis of Microbial and BMI Communities
The relationship between the relative abundances of benthic microbes and BMIs was visualized using a correlation co-occurrence network map. First, benthic microbial and BMI taxa were taxonomically agglomerated at the genus level and converted to percent relative abundance. Taxa representing <7.5% (benthic microbes) or <5% (BMIs) across the datasets were then merged. Next, Spearman's rank-order correlation coe cients and p-values were calculated between the community relative abundance matrices using the R package Hmisc version 4-4.0 [69] . Correlation co-occurrence network maps were constructed using the R package network version 1.16.1 [70,71] , and then plotted using the R package ggnet version 0.1.0 [72] . Genera were color-coded on the maps by BMI stress tolerance/microbial metabolism or BMI behavior/microbial attachment. Stress tolerance values (range: 0-10) were assigned to BMI taxa as previously characterized with intolerant taxa denoted by low stress tolerance values and tolerant taxa denoted by high stress tolerance values [73] .

Paleozoic Roof Pendants Impart Subtle Geochemical Differences in Similar Springs
The ten springs along the foot of the eastern Sierra Nevada provide a broad geographic and geologic representation of springs along the west side of Owens Valley ( Figure 1A). These rheocrene type springs [3] were physicochemically similar, being characterized by cool, fresh water (temperature: 11.2-18.7 ℃; speci c conductance: 66.1-470.7 µS/cm), circumneutral pH (pH: 6. 94-8.16 An HCA incorporating major ions, Si, and speci c conductance data divided spring waters into the two major geochemical groups ( Figure 1B). Student's t-tests were utilized to assess speci c geochemical differences between the groups for ~60 physical and chemical analytes. This revealed that the two spring groups possessed signi cantly different Ca 2+ /Na + molar ratios (p-value = 3.1 * 10 −5 ), Ca 2+ /Mg 2+ molar ratios (p-value = 0.05), divalent to monovalent cation ion ratios (D/M; p-value = 0.001), and slopes of the spring emergences (p-value = 0.004) (Supplementary Table 1). Spring geochemistry in both groups re ected recharge through granitoid rock, which predominates the high-elevation recharge zones in the Sierra Nevada. All spring waters also contained an excess of calcium (Ca 2+ /Na + > 0.6) over that expected for plagioclase weathering alone (Ca 2+ /Na + of 0.4-0.6; Figure 1C), consistent with dissolution of disseminated calcite during recharge [41,74,75,76,77] . However, only one spring group contained an "excess of excess calcium", with signi cantly higher Ca 2+ /Na + ratios (Ca 2+ /Na + > 1.5). This "excess of excess calcium" signature has recently been interpreted to re ect roof pendant weathering [41] , which is consistent with the location of the ve springs with these signatures. Two emerge downgradient from the Pine Creek Pendant units near Mt. Tom (Elderberry Canyon spring, IES-033) and Wheeler Crest (Wells Meadow B spring, IES-038), whereas the three remaining springs emerge downgradient of unnamed roof pendants near Lookout Point (Supplementary Figure 8). For simplicity, the spring group with "excess of excess calcium" signatures is hereafter referred to as "roof pendant springs", whereas those with lower Ca 2+ /Na + ratios are hereafter referred to as "granitic springs".

Roof Pendant Springs Host More Diverse Benthic Macroinvertebrate Communities
Following quality ltering, 843,111 total 16S rRNA gene sequences from bacteria and archaea were recovered from 39 benthic samples, with an average of 21,618 reads per sample. A total of 14,523 unique ASVs were recovered, approximating the species richness of prokaryotes in these samples, and accounting for 52 bacterial phyla and nine archaeal phyla. Rarefaction curves for all benthic microbial community samples reached an asymptote between 10,000-15,000 sequences, indicating adequate sequencing depth for robust comparisons among springs and spring groups (Supplementary Figure 9). A total of 2,660 individual BMIs were identi ed across ten composite benthic samples with an average of 266 BMIs per sample. From the BMI community data, four phyla, seven classes, 16 orders, 47 families, and 85 genera were identi ed across the ten springs.
Several α-diversity metrics were calculated to determine whether geochemical grouping in uenced diversity ( Figure 2; Supplementary Tables 5-7). No signi cant differences in benthic microbial diversity were observed between roof pendant and granitic springs ( Table S8).
Covariates of community dissimilarity (p-value ≤ 0.05) identi ed among the ~60 physical and chemical analytes were tted to the NMDS ordinations to visualize correlates of community dissimilarity ( Figure 3). Ca 2+/ Na + , Ca 2+ /Mg 2+ , and D/M ratios correlated with microbial and BMI communities from roof pendant springs, re ecting weathering of the roof pendants during recharge. Additional vectors corresponding to spring or sediment characteristics and landscape placement were also signi cant. In roof pendant springs, current velocity correlated with both communities, whereas slopes of the spring emergences and cobble percentage correlated only with microbial communities and gravel percentage only with BMI communities. In granitic springs, spring water residence time and muck percentage correlated with both communities.
LEfSe analysis of benthic microbial communities revealed 408 taxa (8 phyla, 33 classes, 69 orders, 111 families, and 187 genera) that were signi cantly (p-value ≤ 0.05) enriched in either roof pendant springs or granitic springs (Figure 4). Several phyla and one class with multiple child taxa were enriched in roof pendant springs: phyla Acidobacteria (four classes, three orders), Planctomycetes (three classes, four orders), Gemmatimonadetes (three classes, three orders), and Actinobacteria (one class, four orders), and the class Alphaproteobacteria (six orders) ( Figure 4A). In granitic springs, only the archaeal phylum Euryarchaeota (two classes, three orders) had multiple classes or orders enriched. Granitic springs were also enriched with several bacterial classes with fewer child taxa, namely Clostridia (one order), Endomicrobia (one order), and Gracilibacteria (one order), and bacterial orders Bacteroidales, Beggiatoales, Chromatiales, Desulfobacterales, Desulfuromonadales, and Methylococcales.
The enriched families and genera suggest a unique set of physiological traits associated with each geochemical group ( Figure 4B-C). Roof pendant springs were enriched with aerobic/facultatively anaerobic genera predominantly from Alphaproteobacteria (eight families, 11 genera), Actinobacteria (eight families, eight genera), Planctomycetes (three families, two genera), and Acidobacteria (two families, two genera). These taxa included those known to produce prosthecae/stalks and holdfasts in the Alphaproteobacteria and Planctomycetes, such as Bauldia, Hirschia, Gemmata, Pedomicrobium, Pir4 lineage (family: Pirellulaceae), Schlesneriaceae, Planctomyces, Pedomicrobium, Hyphomicrobiaceae, Micropepsaceae, SWB02 (family: Hyphomonadaceae), and Rhodomicrobium. Granitic springs were instead enriched with anaerobes, including fermenters in the Clostridia (one family, four genera), sulfate reducers in the Deltaproteobacteria ( ve families, eight genera), and methanogenic archaea within the Euryarchaeota (three families, three genera). Aerobic bacteria included methylotrophic bacteria in the Gammaproteobacteria (six families, four genera) and sul de oxidizers in the Beggiatoaceae ( Figure 4C). A full list of signi cant results from the LEfSe analysis can be found in the supplement (Supplementary Table 9). BMI communities differ in composition at multiple taxonomic levels (Supplementary Figure 14; Figure 5). At the genus level, communities in roof pendant springs were relatively even and dominated by vagile insect taxa, such as Lepidostoma, Optioservus, Gumaga, Malenka, Ironodes, Baetis, Argia, and Enchytraeidae with additional populations of lower abundance taxa. These taxa are primarily shredders, collector-gatherers, and predators. Pyrgulopsis was present in only one roof pendant spring, Elderberry Canyon Spring (IES-033), where it was the most abundant taxon. In contrast, granitic spring BMI communities were dominated by Pyrgulopsis and stress-tolerant, non-insect taxa, mainly Pisidium and Hyalella, with less abundant populations of the insect taxa Argia, Optioservus, and Tricladida.

Discussion
We documented distinct communities of benthic microorganisms and macroinvertebrates in springs that recharge through Paleozoic roof pendants thousands of meters upslope from the spring emergences.
Weathering of these roof pendants imparts an "excess of excess calcium" and an increased proportion of dissolved calcium compared to other ions [41] that may promote initial microbial attachment to the negatively charged granitic substrate and enhance bio lm maturation within the spring benthos. Bio lms are multispecies surface-attached microbial communities embedded in a three-dimensional matrix of endogenously produced EPS composed of polysaccharides, proteins, extracellular DNA, and lipids [78] .
Bio lm formation begins with the development of a surface coating comprised of biological material from the surrounding environment, known as a conditioning lm, that enhances the initial adhesion of microbial cells [79] . Divalent cations are known to reduce electrostatic repulsion by forming cross-bridges between negatively charged surfaces and negatively charged cell surface structures, such as phospholipids, glycolipids, and teichoic acids [80,81] . The presence of divalent cations can also promote the expression and functionality of speci c adhesins in a variety of bacterial species [80,81] . During bio lm maturation, the formation of non-speci c cation cross-bridges between cells and EPS enhances the thickness and structural stability of the bio lm and traps nutrients [81,82] . Most research on the role of Ca 2+ or other divalent cations in microbial attachment and bio lm formation focuses on ionic strength, yet neither Ca 2+ concentration nor total ionic strength accounted for the community structures observed in the roof pendant springs. Instead, these patterns were attributed to differences in valence, speci cally Ca 2+ /Na + ratios.
A past study showed that Ca 2+ /Na + ratios >1.0 resulted in larger microbial ocs (aggregations) in activated sludge compared to those with ratios <1.0 [83] . Although our study focuses on benthic bio lms rather than suspended aggregates, we also identi ed Ca 2+ /Na + ratios as the strongest correlate of community structure and propose that high Ca 2+ /Na + ratios promote microbial attachment to granite in these fast-owing springs. Only springs in the roof pendant group possessed Ca 2+ /Na + ratios above 1.0 (average: 1.5 ± 0.12), and these springs were enriched with unrelated lineages of bacteria often associated with bio lms, including members of the phyla Acidobacteria, Planctomycetes, Gemmatimonadetes, and Actinobacteria, and the class Alphaproteobacteria [84] (Figure 4). Several of the enriched Alphaproteobacteria and Planctomycetes also produce prosthecae and/or stalks, many of these structures are tipped with holdfast adhesins, which are secreted polysaccharides used for attachment and bio lm formation [85,86,87] . Previous work with a puri ed polysaccharide holdfast adhesin (fr2ps) from Hyphomonas MHS-3 demonstrated a critical role for Ca 2+ cross-bridging in the initial binding of this adhesin to germanium oxide surfaces [88] . We postulate that the higher abundance of bio lm-associated and prosthecate/stalked bacteria in roof pendant springs is a response to the molar excess of Ca 2+ , which may promote microbial surface attachment and bio lm development. In owing aquatic systems such as spring runs or streams, surface attachment and bio lm formation are necessary for benthic microorganisms as it prevents cells from being ushed out of the system. As cells attach and settle in the system, endogenous nutrient cycling is increased, promoting the development of metabolically diverse assemblages and enhances trophic interactions [84,89] .
Calcium may also be important for structuring BMI communities, as observed in other aquatic systems [90,91] , although most studies do not address mechanisms for this relationship. We propose that Ca 2+ may structure these BMI communities indirectly through trophic interactions with microbial bio lms (as discussed below) or directly by in uencing BMI physiology. For example, Ca 2+ is present at a 1:1 molar ratio with phosphoserines in the major silk adhesin of caddis y larvae, where it forms stabilizing cation cross-bridges between antiparallel β-sheets [92] . This adhesin is initially secreted from the silk gland as a uid that is then solidi ed and stabilized into metallo bers by Ca 2+ absorbed from the surrounding environment [92] . When Ca 2+ is replaced with Na + , the β-sheets in the stabilized silk adhesin revert to random coils, destroying the structure and adhesive properties of the silk [93] . Ca 2+ -stabilized silk adhesins are necessary for substrate attachment and construction of composite cases comprised of stones and detritus that are utilized by the larvae for protection and in food gathering [94] . In our study, the caddis ies Gumaga and Lepidostoma were dominant BMI taxa in roof pendant springs but were nearly absent from granitic springs, which may be a direct response to Ca 2+ /Na + ratios. Given their roles as shredders, Gumaga, Lepidostoma, and other enriched taxa, such as the stone y Malenka, likely interact extensively for the synergistic degradation of allochthonous organic matter, explaining the well-connected microbe-BMI network in roof pendant springs ( Figure 6).
Allochthonous organic matter is the dominant source of organic carbon in most springs and streams [95,96,97] and is typically derived from the debris of riparian vegetation, collectively referred to as CPOM (>1 mm in particle size) [32] . BMI shredders and microorganisms are responsible for the synergistic breakdown of CPOM into FPOM (<1 mm in particle size) that is then ingested by BMI collector-gatherers [32,97] . Both shredders and collector-gatherers preferentially feed on organic matter that has been colonized by microbial bio lms [32,98] due to the stoichiometrically favorable nutrient ratios (C:N and C:P) of microorganisms over plant biomass [32,97] . EPS, which often comprise more biomass than cells in bio lms, are an additional food source for BMIs, further facilitating the transfer of energy to higher trophic levels [99,100] . Thus, we propose that Ca 2+ -enriched populations of bio lm-associated and stalked/holdfast bacteria in roof pendant springs exert a bottom-up control on shredders and collectorgatherers as a food source and metabolic partner. Subsequently, shredders may enforce top-down controls on CPOM-associated bio lms by increasing CPOM surface area and the availability of microbial attachment sites [32,97,101] or potentially through selective grazing of larger microbial cells/ laments in the bio lm.
An alternative top-down control on benthic microbial bio lms could potentially be driven by the hydrobiid gastropod genus Pyrgulopsis, given their high population density in all granitic springs and their feeding behavior as scrapers. Our results, however, do not support this hypothesis. Only the roof pendant Elderberry Canyon Spring (IES-033) hosted Pyrgulopsis ( Figure 5). Despite the presence of Pyrgulopsis, the microbial and BMI communities in the Elderberry Canyon Spring clustered with the other roof pendant spring communities on the NMDS plots ( Figure 3). Furthermore, the microbial communities in springs with or without Pyrgulopsis were not signi cantly different (ANOSIM: R = 0.06, p-value = 0.12). Additionally, Elderberry Canyon Spring hosts enriched populations of bio lm-associated and stalked/holdfast bacteria, like other roof pendant springs. Altogether, we did not identify any evidence consistent with top-down control of bio lm-associated and stalked/holdfast bacteria through grazing by Pyrgulopsis.
In aquatic environments, current velocity is one of the most important factors determining the grain size of benthic substrate [102] , along with microbial [103] and BMI [26,27] communities. Fine grain sizes that dominate in low-ow springs slows oxygen diffusion, selecting for anaerobic microorganisms and spatially controlled redox cycles of C, N, S, and other elements [104] . The roof pendant-recharged springs tended to have higher currently velocities and were appropriately enriched with aerobic and facultatively anaerobic bacteria. This spring group also contained diverse, low-moderate stress tolerant insect taxa, in line with the coarser substrates found in these springs. In contrast, the generally slower owing graniterecharged springs were enriched with anaerobes such as obligate fermenters, sulfate reducers, and methanogens, as well as the corresponding chemolithotrophs-sulfur oxidizers and methylotrophs ( Figure 4). The granite-recharged springs also hosted highly uneven BMI communities dominated by moderate-high stress tolerant non-insect taxa. Despite these trends suggesting the importance of current velocity, there were a couple exceptions. The current velocity at South Harry Birch Spring (IES-043) was among the lowest of any of the springs, whereas Unnamed Spring north of Red Mountain (IES-029) was among the fastest. Even with their exceptional current velocities, both microbial and BMI communities in these springs clustered with their respective geochemical group (Figure 3), underscoring the hypothesis that the roof pendant-driven molar excesses of Ca 2+ is the primary factor structuring these communities.
A third alternative environmental control was spring water residence time, or the time elapsed between recharge and discharge. In general, older residence times are associated with warmer, more geochemically evolved water [105,106] . The speci c physicochemical effects associated with increased residence time are highlighted as important factors structuring microbial communities in groundwater-fed ecosystems, such as streams [16] . and aquifers [107,108] . Additionally, temperature and salinity are wellrecognized drivers of microbial community ecology across multiple biomes [18,19,20] . In our study, spring residence time estimates varied from ~40-200 years for most springs to over 1,200 years for Lubken Canyon Spring 2 (IES-054) and Spring along Hogback Ck. A (IES-024). Residence time correlated with both communities and trended with vectors for speci c conductivity and many individual ions, diagnostic of increased water-rock interactions ( Figure 3). These relationships, however, were subordinate to the geochemical groupings. Ultimately, although residence time appears to play some role in structuring the ecology of these communities, the molar excess of Ca 2+ imparted by the roof pendant remains the dominant factor.

Conclusions
This study answers the call for implementation of an integrated ecohydrogeology approach to better understand spring ecosystems [2,3,4] . By using this approach, we demonstrate, for the rst time, a strong ecological response of desert spring benthic communities across multiple trophic levels to subtle changes in spring geochemistry imparted by recharge through roof pendants. This response is striking given the overall geochemical similarity of these springs, as the recharge path for both groups is predominantly through granitoid rock. The tightly connected community responses to roof pendantderived Ca 2+ /Na + ratios-namely enhanced benthic bio lms and abundant, diverse populations of shredder and collector-gatherer insect taxa in roof pendant springs-likely increase the ux of biomass and energy through the benthos and beyond to higher trophic levels in these springs as well as associated riparian areas. These and other ecological responses to subtle differences in spring geochemistry should be further explored in the Sierra Nevada and beyond, preferably within an ecohydrogeological framework, to re ne our understanding of the connections between spring hydrogeochemistry and ecology.  Benthic macroinvertebrate diversity, but not microbial diversity, is higher in roof pendant springs. Boxplots of A) Observed sequence variants, B) Gini-Simpson Index, C) Shannon Diversity Index, and D) Faith's phylogenetic diversity (Faith's PD) values for benthic microbial communities. Boxplots were also made of E) Observed counts, F) Gini-Simpson Index, and G) Shannon Diversity Index for BMI communities. All boxes are colored according to geochemical grouping: roof pendant springs (blue) and granitic springs (red). IQRs, medians, and data points are shown. For the benthic microbial community, no statistically signi cant differences were observed between the hydrogeochemical groups. For the BMI community, richness and diversity were higher (p-value < 0.05) in roof pendant springs.

Figure 3
Both microbes and BMIs roof pendant spring communities are unique. NMDS ordination of Bray-Curtis dissimilarity between ten springs for both A) benthic microbes and B) BMIs. Each point is colored by geochemical grouping: roof pendant springs (blue) and granitic springs (red). Physicochemical factors deemed to be signi cantly (p-value < 0.05) correlated with community dissimilarity were overlaid on the NMDS plot. Names of some environmental parameters have been abbreviated on the plot: current velocity (CV), wetted width of the spring run (WW), the divalent to monovalent cation ratio (D/M), residence time (RT), and speci c conductance (SPC). The in uence of hydrogeochemical grouping on βdiversity was found to be signi cant via ANOSIM for both benthic microbial communities (R = 0.28, pvalue = 0.001) and BMI communities (R = 0.65, p-value = 0.009).

Figure 4
LEfSe analysis reveals physiological differences between microbial communities in the spring groups.
Linear discriminant analysis (LDA) effect size (LEfSe) A) cladogram and B+C) barplots depicting microbial taxa enriched (p-value < 0.05) in either B) roof pendant springs (blue) or C) granitic springs (red). Results shown in the barplots are plotted at the family and genus level. Relevant physiological and ecological traits for families and genera are depicted with colored triangles. Taxa with alpha-numeric names that had no inferred physiology/ecology and unclassi ed taxa were removed at all taxonomic levels for both the cladogram and barplots; the full list of signi cantly enriched taxa can be found in the supplement (Supplementary Table 9).

Figure 5
Benthic macroinvertebrate communities are distinct between geochemical groups. Relative abundance of BMI communities identi ed microscopically. Roof pendant springs host diverse insect taxa, whereas granitic springs are dominated by Pyrgulopsis. Taxa calculated to be <5% across the dataset were agglomerated and plotted as such. Most taxa were identi ed at the genus level, with the exception of Enchytraeidae (family) and Tricladida (order).

Figure 6
High connectivity of a BMI-microbe network in roof pendant springs. A co-occurrence network map depicting the positive correlations identi ed following the Spearman's rank-order correlation analysis performed on relative abundance data from benthic macroinvertebrates (circles) and microbes (triangles) in the ten springs. Nodes are colored by A) BMI functional feeding group (FFG) and microbial attachment or B) BMI stress tolerance and microbial metabolism. Only signi cant (p-value < 0.05) comparisons and nodes with more than one edge are displayed. Major network clusters have been marked with numbers that have been colored to demonstrate the corresponding geochemical grouping.

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download.