Comparisons of Twelve Freshwater Mussel Bed Assemblages Quantitatively Sampled at a 15-year Interval in the Buffalo National River, Arkansas, USA

: Historically, 23 freshwater mussel species have been documented from the Buffalo National River (BNR), a 246 km, free-flowing river in northern Arkansas. The potential threats to BNR include land use/land cover changes, eutrophication, recreation, physical habitat changes, and various climate change-related effects. Twelve quantitative mussel bed sites were established and then sampled using a stratified random sampling protocol to evaluate the long-term changes between 2006 and 2020–2021 in population and assemblage characteristics. We compared (1) over-all mussel bed persistence, sampling confidence levels and study-wide relative abundances, and compared species’ size and size-frequency distributions; (2) 10 overall site assemblage variables using paired t -tests; (3) site-level mean density, richness, and diversity indices using pair-wise Mann–Whitney U statistics; and (4) assemblage composition using Non-Metric Multidimensional Scaling. The major findings included the following: (1) sampling efforts based on a targeted 80% confidence level appears relatively robust, (2) BNR mussel assemblage composition and structure were relatively stable (however, small mussel bed persistence is a concern), (3) 7 of 23 sites were outliers based on freshwater mussel composition and habitat characteristics, and (4) assemblage composition changed with three species declining ( Actinonaias ligamentina , Lasmigona costata , and Ptychobranchus occidentalis ) and four species increasing ( Cambarunio hesperus , Cyclonaias tuberculata , Eurynia dilatata , and Venustaconcha pleasii ) between monitoring events.


