Herpetofauna Occupancy and Community Composition along a Tidal Swamp Salinity Gradient

Occupancy patterns of herpetofauna in most tidal freshwater swamps are unknown. Tidal freshwater swamps currently face multiple threats, including salinization, which can influence their associated plant and animal communities. The impacts of salinization to herpetofauna communities in tidal freshwater swamps have not been assessed. To improve predictions regarding these herpetofauna, we conducted surveys in tidal freshwater swamps of the Savannah National Wildlife Refuge located in South Carolina, USA, from March to June, 2016 and 2017, using a variety of methods. Goals included inventorying species, determining communities, examining microhabitat associations, and modeling occupancy to predict the impacts of salinity changes. We detected 8 species of amphibians and 12 species of reptiles in our surveys. Community analyses failed to detect patterns related to measured environmental variables. Species richness and diversity declined along the salinity gradient, but the observed patterns did not match our predictions and may instead be related to site-level heterogeneity. Microhabitat associations were detected for two amphibian species via occupancy analyses. Occupancy and regression analyses indicated soil salinity may be a factor affecting nine species’ occurrences. Amphibian detections may be affected by water depth, pH values, and weather conditions. These results expand our understanding of herpetofauna within an understudied, and threatened, wetland type.


Introduction
Tidal freshwater forested wetlands (herein referred to as 'tidal swamps') are one of the wetland classes that occupy the upper reaches of estuaries, with extent directly related to rates of river discharge ). Tidal swamps occupy over 200,000 ha of the United States' coastal areas and range from Maryland to Texas (Field et al. 1991;Doyle et al. 2007). They typically occur in freshwater conditions (< 0.5 practical salinity units (psu)) from the upstream edge of tidal influence to the downstream boundary with oligohaline marsh (0.5 to 4 psu; Odum et al. 1984). Trees in tidal swamps die if exposed to chronic salinity levels of 2 psu or greater, with subsequent conversion to freshwater or brackish marsh occurring depending on the salinity regime (Conner et al. 2007;Hackney et al. 2007;Krauss et al. 2009a). The dominant plant species typically include swamp tupelo (Nyssa biflora; Walter), baldcypress (Taxodium distichum; Rich) and water tupelo (Nyssa aquatica; Linnaeus), with a variety of shrub and herbaceous understory plant species depending on salinity levels and canopy cover (Odum et al. 1984;Duberstein et al. 2014). This wetland class has been understudied until recently (Krauss et al. 2009b).
Despite recent advances, more research is needed on the function of tidal freshwater swamps as wildlife habitat. There are relatively few studies involving herpetofauna in estuaries or estuarine floodplains (Odum et al. 1984;Dunson and Seidel 1986;Kinneary 1993;Rubbo and Kiviat 1999) and no research specifically assessing herpetofauna ecology in tidal freshwater swamps. Swarth and Kiviat (2009) gave two reasons for this lack of research: tidal wetlands occur in a relatively small extent of coastal rivers, and the soft sediments in tidal wetlands make field work difficult. Odum and others' (Odum et al., 1984) foundational review of tidal freshwater wetlands listed 102 possible species of herpetofauna. However, they based their list on known geographic distributions and the assumption that herpetofauna occurring in nontidal wetlands can also occur in tidal wetlands. The Nerodia genus of snakes, and lizards as a whole, were the only herpetofauna Odum et al. (1984) specifically mentioned as occurring in tidal swamps. Their review was primarily focused on tidal freshwater marsh habitat. Another review by Odum (1988) comparing freshwater and salt marshes further identified species using tidal freshwater marshes.
It is assumed that these herpetofauna use both tidal freshwater marsh and swamps because of their proximity and connectivity, but this assumption is untested. Odum noted in both reviews that herpetofauna species richness declined from tidal freshwater marsh to salt marsh. Marsh and tidal creek surveys in a New York estuary detected low densities of herpetofauna, with few turtle and snake species and only one frog species (Rubbo and Kiviat 1999). No salamanders were detected below the mean water level, although some species were found on elevated terrain (Rubbo and Kiviat 1999). Based on an exhaustive review of the literature, we hypothesized that either tidal swamps are poor habitat for herpetofauna and relatively few occur there, or that no herpetofauna research has been published for these wetlands because of sampling difficulties.
Salinization, from sea level rise or physical changes to the environment, is a growing threat to coastal wetlands because it causes changes to the physical and chemical characteristics of the wetland environment, resulting in shifts in wetland vegetation types (Hackney et al. 2007;Herbert et al. 2015). Salinization of coastal freshwater wetlands results in a transition of freshwater habitats, including marshes and tidal freshwater swamps into oligohaline or brackish marshes, with a subsequent decrease in species diversity and changes in ecosystem function (Conner et al. 2007;Hackney et al. 2007;Wetzel and Kitchens 2007;Krauss et al. 2009a, b;Herbert et al. 2015). These impacts have been documented in numerous coastal wetlands both in the United States and abroad (Conner et al. 2007;Hackney et al. 2007;Wetzel and Kitchens 2007;Krauss et al. 2009a, b;Herbert et al. 2015). Additionally, impacts on tidal swamp microhabitat availability could result from the loss of soil stabilization provided by trees, negatively affecting wildlife species both through loss of habitat, and increased salinity. However, the impacts to the fauna of these wetlands are poorly documented, and, in some cases, the faunal communities are unknown.
The ongoing Savannah Harbor Expansion Project, located on the boundary between Georgia and South Carolina, USA, may result in elevated saltwater intrusion, altered hydrology, and decreased dissolved oxygen levels in parts of the Savannah River estuary (USACE 2012), and is the basis for this study. Freshwater diversions, flow alterations, dissolved oxygen aeration systems, and upstream mitigation land acquisitions have been proposed as management actions to mitigate the impacts of the harbor expansion on the Savannah River estuary. These management actions are predicted to result in a net loss of tidal freshwater swamps in the Savannah River estuary, although the total area of tidal freshwater wetlands is expected to remain stable as these tidal swamps are converted to tidal freshwater marsh. Additional harbor expansions are being proposed at various ports along the eastern coast of the United States as other harbors expand to meet increased shipping traffic and ship sizes that have resulted from the recent Panama Canal expansion project. Since so little is known about the fauna of tidal swamps, we do not have an adequate baseline to evaluate to what degree these management actions will successfully mitigate the impacts to wildlife in tidal swamps. Establishing a baseline for the Savannah River estuary could be a useful guide for stakeholders involved with coastal wetlands in other areas that may be impacted by harbor expansions.
We sought to address these issues by studying reptile and amphibian species in tidal freshwater swamps along a salinity gradient within the Savannah National Wildlife Refuge in South Carolina, USA. We specifically focused on reptiles and amphibians in tidal freshwater swamps because: (a) they are important components of most trophic webs (Deutschman and Peterka 1988;Regester et al. 2006); and (b) most amphibians display a biphasic life cycle with an aquatic larval stage that is sensitive to environmental changes (Rowe et al. 2003). Therefore, herpetofauna occupancy can be a useful proxy for assessing impacts on other aquatic freshwater wildlife.
We compiled a herpetofauna species inventory, assessed possible microhabitat associations, and compared herpetofauna occupancy and community composition along a salinity gradient in tidal swamps. Our hypotheses were: 1) herpetofauna species richness and diversity will decrease linearly with increasing salinity; 2) herpetofauna occurrence will decrease with increasing salinity; 3) herpetofauna richness and diversity will be greater in areas with more hummock microtopography; and 4) there are distinct communities of herpetofauna associated with changes in salinity related to proximity to the river mouth.