Introduction
Freshwater mussels (Unionida) are some of the most imperiled taxa in the world [1][2][3], with the largest number of recent extinctions of any taxonomic group [4,5].North America retains the greatest unionid diversity with approximately 330 recognized species; however, more than 70% are considered endangered, threatened, or of special concern [3,6,7].Since the late 20th century, this faunal group has received a concerted effort from federal and state natural resource agencies to compile species status information [2,8,9].Conservation efforts are hampered by our understanding of this faunal group in general and the diffuse and chronic underlying impacts that have been difficult to identify and remediate.Knowledge gaps exist in mussel systematics, life-history characteristics, fish-host relationships, and the mechanisms that drive greater sensitivity to stressors and pollutants than other organisms [1,3,10].
Ecologies 2024, 5 2 The most cited reasons for mussel imperilment are impoundments, invasive species, pollution, habitat destruction, and sedimentation [1,[11][12][13].Extinction gradients are observed downstream of large-scale impoundments due to altered temperature regimes, flow seasonality, food availability, and fish-host migration, which impedes the dispersal of glochidia [13].These factors, combined with drought from global climate change, have created smaller isolated populations susceptible to genetic bottlenecks [14].Extinction gradients are observed downstream of large-scale impoundments due to altered temperature regimes, flow seasonality, food availability, fish-host migration, and sediment deposition [13,15].Declines in mussel abundance and species richness may permanently alter ecosystem functioning and resilience to disturbances and environmental change [16][17][18].Long-lived organisms like unionids tend to have reduced inter-annual variation, which makes them ideal for assessing long-term ecosystem health.Identifying shifts in mussel assemblage composition and abundance over time will provide insight into the temporal scales of habitat degradation and potential environmental factors limiting population growth within a river system.Biodiversity losses and shifts in species dominance patterns include the decline in the abundance of common species [19].Dominant mussel species are drivers of ecosystem processes and food web dynamics through their role in nutrient cycling and primary production [20].Maintaining mussel communities will influence a cascade of ecosystem services throughout the trophic web.Few studies have documented the long-term quantitative changes in mussel communities [20,21], especially those of the "least disturbed" or minimally impacted reference sites such as designated wild and scenic or extraordinary resource waters.
The Buffalo National River (BNR) is a 246 km long, free-flowing river in northern Arkansas, which flows through the Boston Mountains and Springfield and Salem plateaus to its confluence with the White River near Buffalo City, Arkansas (Figure 1).The BNR was established as the first National River in 1972 (Public Law 92-237) and later designated as both a Wild and Scenic (for a headwaters segment; Arkansas Wild and Scenic Rivers Act of 1992) and an Extraordinary Resource Water [22,23].However, only 11% of the BNR watershed is managed by the National Park Service (NPS) [21] and the watershed has been affected by a variety of threats [22,[24][25][26][27].In the early 1900s, BNR mussels were commercially harvested, and livestock were able to access the riverbed for stream crossing and drinking water, both of which resulted in population alteration and water quality reduction.More recently, sections of the mainstem and tributaries in the BNR watershed were 303 d listed as impaired under the Clean Water Act for bacteria and low dissolved oxygen [24], with increasing temperatures [28,29]), nutrients (in particular inorganic nitrogen), erosion/sediments, and trash identified as pollutants of concern [25].The National Park Service [29] recently reported increased visitor use based on entry-point monitoring and anecdotal evidence of watercraft use.
Freshwater mussel assemblages in the BNR were first documented by Meek and Clark [27] in which they identified 22 species across 26 sites and reported the mussel fauna as "neither large nor plentiful" compared to that of the White River.The BNR was revisited by Harris [30] and Matthews [31], as summarized in Matthews et al. [32], and they evaluated the species composition and the physical locations of mussel assemblages.Harris [30] recorded two additional species (Pleurobema sintoxia and Eurynia dilatata) but did not observe three species (Lampsilis siliquoidea, Potamilus purpuratus, and Ligumia recta) previously listed by Meek and Clark [27].Matthews et al. [32] expanded the survey area and distribution of freshwater mussel assemblages by approximately 33%, quantitatively sampled 23 assemblages, documented 23 species observed, and initiated a long-term monitoring program at 12 sites across three of the four upstream-to-downstream assemblage types observed in their study.Matthews et al. [32] determined 78% of the 23 species identified were state heritage ranked S1-S3 (including five of the six most abundant species) and considered the BNR a mussel diversity refuge in the Ozark Highlands.The definitions for status rank categories follow Master et al. [33]: S1, Critically Imperiled-at very high risk of extirpation in the jurisdiction due to very restricted range, very few populations or occur-rences, very steep declines, severe threats, or other factors; S2, Imperiled-at high risk of extirpation in the jurisdiction due to restricted range, few populations or occurrences, steep declines, severe threats, or other factors; S3, Vulnerable-at moderate risk of extirpation in the jurisdiction due to a fairly restricted range, relatively few populations or occurrences, recent and widespread declines, threats, or other factors.
abundant species) and considered the BNR a mussel diversity refuge in the Ozark Highlands.The definitions for status rank categories follow Master et al. [33]: S1, Critically Imperiled-at very high risk of extirpation in the jurisdiction due to very restricted range, very few populations or occurrences, very steep declines, severe threats, or other factors; S2, Imperiled-at high risk of extirpation in the jurisdiction due to restricted range, few populations or occurrences, steep declines, severe threats, or other factors; S3, Vulnerable-at moderate risk of extirpation in the jurisdiction due to a fairly restricted range, relatively few populations or occurrences, recent and widespread declines, threats, or other factors.1.The 2006 and 2020-21 BNR freshwater mussel sampling site codes, river kilometer (Km), drainage area, slope, percent development (Dev), impervious cover (Imp), pasture (Pas), mean annual precipitation (MAP), and sampling efforts.The symbol (X) denotes sites that were sampled or pairwise comparisons were conducted.This study aimed to evaluate the changes in freshwater mussel population and community characteristics between 2006 and 2020-2021 at defined mussel assemblages (or mussel beds) in the BNR using a stratified random sampling protocol.Our first objective was to evaluate the overall characteristics of mussel bed persistence, sampling confidence levels, and study-wide relative abundances between sampling events and to compare size and size-frequency distributions of species with study-wide relative abundances ≥1%.We expected that (1) all sites would contain mussel beds, that (2) sites would be adequately sampled at ≥80% confidence levels, and that (3) study-wide species' relative abundances, sizes, and size-frequency distributions would be similar between events.Our second objective was to compare 10 overall site assemblage variables between events using paired t-tests.We expected no significant differences in any variables between events.Our third objective was to compare site-level mean density, richness, and diversity indices between events using pair-wise Mann-Whitney U statistics.We expected no significant differences between events.Our fourth objective was to compare assemblage composition between events using Non-Metric Multidimensional Scaling.We expected sampling sites between events would cluster together, and sites would separate across an upstream-to-downstream gradient with sites in close geographic proximity clustering together.

Overview
Quantitative site sampling in 2006 and 2020-21 was proceeded by large-scale qualitative sampling in 2004-2005 [31] and 2019-2020 [34], respectively.After 12 monitoring sites were identified in 2006 (Table 1; Figure 1), quantitative sampling was conducted via a six-step process of qualitative reconnaissance, mussel bed delineation, mussel bed mesohabitat determination and stratification, stratified random sampling, sample processing, and in-field sampling confidence level determination with additional sampling if needed.The 2020-21 quantitative sampling followed the same methods; however, two 2006 monitoring sites no longer persisted, and three additional sites were added (Table 1).

Study Area
Stream geomorphological characteristics across the upstream to downstream gradient of sampling sites ranged from stream orders 4 to 6, drainage areas 365-3357 km 2 , and slopes 5.097-1.257m/km, respectively [35].Meanwhile, land use land cover characteristics (i.e., % development, impervious cover, and pasture land) generally increased from the upstream to downstream gradient of monitoring sites and ranged from 6.68 to 14.80% pasture, 2.38-3.38%development, and 0.14-0.27%impervious cover [35].The entire BNR watershed mean annual precipitation is 121 cm/year [35].

Monitoring and Sampling Techniques 2.2.1. Qualitative Reconnaissance
Qualitative reconnaissance and sampling consisted of time-constrained walking visuals and snorkeling and/or diving searches that have been shown to maximize species richness determinations while remaining cost-effective [36].Mussel assemblage sample sites were searched for no less than one person-hour, and all native mussel individuals were enumerated and identified using nomenclature of Williams et al. [37], as modified by the FMCS Common and Scientific Names subcommittee [38].

Mussel Bed Delineation
Field personnel located areas inhabited by mussels during qualitative reconnaissance, estimated densities and species richness, and reviewed longitudinal distribution of assemblages within the BNR prior to site selection for the mesohabitat determinations.The area inhabited (length X width) by assemblages with estimated mean densities >1 mussel/m 2 (i.e., mussel beds) was determined, GPS coordinates were recorded at upstream, midsection, and downstream locations, and the bed perimeter was marked with weighted floats or survey flags/tape.Sites selected were distributed along the length of the study area.

Mesohabitat Determinations and Stratification
Mussel bed stratification was based on habitat characterization and the resulting mesohabitat determination was performed at each site following the Basin Area Stream Survey (BASS) protocol [39,40].This methodology included delineating and recording each mesohabitat type (e.g., riffle, run, pool, backwater), qualitatively sampling microhabitat features such as length, width, depth, substrate, embeddedness, and instream cover, and riparian cover characteristics of bank stability and canopy coverage.We also conducted a 100-sample Wolman pebble count [41] to characterize the reach with the percentage of fine sediment, sand (fine sediment ≤ 2 mm), gravel (3-64 mm), cobble (65-256 mm), boulder (>256 mm), and bedrock.However, for the purposes of this study, we only used the microhabitat type information and associated stream length, width, and depth information for the next step of stratified random sampling.

Stratified Random Sampling
Sampling used simple random stratified methods [36,42].Once assemblage areas were defined, an x, y coordinate grid was established, and quantitative sample locations were selected utilizing random number tables.If a mussel assemblage occupied diverse habitat types (i.e., riffle/pool, different substrates) or had substantially differing morphology (e.g., widths), the sampling area was stratified with each stratum sampled proportionally to area.A minimum of 10 1 m 2 quadrat samples were taken at each site (unless precluded by site disturbance concerns, e.g., BR0430 in 2006), with additional quadrats added to reach an 80% confidence interval [43] for mussel density.Additionally, a 25-quadrat maximum was implemented to address substrate and mussel assemblage disturbance concerns.Sampling was performed by placing the 1 m 2 quadrat at the randomly selected location with a 0.25 m 2 quadrat randomly placed within.Water depth and velocity were measured in the quadrat center.Boulders and cobble were removed from the sampling area, the remaining substrate of the 0.25 m 2 was extracted to a 10 cm depth and sieved to detect juvenile mussels, and the remaining portion of the 1 m 2 quadrat was searched by hand to a depth of approximately 10 cm.

Sample Processing
All mussels encountered within the 1 m 2 quadrat were bagged, brought to the surface, and identified to species [37], and updated in common data sets using FMCS [38].The length for each mussel was measured to the nearest 0.1 mm with digital calipers.Native mussels were then replaced into the substrate from which they were taken.Voucher specimens were curated in the Arkansas State University Museum of Zoology and cataloged with the NPS.Pictures were taken of each site to document upstream and downstream reaches.Site codes were assigned using BNR standardized codes, and river kilometers (rkm) were determined from the USGS National Hydrography Dataset (NHDPlus HR) with rkm 0 at the mouth.