Site Description
The Savannah River acts as the state line between Georgia and South Carolina along its 476-km length, and its watershed is approximately 27,414 km 2 (Smock et al. 2005). The average annual flow rate, as recorded by a United States Geological Survey river gauge about 98-km north of the river mouth, is approximately 340 cubic-meters per second. The river is currently used for recreation, hydroelectric power generation, thermoelectric cooling, drinking water, and commercial shipping and navigation (Smock et al. 2005;USACE 2012). Flows are regulated to facilitate flood control and maintain water sources during drought conditions by preventing water releases once water levels drop below pre-determined levels (USACE 2012). The National Weather Service classifies the regional climate as humid subtropical, with warm, humid summers, mild winters, and few days of freezing or below freezing temperatures.
The Savannah National Wildlife Refuge is located on the Savannah River in the lower coastal plain of the South Carolina lowcountry and has 11,736 ha of freshwater marshes, tidal rivers and creeks, and bottomland hardwood forests (USACE 2012). The lower portion of the Savannah River undergoes a regular, twice a day, tidal flooding regime and is a salt-wedge type estuary (Hansen and Rattray Jr 1966). The tidal range of this river is greater than 3 m, and tidal influences persist up to 45 km upstream of the river mouth (Duberstein and Kitchens 2007). The range and consistency of the tides keep most tidal swamp soils constantly saturated, even during extended droughts (Duberstein and Kitchens 2007). Tidal flooding occurs in areas closer to the river, but the more remote areas are likely more influenced by tidal forcing of the groundwater table, which best explains their persistently saturated soil conditions (Duberstein and Kitchens 2007). Krauss et al. (2018) determined that there are about 7900 ha of tidal swamps within the Savannah National Wildlife Refuge. Approximately 400 ha of these swamps are considered saltstressed, wherein tree species reduce or suspend growth and reproduction, with noticeable effects observed at soil pore water salinities of~4 psu. The predominant tree species at the Savannah River study sites are flood-tolerant species such as water tupelo, swamp tupelo, water oak (Quercus nigra; Linneaus), and baldcypress. The predominant shrub species are hazel alder (Alnus serrulata; (Aiton) Willd.) and dwarf palmetto (Sabal minor; (Jacq.) Pers.) (Duberstein et al. 2014). Common herbaceous plants include members of the genus Eleocharis, Typha, Pontedaria, Leersia, and Sagittaria.