In-Field Sampling Confidence Level Determination and Additional Sampling
Sampling in 2020-21 employed predefined assemblage areas consistent with the 2006 delineations.The sampling confidence level reflected the extent to which the data accurately represented the true biodiversity characteristics of the mussel assemblage.After preliminary sampling of a minimum of 10 1-m 2 quadrats, a final sample size number was determined utilizing the equation provided by Southwood [43].The desired relative error selected for this survey was 20%.If additional samples were needed beyond the initial 10 1 m 2 samples based on the estimates [43], samples were added up to a maximum of 25 total samples.

Data Analysis
Descriptive summary statistics of sampling sites' (beds') variables were conducted using data analysis tools in Microsoft Excel.Species richness, evenness (E), Shannon's diversity (H'), and Simpson's diversity (D) were calculated in PCOrd [44].Population and community numerical standing crop estimates and 95% confidence intervals were calculated following Huebner et al. [45].We evaluated our sampling efforts at each site and event (event X site) by estimating the relative error and the number of samples to achieve 80% and 90% confidence levels for density and species richness using methods described in Southwood [43].To compare 2006 and 2021-21 median size and size cumulative frequency distributions of each species with study-wide relative abundance ≥1%, we used Mann-Whitney U-signed rank SigmaStats Version 4 [46] and Kolmogorov-Smirnov D tests in ExcelStats [47], respectively.To test for significant differences in bed variables between the two (2006 versus 2020-2021) sampling events, we used a paired two-tailed t-test while also testing for Shapiro-Wilk normality and estimated sample size power using Sigma-Stat, Version 4 [46].Pair-wise comparisons of mean density and mean richness between sampling events of each site were assessed using SigmaStat, Version 4 [46] to conduct signed rank tests (Mann-Whitney U) due to most datasets having non-normal data based on Shapiro-Wilk tests.Finally, to evaluate if event X site communities changed between 2006 and 2020-21, we conducted a Non-Metric Multidimensional Scaling (NMS) analysis on the community-level data using PCOrd [44].Analysis options for the NMS were set as Morisita-Horn distance measure, 500 iterations, random starting coordinates, autopilot, and a slow and thorough speed setting.We used the Morisita-Horn distance measure due to a high number of zeros in our data set and unequal sampling efforts at each site.

Mussel Assemblage Sampling Site Variables' Descriptions
The defined mussel bed area (m 2 ) ranged from a low of 48 m 2 to a high of 622 m 2 in 2006 and a low of 48 and high of 3625 m 2 in 2020-21 (Table 3).In 2006, we used 6 to 25 1 m 2 quadrats to sample mussel beds, representing 4 to 10% of the defined mussel bed area.Meanwhile, in 2020-21, we used 10 to 25 1 m 2 quadrats to sample the mussel beds, representing 1 to 21% of the defined mussel bed.The mean density ranged from 1.In 2006, the number of species with only one occurrence at a site (singlet) ranged from zero to five, while the number of species with only two occurrences at a site (doublet) ranged from zero to three.Meanwhile, the number of species with only one occurrence at a site (singlet) ranged from one to five while the number of species with only two occurrences at a site (doublet) ranged from zero to two in 2020-21.

Sampling Confidence Levels of 2006 and 2020-21
In 2006, the sample size confidence level based on density ranged from 57.2 to 91.1 with 9 of 12 sampling efforts meeting the 80% confidence level target, and the confidence level based on richness ranged from 72.2 to 96.4 with 9 of 12 sampling efforts meeting 80% (Table 4).Meanwhile, the 2020-21 sample size confidence level based on density ranged from 67.5 to 91.6 with 9 of 13 sampling efforts meeting the 80% confidence level, and the confidence level based on richness ranged from 75.6 to 94.3 with 11 of 13 sampling efforts meeting 80% (Table 4).For both sampling density and richness, we provided the number of samples needed to reach 80 and 90% confidence levels to compare the required efforts (Table 4).

Size and Cumulative Frequency Distributions
Mann-Whitney U pairwise comparisons between 2006 and 2020-21 median lengths for 11 species resulted in seven statistically significant comparisons with six showing significantly a higher length in 2020-21 samples (Table 5).Pair-wise Kolmogorov-Smirnov length cumulative frequency distribution analysis resulted in nine significantly different pair-wise comparisons (Table 5).Six species (A.ligamentina, C. hesperus, L. reeveiana, L. costata, P. occidentalis, and V. pleasii) showed significantly different distributions with low to no recruitment (Table 5; Figure 2).Two species (C.tuberculata and E. dilatata) showed significantly different distributions with different recruitment patterns between 2006 and 2020-21 (Table 5; Figure 3).Toxolasma lividum showed significantly different distributions with higher recruitment but fewer larger individuals in 2020-21 (Table 5; Figure 4).Meanwhile, two species (F.ozarkensis and P. sintoxia) showed similar distributions between the two sampling events (Table 5; Figure 5).5) with more recruitment but fewer larger individuals in 2020-21.Bins represent a dedicated 5 mm size range in which the size frequencies were summed.