Field Methods
Five study sites were selected to sample two 140 ha areas within the floodplain of the Savannah National Wildlife Refuge approximately 42 river km (26 river mi) from the mouth of the river (Fig. 1). These study sites were chosen to capture the current extent of an existing tidal swamp salinity gradient and measure observed differences in herpetofauna species between the [Front] Savannah River and the Back River, one of its distributaries. Streamside study sites (i.e., Swamps 1 to 4) were next to the main channel of the Savannah River and several tidal streams (Fig. 1). The backswamp study site was adjacent to the Little Back River, a distributary of the Savannah River (Fig. 1). The streamside study sites were created along a salinity gradient to assess the impact of increasing salinity on herpetofauna occupancy and community composition. The backswamp study site was created to expand the spatial scope of our study into more remote areas of the estuary and assess differences in herpetofauna community composition and occupancy related to hydrology and salinity.
Stratified-random, 10-× 10-m (100 m 2 ) sampling grids (N = 82) were generated to collect data for herpetofaunal occurrence and environmental parameters in the streamside and backswamp study sites over the course of two field seasons. The sample grids within the streamside area were generated inside the four study sites along the salinity gradient, whereas the sample grids in the backswamp area were generated inside the backswamp study site, which was entirely freshwater (Fig.  1). Sample grids were sampled from 1 March -1 June, 2016 and 2017, to maximize species detections as herpetofauna emerged from winter hibernacula early in the season and exhibited increased activity as the season progressed. There was a risk that we could have failed to detect winter or summeractive species during this sampling period, but we feel that this risk was minimized due to our use of a variety of sampling methods and the degree of temporal overlap in their activity periods and our sampling periods.
During 2016, the sample grids (N = 52) were repeatedly sampled once per month for a total of three visits per grid. In 2017, we adjusted the sampling regime to allow longer residence times for sampling equipment, in hopes this would lead to higher detection rates (especially for cover boards). The 2017 sample grids (N = 30) were repeatedly sampled three days per month for a total of nine visits per grid. All sample grids were placed greater than 100 m from the nearest river or large tidal creek to minimize edge effects and avoid areas that might have great differences in soil composition and hydrology (e.g., high flow rates during ebb conditions).
We surveyed the 10-× 10-m sample grids for herpetofauna using multiple methods outlined by Heyer et al. (1994) and Graeter et al. (2013): area-constrained visual surveys, anuran vocalization surveys, aquatic traps (pyramid crawfish traps, Gee minnow traps, turtle hoop nets, and trash can traps), and cover boards. We used dip nets as outlined by Shaffer et al. (1994) and automated recording devices to collect supplementary information on amphibian reproduction. All detected animals were field identified to species level, captured (if possible), measured for snout-vent length and tail length (if applicable), and released. Anuran vocalization surveys lasted five minutes and were conducted at the same time as the visual encounter surveys at each sampling grid. Individual anurans were counted as 'in' if they were within 25 m of the grid and counted as 'out' if they were over 25 m from the grid. Data analyses only utilized calls that were counted as 'in'. This system prevented distant groups of chorusing frogs from being counted. Due to issues determining the ranges of anuran vocalizations detected via the automated recording devices, the recording device detections were not used for data analyses. Instead they were used to supplement the species inventory and provide information on reproductive effort within each site.
Aquatic traps were placed in suitable areas that retained standing water between tidal phases. One pyramid crayfish trap and one minnow trap were set at each 10 × 10 m grid. Hoop nets and trash can traps were set in locations that had sufficient depth and/or were located along movement corridors (e.g., tidal creeks and rivulets). The traps were checked daily, in accordance with university Institutional Animal Care and Use Committee protocols. Four 0.91-× 0.61-m cover boards were arrayed in a grid pattern at each 10-× 10-m sampling grid and were checked during the visual encounter surveys.
A total of six automated recording devices (Song Meter Model SM1, Wildlife Acoustics, Maynard, Massachusetts, USA) were randomly placed in the streamside and backswamp study sites. Five recording devices were deployed within the four salinity gradient study sites (two recorders were placed in the Swamp 1 study site due to its larger size), and the sixth recording device was deployed in the backswamp study site. The recording devices were spaced ≥800 m apart to ensure sampling independence. All automated recording devices were programmed to record daily for three minutes at the start of each hour from 8:00 P.M. to 1:00 A.M. Eastern standard time and were deployed for a minimum of ten days each month.