The Paired t-Test of Common 2006 and 2020-21 Sites
We conducted 10 paired t-test comparisons for 10 overall variables (Table 3) associated with the common 2006 and 2020-21 sampling sites, and all 10 datasets had normal distributions based on Shapiro-Wilk Normality (p ≤ 0.05); however, almost all datasets had low statistical power (Table 6).While none of the paired t-test comparisons were statistically significant, four variable means were higher (sample size, % bed sampled, S1-S3 richness, and S1-S3 abundance) and six variable means were lower (density, total richness, evenness, Shannon's diversity, Simpson's diversity, and number of singlets and doublets) in the 2020-21 monitoring efforts.5).Bins represent a dedicated 5 mm size range in which the size frequencies were summed.

The Paired t-Test of Common 2006 and 2020-21 Sites
We conducted 10 paired t-test comparisons for 10 overall variables (Table 3) associated with the common 2006 and 2020-21 sampling sites, and all 10 datasets had normal distributions based on Shapiro-Wilk Normality (p ≤ 0.05); however, almost all datasets had low statistical power (Table 6).While none of the paired t-test comparisons were statistically significant, four variable means were higher (sample size, % bed sampled, S1-S3 richness, and S1-S3 abundance) and six variable means were lower (density, total richness, evenness, Shannon's diversity, Simpson's diversity, and number of singlets and doublets) in the 2020-21 monitoring efforts.5).Bins represent a dedicated 5 mm size range in which the size frequencies were summed.

The Pairwise Density and Richness Comparisons of 2006 and 2020-21
For the 10 commonly sampled sites, we conducted pairwise comparisons for sampling density and richness using the non-parametric Mann-Whitney U statistic because most of our datasets were found to have non-normal distributions based on Shapiro-Wilk Normality (p ≤ 0.05).Of the 20 pair-wise comparisons, only 6 comparisons were statistically significant (Table 7).The 2020-21 sampling was found to be significantly lower for BR0260 richness, BR0626 density and richness, and BR0820 density and richness but was statistically higher for BR0579 density (Table 7).Unlike the paired t-test comparisons, no general nonsignificant trends in higher or lower density and richness were observed between the two sampling events.

Overview
We identified four major findings while comparing our 2006 and 2020-21 BNR monitoring sites.First, our sampling protocol to achieve a targeted 80% confidence level appeared to be relatively robust, but additional samples were required in some cases.Second, BNR mussel assemblage compositions and structures were relatively stable through the ~15-year sampling window; however, there are concerns for mussel bed persistence due to ~17% loss of sampled beds.Third, in addition to two 2006 sampling sites no longer persisting in 2020-21, seven of the remaining twenty-three sites were outliers based on their freshwater mussel composition and characteristics.Finally, based on both relative abundance and the NMS, assemblage composition had changed with three abundant species declining (A.ligamentina, L. costata, and P. occidentalis), and four species increasing (C.hesperus, C. tuberculata, E. dilatata, and V. pleasii) between sampling events.

Sampling Effort
To provide guidance for gathering statuses and monitoring information on freshwater mussels, Smith and Strayer [36] proposed adopting standardized sampling designs and methods based on the needs of the study.For the BNR, we used a stratified random sampling design which has been commonly used for surveying mussels in Arkansas rivers since the early 1990s [31,[48][49][50][51][52][53].Overall, our sampling efforts during both monitoring events were robust, with 18 of 25 event X site sampling efforts achieving 80% or higher confidence levels for site mean density and 20 of 25 event X site sampling efforts achieving 80% or higher confidence levels for site mean richness.The 80% confidence level for sample density and richness provided valuable information on the status of freshwater mussels, sufficient to inform management decisions at a reasonable expense.For example, achieving 90% confidence levels for both density and richness would require doubling the sample efforts in most cases.These sampling efforts and confidence level results are consistent with other studies that used similar sampling designs across small to large streams [42,52].

Species and Assemblage Composition and Structure Changes
The paired t-tests of 10 variables comparing 2006 and 2020-21 sampling efforts for 10 common sites showed no significant differences.Only 6 of 20 Mann-Whitney U tests comparing the mean density and richness of 2006 versus 2020-21 sites were statistically different.However, in addition to 5 of 6 Mann-Whitney U tests indicating that 2020-21 richness or density was significantly lower than 2006 samples, most assemblage structure variables in the paired t-test were lower in 2020-21, though not statistically significant.The exceptions were that sample size, % bed sampled, S1-S3 richness, and S1-S3 abundance were slightly higher in 2020-21 than in 2006.Furthermore, among the dominant taxa (>1% relative abundance), five taxa declined, four taxa increased, and one species stayed the same in relative abundance between 2006 and 2020-21.Declining relative abundance taxa included A. ligamentina, F. ozarkensis, L. costata, P. sintoxia, and P. occidentalis, with A. ligamentina decreasing from 9.3 to 0.3% and P. occidentalis decreasing from 31.2 to 22.5%.Meanwhile, A. ligamentina, L. costata, and P. occidentalis also showed statistically different size-frequency distributions with lower recruitment.Conversely, the increasing relative abundance of dominant taxa included C. hesperus, C. tuberculata, E. dilatata, and V. pleasii.In contrast, C. hesperus and V. pleasii size-frequency distributions indicated different size frequency distributions between sampling events with lower recruitment in the 2020-21 samples while C. tuberculata and E. dilatata showed different distributions with different recruitment frequencies between sampling events.Cambarunio hesperus increased from 8.8 to 18.3%, C. tuberculata increased from 5.8 to 9.2%, E. dilatata increased from 1.5 to 4.5%, and V. pleasii increased from 3.9 to 13.5%.Furthermore, L. reeveiana maintained its relative abundance at 15.1% between the ~15-year sampling events.In terms of species observations, five species were observed in 2006 (C.aberti, E. triquetra, F. flava, T. verrucosa, and T. parvum) but not in 2020-21, while two species were observed in 2020-21 (A.viridis and L. recta), but not in 2006.In all cases, these observations were low (<11) with A. viridis, C. aberti, E. triquetra, L. recta, and T. verrucosa all represented by one-three individuals.The 2006 E. triquetra record was from one individual and could have been a misidentified A. viridis [54].
Stream size and life history strategies influence the occurrence and distribution of North American freshwater mussels [55].Haag [55] identified dominant (10% or more of the relative abundance) mussel species across small, medium-sized, and large Mississippian region streams.Three life history strategies have been proposed for North American freshwater mussels that include opportunistic (disturbed and unstable but productive habitats), periodic (unproductive or high variability/stress habitats), and equilibrium (stable, productive habitats) strategies.Species that declined between 2006 and 2020-21 were mostly equilibrium-strategy species (A.ligamentina, F. ozarkensis, P. sintoxia, and P. occidentalis) that are known to inhabit mid-sized streams, while L. costata is considered a periodic strategist that mostly inhabits small streams and to a lesser extent mid-sized streams.Two of the four species that increased between 2006 and 2021-22 do not have a clear characterized life history strategy (C.tuberculata and E. dilatata) but C. hesperus is an opportunistic strategist that inhabits mostly small streams and V. pleasii can be classified as a periodic-strategy species and is not known as a dominant species in any stream size category.For the rare species observed (C.aberti, E. triquetra, F. flava, T. verrucosa, and T. parvum), none were expected to be dominant (>10% relative abundance) small, medium, or large river species except for F. flava, but they represented each life history strategy: opportunistic (T.parvum), periodic (E.triquetra), equilibrium (C.aberti, F. flava, T. verrucosa).
This suggests that while there have been small changes, mostly declines, in assemblage characteristics such as diversity metrics, there have been modest to large changes in recruitment and assemblage composition in BNR mussels.While some of these compositional changes may be attributed to the small number differences in sampling sites between the ~15-year sampling events, the nearly 10% relative abundance declines of A. ligamentina (G5 and S4) and P. occidentalis (G3/G4 and S3) are noteworthy, as are the nearly 10% relative abundance increases of C. hesperus (G5 and S3) and V. pleasii (G3/G4 and S3).For those taxa that declined between the ~15-year sampling events, it is possible that stream conditions have become less stable than in the past.For example, both A. ligamentina and P. occidentalis have been reported to be thermally and drought sensitive [56,57].Given that rising temperature is a growing concern in the BNR [26,28], increases in the magnitude, duration, and frequency of these events pose major risk to long-lived species.Meanwhile, for the species that increased in relative abundance, they are classified as non-equilibrium-strategy species that are better suited to higher disturbance and stress.These overall relative abundance changes, plus the observation that S1-S3 richness and relative abundance increased from 2006 to 2020-21, suggest that freshwater mussel assemblage changes are occurring in the BNR.These changes warrant further monitoring of mussel assemblages and environmental variables (e.g., geomorphological channel shifts and erosion, nutrients, climate and land use land cover change).Continued habitat monitoring within mussel bed areas in conjunction with quantitative sampling events may provide more insight on these shifts in composition.

Mussel Bed Persistence and Outliers
Between the two sampling events, we sampled 25 event X site combinations with 10 common sites across the ~15-year interval.In 2020-21, we added three additional longterm quantitative monitoring sites with BR0256 added to the upstream extent of monitoring and BR0578 and BR0581 added due to their bed size, density, and richness.Sites BR0350 and BR0810 were established in 2006, but neither site was located during the 2020-21 sampling effort, and these two sites are no longer believed to be extant.The loss of these two mussel bed sites across the 15-year sampling interval represents a 17% loss of sampled mussel beds.Anecdotal observations at BR0350 and BR0810 suggest substantial geomorphological changes occurred in the ~15-year time interval.The BNR is a relatively high gradient stream system (with steeper gradients and rapid elevation changes in the upper reaches in the Boston Mountains) and channel-changing flood events are known to occur [31].Discharge measurements at USGS monitoring stations revealed two moderate flood stages (<10.7 m) and three minor flood events (<8.2 m) occurred between sampling events and satellite imagery at each site documented channel changes.Across the two mussel bed sites along the BNR, varying flood-event intensities may have collectively degraded the habitat, rendering it unsuitable due to physical disruption, sediment accumulation, bank erosion, water quality fluctuations, and food availability.Dedicated monitoring such as the current study paired with side-scan sonar surveys can document mussel bed persistence and habitat stability in the BNR.
For our NMS analysis, we expected that event X site samples would show an upstreamto-downstream distribution gradient with spatial and temporal assemblage similarities.However, aside from the downstream sites (BR800s) being distributed at higher values of Axis 2 and lower values of Axis 3, there was little evidence of an upstream-to-downstream distribution for our NMS analysis.In retrospect, because our sampling sites were distributed within the fourth to sixth stream order section of the Buffalo River, our pattern is consistent with patterns observed by Haag [55] in which certain genera and species are known to differentially inhabit small, medium, or large streams where our BNR sampling occurred withing Haag's "medium sized streams" assemblage structure.Furthermore, our NMS analysis resulted in four event x site samples (i.e., 06BR0810, 2021BR0256, 2021BR0304 and 2021BR0408) being outliers.The 06BR0810 was characterized by high numbers or the presence of A. ligamentina, A. marginata, A. plicata, C. aberti, F. flava, F. ozarkensis, L. costata, L. cardium, L. reeveiana, P. sintoxia, and P. occidentalis, and was the only site to include T. verrucosa.It should also be noted that the 06BR0810 mussel assemblage (i.e., bed) was one of two sites that was not found in the 2020-21 sampling event and is most likely extirpated given no mussels were found.The 2021BR0256 assemblage was influenced by only having five species present (F.ozarkensis, L. cardium, L. reeveiana, L. costata, and P. occidentalis), which were low in abundance, and none of which were unique to 2021BR0256.The 2021BR034 assemblage was influenced by being represented by only six species with C. hesperus dominating (n = 50) the abundance while the other four species had abundances of one-three individuals.The 2021BR0408 assemblage was influenced by relatively low abundances of 10 commonly sampled BNR species.

Conclusions
While our analysis suggests that BNR mussels were relatively stable between the 15-year sampling interval, mussel bed persistence and trends in population and assemblage characteristics warrant continued monitoring.Investigation into areas such as hydrodynamic modeling, changes in fish host communities, and dissolved oxygen and pollutant load dynamics should be considered.With only a relatively small portion of the BNR watershed managed by the NPS, a catchment approach related to tributaries feeding into the mainstem may be necessary to include in a monitoring plan that will involve multiple entities such as local stakeholders and state agencies.Long-term quantitative studies of mussel assemblages paired with monitoring environmental drivers in the BNR are essential in documenting changes and providing timely data points for conservation and management decisions.
3 to 22.6 individuals/m 2 in 2006 and from 2.4 to 23.5 individuals/m 2 in 2020-21.The community numerical standing crop ranged from 213 (±124) to 14,801 (±2227) individuals in 2006 and from 298 (±127) to 83,558 (±15,369) individuals in 2020-21.The total richness ranged from 4 to 16 in 2006 and from 5 to 12 in 2020-21.The number of S1, S2, and S3 richness ranged from three to six in 2006 and from three to seven in 2020-21.Furthermore, the combined S1, S2, and S3 relative abundance ranged from 39.5 to 93.4% in 2006 and from 33.8 to 95.7% in 2020-21.The mussel bed assemblage evenness ranged from 0.560 to 0.894 in 2006 and from 0.363 to 0.875 in 2020-21.Shannon's diversity (H') ranged from 0.892 to 2.260 in 2006 and from 0.651 to 2.015 in 2020-21.Simpson's diversity (D) ranged from a low of 0.4930 to a high of 0.8741 in 2006 and from a low of 0.2758 to a high of 0.8423 in 2020-21.