Environmental Covariate Data Collection Methods
Meteorological and water quality data were collected prior to the start of all surveys during both sample seasons. Air temperature, relative humidity, and average wind speed were measured using a handheld weather meter (Kestrel 3000, Kestrel Instruments, Boothwyn, Pennsylvania, USA). Water temperature, surface water salinity, pH, and conductivity were measured using a low-range combination pH, conductivity, and total dissolved solids meter (HI98129, Hanna Instruments, Smithfield, Rhode Island, USA). Water depth was measured at the deepest point in each sample grid using a folding ruler; values were rounded to the nearest half a centimeter. We measured dissolved oxygen levels prior to the start of each survey during the 2017 season using a portable dissolved oxygen meter (MW600, Milwaukee Instruments, Inc., Rocky Mount, North Carolina, USA).
We measured tree basal area during April of each sampling season using a 10-factor wedge prism, which was attached to a 63.5-cm string to ensure that the prism was held at a uniform distance from the body. Observers picked a beginning reference point and, while holding the prism at chest level and 63.5 cm away from their body, moved in a circle back to the reference point. Observers determined which trees, if any, had stems wider than the prism. Trees with stem diameters wider than the prism were counted for the basal area calculation, and the total number of trees with stems wider than the prism were multiplied by a basal area conversion factor of 10 to generate the final basal area values.
We measured tree canopy cover using during April of each sampling season using a convex spherical crown densitometer (Convex Model A, Forestry Suppliers, Inc., Jackson, Mississippi, USA), which had a reflective mirror face with 25 etched grids. These grids were visually subdivided into four sub-grids to generate a total of 100 sub-grids. Observers stood in the middle of the sample grids and used the densitometer and a compass to collect canopy cover values in four cardinal directions. To measure the amount of canopy cover, observers held the densitometer level and away from the body, and they counted the number of sub-grids on the densitometer which were filled with canopy cover. The four canopy cover values were averaged to generate one tree canopy cover value for each sampling grid.
We estimated hummock and hollow habitat microhabitat cover once per season within the 10-× 10-m sample grids by subdividing the grids into sixteen 2.5-× 2.5-m quadrats and visually estimating percent cover of hummocks and hollows within each quadrat. Hummocks were defined as areas of: raised soil topography, at least 10 cm in height that were at least 1 m 2 in total area and not covered by any trees with a diameter at breast height > 10.0 cm. The average height of most hummocks is 15 to 20 cm, so our 10-cm height included most below-average height hummocks (Duberstein and Conner 2009). Hollows were delineated as lower elevation areas, less than or equal to the base elevation of the floodplain (Duberstein and Conner 2009).
We measured soil compaction once per season using a static cone soil penetrometer (Static Cone Penetrometer Complete, AMS, Inc., American Falls, Idaho, USA). The cone of the penetrometer was inserted into the ground to a depth of 15 cm, and the pressure gauge reading at that depth was recorded. We collected four soil compaction readings at four equidistant points within each 10-× 10-m sample grid, with one measurement in each of the four cardinal directions. These four soil compaction values were averaged to generate one soil compaction value for each sample plot.
We measured soil salinity using water quality well stations, and we maintained one well station per study site. Well stations consisted of a platform stabilizing two wells that each housed an autonomous water sensor, one for monitoring aboveground (floodwater) salinity and one for monitoring belowground (soil) salinity. Water was free to move between belowground and aboveground sections of all wells. Sensors continued to operate accurately if water levels exceeded the tops of the wells, such as during a hurricane storm surge. Tops of the wells were sealed with locking caps, to which the suspension wires were attached. The locking caps help ensure better precision on water level readings between periods of interfacing with the sensors (e.g., downloading data, cleaning, calibrating). We did not remove caps to belowground sensor wells if the estuary flood stage exceeded the top of the well. This was to prevent artificial rapid exchange of belowground and aboveground water. Belowground salinities were measured with Aquatroll 200 sensors (In-Situ Inc., Fort Collins, Colorado, USA), which were suspended from the top of the well to a depth of approximately 36-cm below the soil surface. The aboveground sections of the belowground wells were made of schedule 40 PVC that did not allow water exchange. Data from the autonomous sensors were downloaded monthly, and we averaged the soil salinity values to generate one monthly value per study site. The monthly, sitelevel soil salinity values were extrapolated to each of the sample grids within the study sites, since previous exploratory work has suggested that soil salinity would not vary substantially between our plots within each study site.

Statistical Analysis
We ran single-season occupancy models using the functions within the 'unmarked' package of the 'R' statistics software (MacKenzie et al. 2002;R Core Team 2013). We modeled occupancy as a function of seven site covariates and five observation covariates (Table 1) using a logit link function. The observation covariates were evaluated separately, after which the top selected observation covariate model was combined with the site covariates to create multi-covariate candidate models. This resulted in a set of approximately 20 candidate models per species for Aikake Information Criterion (AIC; Akaike 1973; Anderson and Burnham 2002) model selection.
Models that failed to converge were discarded from the model selection and subsequent interpretation. Occupancy models were assessed via their AIC scores, with the lowest scores indicating the model that had the highest likelihood of being selected among the candidate models. We assessed the models via their AIC scores, ΔAIC values, and AIC weights. Next, we backtransformed the occupancy and detection probability estimates from the top models and calculated their 95% confidence intervals. Then, we calculated the effect sizes of the site and observation covariates for each of the top models and determined the significance and predictive power of the models by comparing their associated standard error and p values. Lastly, we calculated the estimated proportion of sites occupied for each of the species along with 90% confidence intervals.
Detections/non-detections were analyzed as a 1/0 binary response variable. Site covariates included all environmental variables listed in Table 1. Observation covariates included: Julian calendar date, survey starting time, air temperature, wind speed, and the weather condition during sampling. When necessary, the site and observation covariate data were standardized to have a mean of 0 and a standard deviation of 1 to account for the variation in measurement units.
The backswamp study site was not monitored for soil salinity changes, in contrast to the streamside sites. However, measurements collected by multi-parameter water quality sondes every 15 min (YSI model 600 XLM, Yellow Springs, OH, USA) between 19 October 2001-07 July 2003 (i.e., during an extended drought) never exceeded 1.01 psu (Duberstein, unpublished data). To remedy the lack of contemporary salinity monitoring in the backswamp zone, we extrapolated the soil salinity values for the Swamp 1 study site to the backswamp study site. These soil salinity values were most similar to those measured within the backswamp study site, and they were within 0.90 psu of the range measured from 2001 to 2003. We hypothesized that some of the observation covariates may have had quadratic relationships with the occupancy and detection probabilities. Therefore, we modeled the air temperature, date, and start time covariates with both a linear and a quadratic effect to test this assumption.
We used PC-ORD Version 6 (McCune and Mefford 2011) software as well as the 'vegan' package in the 'R' statistics software to conduct community analyses. These analyses were used to test the hypothesis that species were distinctly grouped according to one or more of our measured environmental covariates, and they included: species richness and diversity calculations, indicator species analysis, non-metric multidimensional scaling, and cluster analysis. Sample plot by species matrices and sample plot by environmental variable matrices were created as the bases for analysis. Prior to community analysis, 13 of the 82 grids (~16%) were removed due to zero detections for all species.
We evaluated site-level differences in species richness and Shannon diversity for each survey year using a one-way Analysis of Variance (ANOVA) to test for overall differences, followed by a post-hoc Tukey test for pairwise comparisons. We also used a series of a priori contrasts to evaluate how observed site-level species richness patterns along the salinity gradient differed from our predicted pattern. We calculated the mean dissimilarity between sample grids using the Bray-Curtis dissimilarity statistic. Values close to zero suggest that areas share the same species, whereas values close to one suggest that sites do not share any species (Bray and Curtis 1957). We also estimated mean Whittaker's species turnover between sample grids and study sites by dividing the gamma diversity by the mean alpha diversity of sample grids and study sites and subtracting one from the quotient. We estimated the total species richness for the areas we sampled by using the Chao, first order jackknife, second order jackknife, and bootstrap species richness estimators.
Non-metric multidimensional scaling was conducted using a Sorenson distance measure, 150 runs of the data, and a stability criterion of 0.0000001 which was achieved by running at least 50 iterations to evaluate stability. We used vector scaling of environmental covariates to assess relationships between species ordinations and environmental covariates. Cluster analyses were  conducted using a hierarchical agglomerative polythetic process, and we used the Sorenson distance measure, the flexible beta linkage method (beta value −0.25), and five group membership to define community clusters. Community clusters were compared to environmental covariates via boxplots comparing the frequency of clusters by site and select covariates. We compared community clusters to surface water salinity, tree canopy cover, pH, water depth, hummock cover, soil salinity, and tree basal area values. Indicator species analyses were grouped using the quantitative response method of Dufrene and Legendre (1997), and the analyses included a randomization test using time of day as a random number seed running 4999 permutations. We defined indicator species according to surface water salinity, tree canopy cover, pH, water depth, hummock cover, soil salinity, and tree basal area. Finally, we conducted four standard least-squares regression analyses with environmental covariates as the explanatory variables and variations of binary species-level detection data (1/0) as the response variables. We chose to do a series of analyses to ensure that we did not fail to detect any possible relationships between species-level detections and environmental covariates. Each of these regressions has implicit and explicit assumptions that are often difficult to verify, and violations of these assumptions could lead to questionable results and reduced statistical power. Thus, we ran the four variations of the least-squares regression analysis to ensure that we could detect any patterns that may be present in the data, however subtle.
The first least-squares regression analysis used the raw binary detection data as the response variable (N = 156 for the 2016 season, N = 90 for the 2017 season). The other three regressions used detection proportions as the response variable. The proportions were created by averaging the detections based on different groupings of the data. The second regression analysis focused on the variations in years and study sites by averaging monthly samples and plot-level samples to create one proportion per site (N = 5, per year). The third regression analysis focused on the variations in years, sites, and months by averaging plot-level samples (N = 15, per year). The fourth regression analysis focused on the variations in years, sites, and plots by averaging monthly samples (N = 52 for the 2016 season, N = 30 for the 2017 season). We then evaluated the statistical significance of the relationships between species detection proportion and the environmental covariates for each regression.
Detections of species varied among methods ( Table 2). We detected 5 amphibian (4 frog and 1 salamander) and 6 reptile (4 snake, 1 lizard, and 1 turtle) species, primarily via visual encounter surveys. Two reptile (1 snake, 1 turtle) species were only detected a single time during visual encounter surveys. Cover boards resulted in detection of 2 amphibian (1 frog, 1 salamander) and 5 reptile (3 snake, 2 lizard) species, with 2 reptile (both snake) species that were not detected via visual surveys. Two reptile (1 snake, 1 lizard) species were only detected a single time using this method. We detected herpetofauna during approximately 6% of all cover board checks.
Aquatic traps yielded 3 amphibian (2 frog, 1 salamander) and 5 reptile (4 snake, 1 turtle) species. One amphibian (salamander) and 2 reptile (snake, turtle) species that were detected with the aquatic traps were not detected via visual encounter surveys. One amphibian (salamander) and one reptile (snake) species were detected only once using aquatic traps. Traps were surprisingly inefficient; we only captured herpetofauna during approximately 3% of all trap attempts. Anuran vocalization surveys detected 6 frog species, including 2 that were not detected via visual encounter surveys. Dip netting did not detect any amphibian eggs, larvae, or adults despite visual encounter detections of one egg mass, one larva, and adults.