Figure 2 .
Figure 2. Six species (A.ligamentina, C. hesperus, L. reeveiana, L. costata, P. occidentalis, and V. pleasii) that showed significantly different length-frequency and cumulative frequency (CFD) distributions (Table5) with low to no recruitment in 2020-21.Bins represent a dedicated 5 mm size range in which the size frequencies were summed.

Figure 2 .
Figure 2. Six species (A.ligamentina, C. hesperus, L. reeveiana, L. costata, P. occidentalis, and V. pleasii) that showed significantly different length-frequency and cumulative frequency (CFD) distributions (Table5) with low to no recruitment in 2020-21.Bins represent a dedicated 5 mm size range in which the size frequencies were summed.

Figure 3 .
Figure 3. Two species (C.tuberculata and E. dilatata) that showed significantly different lengthfrequency and cumulative frequency (CFD) distributions (Table 5) with different recruitment patterns between 2006 and 2020-21.Bins represent a dedicated 5 mm size range in which the size frequencies were summed.

Figure 4 .
Figure 4. Toxolasma lividum showed significantly different length-frequency and cumulative frequency (CFD) distributions (Table5) with more recruitment but fewer larger individuals in 2020-21.Bins represent a dedicated 5 mm size range in which the size frequencies were summed.

Figure 3 . 12 Figure 3 .
Figure 3. Two species (C.tuberculata and E. dilatata) that showed significantly different lengthfrequency and cumulative frequency (CFD) distributions (Table 5) with different recruitment patterns between 2006 and 2020-21.Bins represent a dedicated 5 mm size range in which the size frequencies were summed.