Occupancy Analysis
The single-season occupancy analyses were completed for the five species for which there were sufficient data (Table 3). Occupancy and detection probabilities of the top models varied considerably (Table 4). A variety of site covariates were selected in the top models for each species, but soil salinity was consistently selected for most species across both years ( Table 4). The top observation covariates included date and the quadratic effect of date, air temperature, weather condition, and wind speed ( Table 4). The confidence intervals of some of the estimates ranged from 0 to 1 (Table 4), indicating a large amount of uncertainty. The observation and site covariate effect sizes also varied and were largely species specific ( Table 5). The proportion of sites occupied by each species were consistently below 50% across both years.

Community Analysis
The indicator species analysis, non-metric multidimensional scaling, and cluster analysis failed to detect any trends between species assemblages and environmental covariates. The cluster and indicator species analyses found significant groupings, but post-hoc fitting of environmental variables did not result in any clear patterns. However, species richness and diversity analyses yielded interpretable results.
Observed amphibian richness varied from 1 to 6 species per site (SD = 1.77) and was greatest in the Swamp 2 and backswamp sites (Fig. 2a). One-way ANOVA failed to reject the hypothesis that there were no differences in amphibian species richness among sites. Observed reptile richness varied from 1 to 7 species per site (SD = 1.91) and also was greatest in the Swamp 2 and backswamp sites (Fig. 2b). One-way ANOVA and post-hoc Tukey tests indicated that the Swamp 3 site had significantly lower reptile species richness than the Swamp 2 site (p = 0.04). Total species richness varied from 4 to 12 species per site (SD = 3.31), with the greatest richness again observed in the Swamp 2 and backswamp sites (Fig.   2c). A one-way ANOVA and post-hoc Tukey test indicated that the total species richness in the Swamp 3 and Swamp 4 sites were significantly lower than the observed richness in the Swamp 2 and backswamp sites (p = 0.01). Contrast analyses indicated that the Swamp 4 study site had significantly lower total richness than Swamp 2 even when accounting for effects of the salinity gradient (p = 0.04). Shannon diversity values ranged from 2.14 to 7.36 (SD = 2.02), with the same trend among sites (Fig. 2d). One-way ANOVA failed to reject the hypothesis that there were no differences in Shannon diversity values among sites.
The mean Bray-Curtis dissimilarity between the 2016 sample grids was 0.84, and the mean dissimilarity between the 2017 sample grids was 0.85. The Bray-Curtis dissimilarity between study sites was 0.77 during 2016 and 0.54 during 2017. Bray-Curtis dissimilarity between sample grids and study sites generally increased along the salinity gradient, though in some instances dissimilarity values for the Swamp 2 study site and its sample grids were relatively high. Whittaker's species turnover calculations yielded average species turnover rates of 3.30 species per sample grid and 2.13 species per study site along the gradient. Turnover rates increased in response to increasing salinity values. We estimated the total species richness for the areas we sampled by using the Chao, first order jackknife, second order jackknife, and bootstrap