Figure 4 .
Figure 4. Toxolasma lividum showed significantly different length-frequency and cumulative frequency (CFD) distributions (Table5) with more recruitment but fewer larger individuals in 2020-21.Bins represent a dedicated 5 mm size range in which the size frequencies were summed.

Figure 4 .
Figure 4. Toxolasma lividum showed significantly different length-frequency and cumulative frequency (CFD) distributions (Table5) with more recruitment but fewer larger individuals in 2020-21.Bins represent a dedicated 5 mm size range in which the size frequencies were summed.

Figure 5 .
Figure 5. Two species (F.ozarkensis and P. sintoxia) that showed similar length-frequency and cumulative frequency (CFD) distributions between 2006 and 2020-21 (Table5).Bins represent a dedicated 5 mm size range in which the size frequencies were summed.

Figure 5 .
Figure 5. Two species (F.ozarkensis and P. sintoxia) that showed similar length-frequency and cumulative frequency (CFD) distributions between 2006 and 2020-21 (Table5).Bins represent a dedicated 5 mm size range in which the size frequencies were summed.

Table 3 .
Mussel bed and sampling summary statistics for 2006 and 2020-21.NS is not sampled.

Table 5 .
Mann-Whitney (M-W) U and Kolmogorov-Smirnov (K-S) D statistics for study-wide length and length cumulative frequency analysis for 11 species with ≥1% pooled relative abundance sampled at 10 common monitoring sites in 2006 and 2020-21.N = number of individuals.

Table 6 .
Ten paired t-test comparisons for 10 overall variables (Table3) associated with the 10 common 2006 and 2020-21 sampling sites.

Table 6 .
Ten paired t-test comparisons for 10 overall variables (Table3) associated with the 10 common 2006 and 2020-21 sampling sites.

Table 7 .
Pairwise comparisons for sampling density and richness using the non-parametric Mann-Whitney U statistic; most datasets had non-normal distributions based on Shapiro-Wilk Normality (p ≤ 0.05).N = m 2 quadrat samples.Bold text indicates a significant difference (p ≤ 0.05).

Table 8 .
Species correlations (r, r 2 , and tau) with axes 1, 2, and 3. Bold indicates significant correlations with each of the axes.