Standard Least Squares Regression
The results of the standard least squares regressions were dependent on which sample groupings were tested. The relationships between the environmental covariates and species detections were strongest in the analysis that focused on variations in years and study sites (N = 5). The other analyses had significant relationships between species detections and the covariates, but their relationships were weak. In total, there were 10 species (1 salamander, 4 frogs, 1 lizard, and 4 snakes) that had significant relationships with the environmental covariates (Table 6). There were 9 species (1 salamander, 4 frogs, and 4 snakes) with significant relationships to environmental covariates for regressions with no groupings and no averaged covariates. Soil salinity was the most common environmental covariate, and it was significant for 7 species (1 salamander, 2 frogs, 1 lizard, and 3 snakes). In this analysis, water depth and pH were only significant for amphibians (1 salamander, 3 frogs).
The three analyses that focused on the effects of sample groupings each had 7 species with significant relationships to environmental covariates. The number of significant species per taxa and the number of significant environmental variables per species changed with each grouping. Regressions evaluating year and site variations by averaging monthly and plot-level samples had 1 salamander, 3 frogs, and 3 snakes with significant relationships to environmental covariates. Soil and water salinity were both the most common environmental covariate, and they were significant for 4 species (3 frogs, 1 snake). Water depth and weather condition were only significant for amphibians (1 salamander, 1 frog). Regressions evaluating variations in years, sites, and months by averaging plotlevel samples had 1 salamander, 3 frogs, 1 lizard, and 2 snakes with significant relationships to environmental covariates. Soil salinity was again the most common environmental covariate, and it was significant for 4 species (1 salamander, 1 lizard, and 2 snakes). The pH level was only significant for amphibians (3 frogs). Standard least squares regressions evaluating year, site, and plot variations with averaged monthly samples had 1 salamander, 2 frogs, and 4 snakes with significant relationships. Soil salinity was again the most common environmental covariate, and was significant (α < 0.05) for 4 species (1 salamander, 1 frog, and 2 snakes). Wind speed and water

Discussion
Our results indicate that herpetofauna are present in tidal freshwater swamps of the Savannah River estuary, and: 1) herpetofauna richness and diversity decreased with increasing salinity, yet not in the pattern that we had predicted; 2) herpetofauna occurrence varied in response to increasing salinity, with different species exhibiting negative and positive responses; 3) herpetofauna richness and diversity did not appear to be influenced by hummock microtopography; and 4) there were no distinct communities of herpetofauna that could be associated with changes in salinity or other environmental variables. Our results also indicate that soil salinity may be an important driver of herpetofauna occurrence in the tidal swamps of the Savannah River estuary. Lastly, our results support trends previously reported by Odum (1988) and Rubbo and Kiviat (1999) that estuarine wetlands exhibit relatively low herpetofauna diversity, which declines in response to increasing salinity. The Swamp 2 study site exhibited the highest total species richness and diversity values along the salinity gradient (Fig. 2). This did not follow our original hypothesis which predicted that the most freshwater site along the gradient, the Swamp 1 study site, would exhibit the highest richness and diversity values. Increased salinity, feral hog rooting activity, and edge effects created by two small drainage creeks contributed to high structural diversity at Swamp 2, including some of the highest hummock microhabitat cover in our study area. In the Backswamp site, which had the second highest richness/diversity, there was high hydrological heterogeneity due to the influence of a nearby tidal creek. The structural diversity and hydrological heterogeneity at these two sites could explain why species richness and diversity did not exhibit a linear decrease along the salinity gradient as we had originally predicted (Fig. 2).
A landscape mosaic (see Angelstam 1992;Debinski et al. 2001) of marsh and forest habitats, as well as variable hydrology, may be factors driving herpetofauna richness and diversity in Fig. 2 Observed amphibian richness (a), reptile richness (b), total species richness (c), and Shannon diversity values (d) among study sites sampled during the spring of 2016 and 2017 in the tidal swamps of the Savannah National Wildlife Refuge tidal swamps. We were unable to test for the effect of structural heterogeneity because we did not collect this data during either of our sampling seasons. It may be difficult to measure this heterogeneity, given the temporal and spatial variations in structural heterogeneity exhibited by tidal freshwater vegetation (particularly herbaceous vegetation; Dusek 2003), which we observed at our sample grids during both sample seasons.
We expected hummock microhabitat cover to be a significant covariate for many species, but this was only the case for two species (Table 5). The Ring-necked Snake (Diadophis punctatus; Linnaeus 1766) was found only on hummocks, yet hummock cover was not selected as a top site covariate for this species. The lack of a statistical relationship between hummock cover and species diversity or richness may suggest that this covariate is not directly influencing herpetofauna occurrence in tidal swamps, though low sample size and statistical power cannot be dismissed.
Our occupancy analysis and standard least squares regression results suggest that soil salinity may be an important driver for the occurrence of herpetofauna in tidal swamps. Soil salinity has a longer residence time than surface water salinity and may reflect longer-term, average salinities, which may also explain why it was a more common significant covariate than surface water salinity. Increasing soil salinity can stress or kill freshwater and oligohaline wetland plants, particularly the woody vegetation of tidal swamps, and chronically elevated soil salinities can lead to shifts in wetland plant communities from tidal freshwater or oligohaline communities to brackish or saltmarsh communities. These vegetation shifts in turn alter wildlife habitat, which is probably why it was the most common statistically significant covariate in our analyses. We observed these salinity-driven vegetation shifts at a number of our sample grids during both sample seasons.
Most species detections exhibited negative relationships with increasing soil and water salinity, though responses were varied. Habitat generalist species, or those that prefer more open habitat, would be expected to exhibit a positive response to increasing salinity. Alternatively, swamp specialist species, or those which tend to occur in habitats with more woody vegetation or canopy cover, would be expected to exhibit a negative relationship with increasing salinity. Our species richness and diversity results, in addition to several species' occupancy results, appear to indicate that there is a threshold for water and soil salinity values above which both herpertofauna occupancy and species richness/ diversity rapidly declined. This threshold appears to occur around salinities of 0.15 psu or greater for both surface water and soil pore water.
We found a somewhat counterintuitive inverse relationship between amphibian detections and surface water pH. Salinity levels and water pH are interrelated and exert influences on each other (Robinson 1929). As salinity increases, pH tends to move toward neutral (Robinson 1929). However, depending on the initial chemistry, the pH may increase with increasing salinity (Robinson 1929). Natural salts present in the water act as bases, which react with hydrogen and increase pH values (Robinson 1929;Millero 1986). Thus, increases in pH may result from a corresponding increase in salinity. Large temporal spikes in salinity corresponded with spikes in pH values in our measured environmental data, supporting this trend. The pH values are also influenced by increasing temperature, exhibiting a negative relationship (Ashton and Geary 2011). This could explain why we observed stable or decreasing trends in average pH values during both seasons despite the increasing average salinity levels.
Several species were influenced by water depth. In particular, the Southern Two-Lined Salamander exhibited an increase in occupancy probability with increasing depths (Tables 3, 5). Table 6 All statistically significant (p < 0.10) standard least squares regression results (all combinations of hierarchical groupings) and their associated r 2 values for reptile and amphibian species from tidal swamps of the Savannah Wildlife Refuge during 2016 and 2017 with sufficient presence data. A '+' denotes a positive statistically significant relationship between a species' presence data and the environmental covariate, and a '-' denotes a negative relationship This is likely because the higher water levels led this species to congregate on available terrestrial habitat, thus becoming more detectable. Conversely, Green and Gray Treefrogs (Hyla chrysoscelis; Cope 1880) exhibited decreasing detections with increasing depths (Tables 3, 5). Green and Gray Treefrogs were usually detected in shallow habitats further away from flooding, which may negate tidal fluxes. Analyses indicate that Green Treefrogs were not detected past water depths of 20 cm, and Gray Treefrogs were not detected past water depths of 8 cm.
Another consideration is that increasing water depths may increase the presence of Green and Gray Treefrog larval predators (e.g., fish), which would negatively affect their occurrence. Lastly, several reptiles and amphibians exhibited relationships with increasing meteorological variables. Increasing wind speeds negatively influenced several amphibian species' detections, which is to be expected given amphibians' sensitivity to desiccation (Tables 4, 5 and 6). Weather conditions also influenced some reptile detections. Green Frogs and Ring-necked Snakes both exhibited a positive relationship with overcast weather conditions. Green Frogs and Green Anoles both exhibited a positive relationship with increasing air and water temperatures (Table 5), which is an expected pattern for ectotherms. Air temperature and water temperature were highly correlated, so the inclusion of both covariates is likely an artifact of collinearity in our data.

Conclusion
Based on our field observations and results, we concur with Swarth and Kiviat (2009) that herpetofauna are probably limited in tidal freshwater swamps and that their presence in such swamps and particular locales is related to tidal regimes. The tidal regime determines water depths, as well salinity, and these changes in salinity consequently influence other factors such as pH. All of these factors potentially influence the occurrence of herpetofauna species in the tidal freshwater swamps of the Savannah National Wildlife Refuge, though soil salinity may have the greatest influence. Species richness was lower than may have been expected if based on Odum et al.'s (1984) assumption that most herpetofauna occurring in non-tidal freshwater wetlands can also occur in tidal freshwater wetlands. Observed species richness and diversity patterns along the gradient did not fit our prediction of a gradual decline in response to increasing salinity levels, and there appears to be a threshold salinity level above which species richness and diversity rapidly declines.
Given the available evidence, we expect the Savannah Harbor Expansion Project to generate mixed impacts on the herpetofauna in the tidal swamps of the Savannah River estuary. Changes in habitat due to saltwater intrusion may exert large impacts on herpetofauna occurrence, richness, and diversity in the Savannah River estuary, though the severity of these impacts could be difficult to predict. Impacts are likely to vary spatially and temporally in response to natural feedbacks, plus future, interacting anthropogenic influences expected for coastal rivers. Freshwater flow diversions and dissolved oxygen pumps that are proposed in the current harbor expansion mitigation plans may benefit the Savannah River estuary's herpetofauna. Yet, some species may initially exhibit negative responses to the alterations in hydrology caused by the proposed flow diversions. These mitigation actions may need to be supplemented by the acquisition of inland wetland corridors, since longer-term impacts caused by sea level rise may negate the short-term benefits provided by the currently proposed mitigation actions. Sea level rise would be expected to negatively impact herpetofauna in tidal swamps unless they, and their habitat, were able to shift further inland.