Biodiversity Assessment of a Mississippi River Backwater Complex Using eDNA Metabarcoding

1. Introduction

Biodiversity assessment of complex integrated freshwater ecosystems provides essential data for assessing and quantifying spatiotemporal patterns and long term changes in freshwater ecosystems (Jarvis et al., 2023). However, freshwater ecosystems present a range of challenges to effective sampling and monitoring (Radinger et al., 2019; Schramm et al., 2016). The Mississippi River, as with many big river systems, has been severely impacted by anthropogenic changes that threaten its biological diversity and ecosystem function (Best, 2019; DuBowy, 2013). In the present study, we employed environmental DNA (eDNA) metabarcoding (Deiner et al., 2017; Ruppert et al., 2019; Schenekar, 2023) to assess the biodiversity in one network of Mississippi lowland habitats.

The northern extent of the Gulf Coastal Plain of North America begins at the confluence of the Mississippi and Ohio Rivers (Schramm & Ickes, 2016). The Mississippi River below this confluence is the center of an expansive alluvial plain ecosystem that historically was dominated by bald cypress and tupelo bottomland forests and wetlands, and was connected to the river by seasonal inundation of floodwaters (Pflieger, 1997). Beginning in the early 1900s efforts were underway to clear and drain these wetlands for agricultural purposes (Pierce et al., 2012). Contemporary major levee systems effectively separate the alluvial plain ecosystem from its major rivers (e.g. the Commerce MO – St. Francis River Levee System, constructed in 1950) (DuBowy, 2013; USACE, 2014), and 96% of bottomland habitats in the region have been drained. Wetland habitats have become fragmented, and few remain connected to river floodplains (Olson et al., 2016b, 2016a).

The Gulf Coastal Plain is an ichthyofaunal diversity hotspot, with an assortment of lowland endemics and large river species (Cross, F. B., Mayden, R. L., & Stewart, 1986; Isphording & Fitzpatrick Jr., 1992; Jenkins et al., 2015; Noss et al., 2015). The endemic fishes of the region include species that may have once been common but have been negatively impacted by the alteration and destruction of their habitats (Sowa et al., 2007; Warren et al., 2000). Annual flooding in rivers make bottomlands biologically productive as water, nutrients, and organisms are exchanged between the main river channels and backwater habitats (Junk et al., 1989). Additionally, vast river floodplains provide critical habitat to a diverse array of large-bodied river fishes which exhibit behaviors and life-history attributes that take advantage of seasonally available floodplain habitats for spawning and juvenile growth (Allen et al., 2020; Phelps et al., 2015; Schramm & Eggleton, 2006). Efforts to convert floodplains to agricultural land, and channelize rivers for commerce have largely eliminated the seasonal floods that many large river species depend on for reproduction and early life stages (Humphries et al., 1999; Sparks, 1999).

Floodplain habitats often support a greater diversity of species than can be reliably observed in main river channels, however, standard field sampling methods may be inadequate to effectively capture such diversity (Dettmers et al., 2001; Phelps et al., 2015). Specifically, fish assemblage surveys are limited by the challenges of sampling these heavily vegetated, complex, fragmented, and often inaccessible lowland habitats (Schramm et al., 2016). Despite traditional sampling difficulties, efforts to assess the composition of the species present in these remnant habitats are critical for effective management.

Although lowland habitats present substantial challenges for traditional sampling efforts, use of eDNA-based survey methods may reduce such challenges. eDNA is a genetically based, non-invasive biomonitoring tool that can be used to assess the species assemblage in an area of interest (Mächler et al., 2019). The method relies on eDNA that has been released into the environment via shed skin, saliva, blood, excrement, or gametes of living or decomposing organisms (Rees et al., 2014). The detection zone - the size of the area surrounding an eDNA source where it can be detected, is dependent on water flow, mixing characteristics and rate of eDNA degradation. Due to lentic conditions that often characterize river bottom wetlands, the detection zone may be on the order of 100 m or less (Harrison et al., 2019). Studies have shown that eDNA bound to fine silt and clay particles persists longer than suspended eDNA, and can accumulate and persist in an environment (Nevers et al., 2020). Accumulation and persistence of eDNA in sediments may be a consideration for interpretation of data, particularly for high biomass species that generate high concentrations of eDNA (Turner et al., 2015).

The deployment of eDNA-based survey methods has become common practice in aquatic systems, and can be used in addition to or as an alternative to traditional sampling techniques (Cilleros et al., 2019). eDNA can be filtered and extracted from a water sample, amplified, sequenced, and aligned to a reference database to infer the species that are present. While the efficacy of eDNA surveys has been well documented in riverine and marine systems (García-Machado et al., 2022; Lecaudey et al., 2019), few studies have implemented this technique in freshwater wetlands (e.g. Kačergytė et al., 2021).

In this study we employed eDNA metabarcoding methods to sample ichthyofauna in a Mississippi River bottom wetland system and assess the distribution of species among multiple habitats with varying vegetation levels, seasonal flow dynamics, and connections to the Mississippi River. We aimed to assess the alpha and beta-diversity of fish assemblages in a relatively intact remnant natural Mississippi River floodplain ecosystem. Our objectives were to quantify and compare species richness among the different habitats, and test for fish assemblage structure across habitats and seasons (i.e. spring and fall). We compared the species detected using eDNA metabarcoding methods to a database of historical records of fish species collected in the studied area from 1940 to the present. These comparisons provide an assessment of the efficacy of eDNA metabarcoding in this system and a standardized evaluation of the ichthyofaunal assemblages across seasons and habitats, which is rarely accomplished in these imperiled habitats due to limits of traditional sampling.

2. Materials and Methods

2.1. Study System

The Black Island Conservation Area is a 2,602 hectare wetland reserve located within the Black Island river bend in southeastern Missouri, and it is managed by the Missouri Department of Conservation. The complex of bottomland habitats within Black Island are bounded to the east by the Mississippi River channel, and to the west by the Commerce MO – St. Francis River Levee System (Figure 1). The Black Island wetland habitats include an oxbow lake bayou complex, an expanse of low relief seasonally flooded lowlands, and a slough complex that is directly connected to the Mississippi River. The change in elevational relief across the entire Black Island complex ranges as high as ~10 m above the base flow of the Mississippi River, with the majority of seasonally flooded lowland areas elevated about 3 to 4 m above the base flow (USGS, 2024).

We sampled the Black Island complex on two dates, May 30th and October 7th of 2022. During the spring sampling event the Mississippi River stage was 7.1 m at the nearby Caruthersville monitoring station (Figure 2). During the fall sampling event the Mississippi River stage was 0.08 m. This was nine days before the monitoring station recorded the second lowest all-time record of -0.55 m.

The Wolf Bayou complex is comprised of three deeply incised oxbow lake bayous (Samples, Hosner, and Wolf) connected by narrow channels (Table 1, Figure 1b). Except where there are stream inlets, most of the shoreline of all three bayous is relatively steep, resulting in minimal seasonal variation in shoreline characteristics and overall total surface area.

Robinson Lake (Table 1, Figure 1c) is a shallow impoundment formed by a low retention dam. In the spring the shoreline was inundated into the heavily vegetated border, composed of shrubs and small trees, to a penetration of tens to hundreds of meters around the entire perimeter. During the fall sampling period the lake surface dimension was reduced to less than 300 m long and 50 m wide. There was no vegetation in the fall as the shoreline had receded substantially relative to the spring.

We sampled a slough along the margin of Black Island that was connected to the Mississippi River (Table 1, Figure 1d). The slough was cut off from the river in the fall and reduced to isolated pools. We also sampled a steep-sided ditch that was originally constructed to facilitate drainage of agricultural lands before the wetland reserve was established. The ditch was connected to the slough via a large culvert under an access road. The water depth in the ditch matched the depth of the slough, and water was flowing from the ditch into the slough when we sampled in the spring.

2.2. eDNA Field Sampling

We utilized eDNA metabarcoding to assess fish diversity in the Black Island complex. Water samples from the bayous, lake, and ditch were collected by launching a small jon boat equipped with an electric trolling motor. The slough was sampled with the use of waders. At each sample site, three 500 ml water samples were collected in separate sterile bottles. A 2 m extension pole was used to collect the water samples to avoid water contamination from the boat or waders, and water was collected in an arcing motion by submerging each bottle approximately 10 cm below the water surface. We did not sample the substrate. Water samples were immediately pressure filtered through enclosed 0.45 µm polyvinylidene fluoride (Millipore Sigma) filters using 50 ml luer-lock syringes until the volume reached 500 ml or the filter clogged. Field negatives consisted of 500 ml of deionized (DI) water which were filtered in the field to monitor potential field contamination. Following filtration, the filter housing was flooded with 95% molecular grade ethanol, sealed with parafilm, and stored on ice for transport to the lab (Williams et al., 2016). Disposable gloves were worn during the water collection and filtration and changed between sites. Water quality parameters including standard conductivity (µS/cm), dissolved oxygen (mg/L), and temperature (°C) were logged at each sample location using a YSI probe. A Secchi disk was used to measure water clarity (i.e., Secchi depth) and water column depth at each sample site was recorded.

2.3. DNA Extraction, PCR Amplification, Sequencing

Sealed filters were stored at 4^oC until they were extracted. Extraction was conducted within 7 days of sample collection using a modified protocol of the Qiagen Blood and Tissue Kit (Qiagen) for enclosed filters outlined in Spens et al. (2017). To minimize cross contamination, extractions were carried out in a dedicated clean lab with separate dedicated equipment from all post-PCR steps, and all work surfaces and equipment were decontaminated using bleach between extractions (Goldberg et al., 2016). In wetland habitats high concentrations of suspended organic and inorganic matter can interfere with water filtration, speed up the degradation of suspended eDNA in the water column, and inhibit subsequent PCR amplification (Kumar et al., 2022). To address the latter challenge, PCR inhibitors were removed using a OneStep PCR Inhibitor removal Kit (Zymo) after DNA extraction.

Extracted eDNA was PCR amplified in separate sets of reactions using two universal fish mitogenome primer sets: Mifish 12s primers (Miya et al., 2015) and the 16s rRNA primers Chord 16s F TagA and Chord 16s R short primers (Deagle et al., 2009). Both primer sets were 5’ tagged with sequences to provide binding sites for the Illumina sequencing primers. eDNA samples were amplified using AmpliTaq Gold 360 DNA polymerase (Thermo Fisher Scientific) in 50 µL total volume reactions. Reactions were separated into six separate tubes for thermal cycling, and then PCR products were recombined (Ruppert et al., 2019). Potential sources of lab contamination were monitored by inclusion of lab-negative water controls. PCR reactions were set up in a UV-sterilized AirClean 600 PCR Workstation (ISC BioExpress). Cycling conditions for Mifish were denaturation at 95° C for 5 min, then 33 cycles of 95° C for 20 sec, 65° C for 20 sec, 72° C for 1 min, then a final extension step of 72° C for 5 min, and a final hold of 4° C. Cycling conditions for Chord 16s were denaturation at 95° C for 5 min, then 45 cycles of 95° C for 25 sec, 55° C for 30 sec, 72° C for 1 min, then a final extension step of 72° C for 6 min, and a final hold of 4° C.

Aliquots of all PCR products and negative controls were electrophoresed on 2% agarose gels to visually confirm successful amplification of all samples and confirm negative outcomes for all negative controls. Cleanup of PCR products and concentration normalization was carried out using the SequalPrep Normalization Plate Kit (Thermo Fisher Scientific). Sets of 95 amplicons were sequenced at the University of Missouri Genomics Technology Core (https://dnacore.missouri.edu/) using the Illumina MiSeq v2 platform with 150 bp paired- end reads. Field negative controls were included along with field samples for Illumina DNA sequencing to establish baseline expectations for contamination levels in sample data, as well as to establish minimum detection criteria for validation of species detection in field-collected water samples (Klymus et al., 2017; Valentini et al., 2016). An exotic species sample positive control (community tank, PetSmart, Rolla MO) was included on each 96-well plate of samples (Klymus et al., 2017; Nakagawa et al., 2018) to monitor for cross-contamination among samples.

2.4. Bioinformatics

Raw sequencing reads were processed using the Barque (v1.7.3) metabarcoding analysis pipeline (https://github.com/enormandeau/barque), which included steps to remove primer sequences, merge forward and reverse reads, and filter chimeric sequences (Mathon et al., 2021). Taxonomic identification was performed using a reference database consisting of published sequences in GenBank. This database was curated to include only species known to occur in the lower- Mississippi drainage basin to control for erroneous identifications (Schramm et al., 2016). Sequences were identified using a 97% sequence similarity between merged sequences and reference database sequences as a threshold for species assignment and a 95% sequence similarity to assign sequences to a genus level (Deiner et al., 2017). We corrected read counts by subtracting the number of reads detected in negative controls from sample read counts, and we applied a five-read-minimum threshold (range 0.003%–0.5% of reads per sample) using MIN_HITS_SAMPLE (Barque pipeline) for identification of taxon-presence in each sample, similar to thresholds considered for other studies focused on ecological diversity (Alberdi et al., 2018; McColl-Gausden et al., 2021).

2.5. Estimating Species Diversity and Assemblage Structure

Our estimates of species richness were based on detection (presence/absence) of each species at each sample site. To determine presence/absence, we combined 12s and 16s reads for each water sample, and we combined three water samples at each sample site (hereafter combined samples). The number of combined samples ranged from three to seven in each habitat (Figure 1). We treated each habitat as a separate assemblage and calculated coverage, a measure of how completely the assemblage has been sampled (Chao & Jost, 2012) as a function of the number of combined samples. We compared Hill diversity (q = 0-2) estimates standardized to 95% coverage by extrapolation (Chao et al., 2014) to account for differences in sampling efficiency (Roswell et al., 2021) and we produced 95% confidence intervals for estimates of richness with 500 bootstrap replicates. Hill numbers included the raw number of species observed or “richness” (q = 0), an estimate of Shannon diversity (q = 1), and an estimate of Simpson diversity (q = 2), which increasingly down-weight the impact of rare species. Calculations of richness were conducted using iNEXT (Hsieh et al., 2024).

Assemblage relationships among the habitats and seasons were visualized and tested to determine whether significant structure among combined samples was present. Spring assemblage relationships were first examined alone due to our inability to resample sites outside of the bayous in the fall (i.e., water levels were too low). A Sorensen similarity matrix was generated for the spring- only incidence data, and non-metric multi-dimensional scaling (nMDS) was used to visualize community relationships in two-dimensions. PERMANOVA was then run using the same Sorensen matrix, to test for a significant effect of ‘habitat’ among all spring combined samples. The model included ‘habitat’ as a fixed effect, was based on Type III sum of squares, and used unrestricted permutation of raw data (9999 permutations). We followed this model with a pairwise PERMANOVA to determine which ‘habitats’ differed significantly in assemblage composition (assuming ‘habitat’ was found to be a significant effect in the global model). Pairwise PERMANOVA was run using the same parameters as the global model; however, output included P(perm) as well as Monte Carlo P(MC) significance values. The latter was included to provide better clarity of significance when comparing ‘habitats’ containing small sample sizes.

To examine potential effects of habitat and season on assemblage relationships, we limited our dataset to only include the three bayou habitats that were sampled during both spring and fall. We again used Sorensen similarity to quantify relationships among combined samples, and non-metric multi-dimensional scaling (nMDS) to visualize assemblage relationships. PERMANOVA was used to test for significant effect of ‘habitat’ and or ‘season’ among samples. The model included ‘habitat’ and ‘season’ as fixed effects, and the potential ‘habitat’ * ‘season’ interaction. All other model parameters were the same as the previous global PERMANOVA model. In the case of a significant main effect in the global model, pairwise PERMANOVA was then used to determine which specific groups differed (following the same parameters of the previous pairwise PERMANOVA).

Similarity percentage (SIMPER) was used to identify the species driving differences between habitat groups in the spring-only dataset, and between habitat and season groups in the bayou spring and fall dataset. SIMPER analyses were run using the Sorensen similarity matrices previously generated for each dataset. To identify species and environmental factors most correlated with the full Sorensen matrices (i.e., variables correlated with the full matrix rather than subset group comparisons), we used the BIOENV matching procedure based on spearman rank correlation. This procedure identifies individual variables that best match the overall assemblage structure pattern via matrix comparisons. Sorensen similarity matrix comparisons were used to identify the best correlated species, and Sorensen similarity was compared with Euclidean distances among environmental variables to identify the best correlated environmental factors. Environmental variables were normalized to make all variable units comparable prior to generating Euclidean matrices. The BIOENV procedure was run for each of our two datasets. All assemblage analyses (e.g., nMDS, SIMPER, BIOENV) were conducted using PRIMER version 7.0.23, and all PERMANOVA models were run using PERMANOVA+ in PRIMER version 7 (Anderson et al., 2008).

3. Results

In spring 2022 we collected 84 water samples from 28 locations across bayou, shallow lake, slough, and ditch habitats (Figure 1). The volume of water we were able to filter before filters were clogged was dictated by presence of algal growth and suspended sediment load. In the spring, the mean filtered volume per sample was 298 ml (range: 150 – 500 ml). In our fall collection we were able to resample 18 of the locations visited in the spring (with remaining sites being dry). The volume of water filtered per sample from the bayou habitats in the fall was similar to spring, with a mean of 277 ml. However, the fall samples collected in Robinson Lake and the slough were highly turbid in the rapidly receding remnant habitats. We were only able to filter 15-20 ml of water volume per sample before filters were clogged. The third combined sample from Wolf Bayou (Figure 1, Wolf Bayou site 3) was in a small tributary creek far from the bayou. This sample location was determined post hoc to be a separate habitat from the rest of the samples in Wolf Bayou and was excluded from all richness and assemblage structure analyses.

The average Illumina sequence read count per water sample was 77,429 for the 12s marker, and 79,175 for the 16s marker. All field and lab negative controls were visually determined to be negative by gel electrophoresis. Species presence/absence was determined based on an average of 400,184 reads (range: 134,947-687,843) per combined sample. A total of 51 taxa (47 resolved species, three unresolved species-pairs, and one unresolved species-trio at the genus-level) were detected, representing 34 genera and 19 families, summarized in Table 2. Four species clusters remained unresolved due to highly similar reference sequences. In total, 45 out of 51 species were detected by both 12s and 16s markers, with three species each being detected by only the 12s or 16s marker. On average across all combined samples, 66% (range: 48-79%) of species detections were confirmed by both markers with the remaining detections based on only one or the other marker. A combined sample-level record of presence/absence of all taxa detected in the study is provided in Table S1.

3.1. Species Richness

The number of species detected per habitat ranged from 34 (ditch) to 42 (both Wolf Bayou and Robinson Lake). However, habitats varied substantially in size and complexity, and the number of combined samples ranged from three to seven among habitats. Sample species coverage was very good in each of the habitats with estimates ranging between 95% and 99%. For purposes of drawing direct comparisons among habitats and between sampling dates, all species richness estimates were standardized to 95% coverage using rarefaction. After standardization, we found that the slough, Robinson Lake, and Wolf Bayou all exhibited similar richness on the spring sampling date, with the other two bayous and the ditch forming a second group with lower richness (Figure 3). In the bayous, fish assemblages were slightly more diverse on the fall sampling date than in the spring, and the seasonal difference was greater in Samples and Hosner Bayous, than in Wolf Bayou (Figure 3).

We did not quantify fall richness in the lake, slough, or ditch habitats, or include the fall samples from these habitats in comparisons of richness between seasons due to high turbidity, low water volume filtration, and the smaller number of sites that could be sampled. Fewer species were detected in the fall samples from Robinson Lake and the slough, but notably, there were no new species detected in fall samples of these habitats that were not also detected in the more robust spring samples.

Centrarchidae was the most represented family, with 13 species (25% of all taxa) (Table 1). Percidae was represented by five darter species, and several other families were represented by four species (e.g., Ictaluridae, Leuciscidae, Lepisosteidae). Catostomidae would likely have been comparable in taxonomic diversity to other families if species-level resolution had been achieved for the carpsuckers (Carpiodes) and buffaloes (Ictiobus), all of which are present in the region (Pflieger, 1997; Schramm et al., 2016) and previously detected on Black Island (C. Rice pers. comm.). Nine families were represented by single species. Over half of species (28 out of 51) were ubiquitous, which we defined as being detected in at least 75% of combined samples, and present in all habitats. These included six of the sunfish species (Lepomis spp.), largemouth bass (Micropterus nigricans), both crappies (Pomoxis spp.), the buffalo species (Ictiobus spp.), the topminnow species (Fundulus spp.), two catfish species (Ictalurus punctatus and Pylodictis olivaris), three of the gar species (Lepisosteus spp.), the eyetail bowfin (Amia ocellicauda), one darter species (Etheostoma asprigene), and the freshwater drum (Aplodinotus grunniens). All five invasive carp species in the families Cyprinidae and Xenocyprididae were also ubiquitous, and although drawing inference of relative abundance or overall biomass based on proportions of Illumina read counts is complicated by technical and biological factors (Beng & Corlett, 2020; Mathieu et al., 2020), it is notable that nearly 40% of all assigned sequence reads were attributable to just these five invasive carp species, with more than 29% of all reads in the study attributable to the bighead and silver carp (Hypophthalmichthys spp.) alone. The only two native species that rivaled the invasive species with the highest total read counts were the bluegill (L. macrochirus - 12.28%), and the combined buffalo species (Ictiobus spp. - 10.36%).

3.2. Assemblage Structure among Spring Samples

Substantial assemblage structure among spring samples was visualized in two dimensions using an nMDS of Sorensen similarities among all habitats (Figure 4). PERMANOVA confirmed significant differences in assemblages across habitats (df = 5, Pseudo-F = 6.89, p <0.001, Table 3), and post-hoc pairwise PERMANOVA confirmed that all habitats differed significantly except for Hosner and Samples Bayous, which were largely overlapping in nMDS space (t = 0.81, P(perm) = 0.86, P(MC) = 0.61) (Table S2).

Species most correlated with the overall assemblage structure (resultant of BIOENV) included six species with correlation coefficient values between 0.402 and 0.629 (Table S3). Species that were almost exclusively detected in the bayou habitats included the brook silverside (Labidesthes sicculus r = 0.629, Figure 5a), and spotted bass (Micropterus punctulatus r = 0.603, Figure 5b). Species that were largely or completely undetected in the bayous, but detected in most of the lake, slough and ditch samples were the pirate perch (Aphredoderus gibbosus r = 0.441, Figure 5c) and inland silverside (Menidia beryllina r = 0.402, Figure 5d). A species that was ubiquitous among all habitats but absent from the expansive shallow reach of the flooded Robinson Lake was the flathead catfish (Pylodictus olivaris r = 0.450, Figure 5e). The white bass (Morone chrysops) was detected in all slough and ditch habitats but detected in only the deepest portions of Robinson Lake and the bayous (r = 0.409, Figure 5f). Moreover, SIMPER identified species most responsible for pairwise habitat distinctions; however, output from the analysis generally mirrored the overall trends determined through the BIOENV procedure. BIOENV of environmental variables (Table 4) identified conductivity (r = 0.586) and water clarity (Secchi depth r = 0.384) as the two variables most correlated with the assemblage similarity matrix (Figure 6). Conductivity was highest in the ditch, slough, and lake habitats, and uniformly low in the bayou habitats, while conversely, water clarity was greatest in the bayou habitats (Table 4).

3.3. Seasonal and Spatial Variation in Bayou Habitats

PERMANOVA including the three bayou habitats and both seasons indicated significant assemblage differences among bayou habitats (df = 2, Pseudo-F = 2.81, p = 0.0025), but not between seasons (df = 1, Pseudo-F = 2.24, p = 0.055; Table 5). Additionally, the interaction was not a significant effect, and since dropping the interaction from the model did not change any outcomes, the interaction is included in the final model presented. Pairwise PERMANOVA confirmed that differences between Samples and Hosner Bayous versus Wolf Bayou held when both spring and fall data were included (Table 6), the same outcome was found when only spring samples were considered. SIMPER analysis identified paddlefish (Polyodon spathula) and threadfin shad (Dorosoma petenense) as the two main species driving assemblage differences between Wolf Bayou and other bayous. Both species were more frequently detected in Wolf Bayou (P. spathula average frequency per site; SB = 0.13, HB = 0.0, WB = 0.75. D. petenense average frequency per site: SB = 0.25, HB = 0.5, WB = 0.83).

Though differences in assemblage composition between spring and fall were not significant, the fall samples did exhibit marginally higher species richness. One species that was responsible for an increase in fall versus spring richness in Hosner Bayou was the threadfin shad (D. petenense). This species was frequently detected in Wolf Bayou but not the other two bayous in the spring, and then frequently detected in both Wolf and Hosner Bayous, but not Samples Bayou in the fall. Another species that was slightly more abundant in the fall than in the spring was the western mosquito fish (Gambusia affinis), though this was only the case in Samples Bayou and not the other two bayous.

3.4. Comparison to Historical and Recent Capture-Based Records of Species Presence

We obtained a compilation of the MDC historical fish capture-based survey records for backwater habitats on Black Island with sampling dates from 1940, 1966, 1979, and 1994 (C. Rice Pers. Comm.), as well as a fyke net survey of Black Island habitats from the early spring of 2022 that included bayou and lake habitats (D. Ostendorf pers. comm.) (Table 2). The combined capture-based record included 40 species, seven of which were taxonomically unresolved by eDNA metabarcoding (Carpiodes (2), Ictiobus (3), and Fundulus (2)). There were four species in the MDC survey records that were not detected in the eDNA survey. Two were represented by single specimens in the full historical record (Hybognathus hayi and Moxostoma macrolepidotum), and two were absent from our reference database for eDNA detection (H. hayi and Paranotropis shumardi). The only remaining species from the historical record of the Black Island complex that was not detected in our eDNA metabarcoding survey was the flier (Centrarchus macropterus), a species of conservation concern listed by the state as ‘vulnerable’ (MDC, 2024) and not recorded on Black Island after 1940. The spring 2022 fyke net survey was conducted on two dates, with a total of 25 sets, and identified 17 species in total (D. Ostendorf, pers. comm.), all of which were detected in our eDNA metabarcoding survey (Table 2).

Our eDNA metabarcoding survey detected eight native species that were not reported in the MDC survey records for Black Island. These included three sunfishes (L. microlophus, L. miniatus, L. symmetricus), logperch (Percina caprodes), American eel (Anguilla rostrata), and two gar species, the spotted gar (Lepisosteus oculatus), and the alligator gar (Atractosteus spatula). All were relatively rare among eDNA samples except for spotted gar. All four invasive carp species in the Xenocyprididae were unreported in the MDC survey records between 1940 and 1994, with only H. molitrix being recorded in the 2022 fyke net survey.

4. Discussion

4.1. Species Diversity

Quantifying the diversity and distribution of species across temporal and spatial scales is crucial towards a greater understanding of ecological community structure and dynamics (e.g. Geheber & Piller, 2012; Lee et al., 2024; Zbinden et al., 2022). The use of eDNA metabarcoding to assess fish biodiversity and assemblage structure has become common in freshwater (Cilleros et al., 2019; Euclide et al., 2021; Hallam et al., 2021) and marine (Miya, 2022; Schenekar, 2023) environments. We assessed alpha (within habitat) diversity and beta (between habitats) diversity among Mississippi River bottomland habitats, and quantified spatio-temporal fish assemblage relationships. Although habitat complexity makes these systems difficult to thoroughly sample using traditional capture-based methods (Schramm et al., 2016), our use of eDNA metabarcoding allowed a more complete characterization of the diversity harbored by these unique habitats than was previously available.

Our eDNA metabarcoding survey detected 51 out of 134 species listed in a recently compiled comprehensive list of species for the entire lower reach of the Mississippi River (Schramm et al., 2016). We detected eighty-three percent of species (33 of 40) recorded in MDC survey records of Black Island Conservation Area between 1940 and 1994 and tripled the number of species detected in a fyke net survey conducted by the MDC in 2022. Moreover, 19 of 51 taxa we detected using eDNA metabarcoding were never recorded in any previous surveys of Black Island habitats. Our eDNA metabarcoding survey increased the list of known species in the Black Island Conservation Area by thirty-two percent.

Two-thirds of the species detected from their eDNA (34 of 51) were classified as “uncommon” in the Mississippi River and backwaters by Schramm et al. (2016). These results suggest that eDNA metabarcoding was either effective at detecting rare species, or those species may be more common than capture-based surveys generally indicate. There were also seven species detected from their eDNA in this study that are listed as “common” or “abundant” in Mississippi River habitats (Schramm et al., 2016), but have never been reported in Black Island capture-based surveys. Some of these species reach large adult size (e.g. I. punctatus, P. olivaris, A. grunniens) that may make them more challenging to record by capture-based means using methods commonly employed in backwater habitats (e.g. fyke nets). Others are river channel species (e.g. A. chrysocloris, N. atherinoides) (Pflieger, 1997) that may be present in backwater habitats only infrequently. Two-thirds of the species detected by their eDNA are dependent on backwater habitats for at least a portion of their life history (Table 2), while only a third of the species are dependent on fluvial habitats for any portion of their life histories (Schramm et al., 2016).

This study demonstrated that Black Island lowland habitats contain a highly diverse ichthyofauna. Estimates of species richness were comparable to or greater than eDNA-based species richness estimates reported in Missouri Ozark uplands stream habitats, which were conducted using the same set of 12S and 16S gene markers (Lee et al., 2024). In a survey of three sites in the St. Francis River drainage within the Ozark-lowland border region of Missouri, species richness (Hill number q = 0 and 95% species coverage) estimates ranged ~25-35 species (Lee et al., 2024) compared to estimates of 32-40 species per Black Island habitat (Figure 3) in this study.

4.2. Assemblage Structure

Variations in backwater habitat characteristics that include vegetation, submerged structure, depth, seasonal persistence, and connectivity to the river channel impact fish assemblages among Mississippi backwater habitats (Dembkowski & Miranda, 2012; Koel, 2004). On Black Island, deeply incised bayou habitats were the most isolated from the Mississippi River channel, and exhibited the greatest depth, steepest banks, highest clarity, lowest conductivity, and the greatest seasonal stability among all sampled habitats. Robinson Lake, the slough, and the ditch were all more directly connected to the river channel, particularly whenever the river was at flood stage (e.g. in spring 2022, Figure 2). These habitats were seasonally more dynamic regarding depth, surface area, and temperature, with many of the spring sample sites presented as dry on the fall sampling date. The contour was very flat in the lake, with aquatic habitats inundated into heavily vegetated surrounding areas in the spring. Conductivity and Secchi depth (a proxy for clarity) were identified as the two environmental factors most correlated with assemblage structure; specifically, conductivity was lower and clarity was higher in bayous compared to all other habitats. The reduced conductivity and greater clarity of the bayous may result from isolation from the river channel in the spring, though this was not tested.

Generally, the overlap in species occupying various habitats was extensive. However, each of the habitat-types had distinctive fish assemblages, as indicated by nMDS and PERMANOVA. The relatively small subset of species identified as main contributors to the assemblage differences among habitats included two silverside species (family Atherinopsidae), with the brook silverside (L. sicculus) frequenting the bayou complex, and the inland silverside (M. beryllina) being more frequent in Robinson Lake as well as the interconnected ditch and slough. Brook silversides are common in both upland Ozark and small lowland stream habitats throughout the region, while inland silversides are a large-river species, common in the Mississippi River (Pflieger, 1997); so these contrasting distributions reflect differences in connectivity between Black Island habitats and the Mississippi River channel. The western pirate perch (Aphredoderus gibbosus) exhibited a distribution among habitats very similar to the inland silverside, and was detected frequently in the lake, ditch, and slough habitats, but was nearly absent from the bayous. Western pirate perch are widely distributed among lowland habitats and also found in the Mississippi River channel and backwaters (Pflieger, 1997). The distribution of the spotted bass (Micropterus punctulatus) matched that of the brook silverside in the bayous, where it may do well in the relatively clear lentic habitat.

We examined spring and fall bayou samples to assess differences between the three bayous as well as the potential seasonal effect on assemblage composition. Significant assemblage structure within the bayou complex was more subtle and primarily contrasted the larger and deeper Wolf Bayou with the narrower and shallower Samples and Hosner Bayous. Differences between these respective bayou habitats were primarily due to the greater diversity of species detected in Wolf Bayou, including the paddlefish (Polydon spatula), and the threadfin shad (Dorosoma petenense). Both are filter feeders, preferring larger open waters (Pflieger, 1997). Though the threadfin shad was only detected in Wolf Bayou in the spring where it was regularly detected, it was also detected frequently in Hosner Bayou and rarely in Samples Bayou in the fall. This pattern of distributional expansion could be attributable to an increase in seasonal abundance due to spawning between the sampling dates in May versus October.

4.3. Taxonomic Distribution of Diversity

Across all families and in all samples, the greatest taxonomic diversity was exhibited by the sunfishes and black basses (family Centrarchidae). Members of this family generally exhibit sedentary habits and are distributed often in slow moving or still waters with submerged cover (Pflieger, 1997), so Mississippi River lowlands provide ideal habitat for many of these species. In fact, nine of the thirteen extant sunfishes in the genus Lepomis (Harris et al., 2005; Near & Koppelman, 2009) were detected in this study. Five of the sunfish species have widespread distributions, while four of the species are rare or exhibit distributions specific to lowland habitats, making their presence notable. These included the dollar sunfish (L. marginatus) which was ubiquitous according to eDNA data, and the redear (L. microlophus), redspotted (L. miniatus), and bantam (L. symmetricus) sunfishes, each of which were detected in only a few samples in either the bayous or Robinson Lake (Table 2).

Equally important as the families that were well represented in the Black Island habitats is a family that was not well represented. There was notably limited diversity of North American minnow species (family Leuciscidae). The Leuciscids are the most speciose fish family in North America (Schönhuth et al., 2018) as well as in the Mississippi River drainage (Schramm et al., 2016). Only four minnow species were detected in our study, and of those, only the pugnose minnow (Opsopoeodus emiliae), was widespread across habitat types. The Mississippi silvery minnow (Hybognathus nuchalis), emerald shiner (Notropis atherinoides), and golden shiner (Notemigonous crysoleucas) were each detected in only a few samples. The minnow family includes numerous large-river species that have been documented in the adjacent Mississippi River (Pflieger 1997, Schramm et al. 2016), but were not detected in this or previous surveys of the Black Island wetland complex, suggesting that these wetland habitats may not be appropriate for these species, at least at the times and places the surveys were conducted. There are also a number of minnow species whose distributions include many of the small streams of the adjacent alluvial plain ecosystem that were not detected in this study (Pflieger, 1997). The absence of these species could be indicative of the impact of the levee system as a barrier preventing access by some wetland species to bottomland habitats trapped along the margins of the Mississippi River by the levee system. A more comprehensive evaluation of this possibility provides a future avenue for conservation research in this system.

4.4. Invasive and Introduced Species, and Species of Conservation Concern

Our results indicated that eDNA metabarcoding is a useful tool for detecting invasive species as demonstrated in previous studies (Pukk et al., 2021). All five invasive Asian carp were detected ubiquitously across all habitats we sampled. They were also frequently observed, particularly in the shallow slough and Robinson Lake habitats in both the spring and fall sampling trips. The fact that forty percent of all eDNA reads in the study were attributable to these five invasive species is an alarming observation, and underlines the scale of invasive species prevalence and impacts in these lowland habitats. The sheer volume of eDNA detected from these invasive species, and particularly from bighead carp (Hypophthalmichthys nobilis) and silver carp (H. molitrix), which together accounted for 29% of all sequenced eDNA reads, does require caution in making inferences regarding when and where these species may be frequenting within the Black Island habitat complex. Studies have shown that eDNA that becomes bound to silt and clay particles can have persistence times that are orders of magnitude longer than eDNA in the water column (Nevers et al. 2020, Turner et al. 2015). It is possible that high abundance, high biomass species, like the Hypophthalmichthys species, have saturated portions of the floodplain environment with eDNA, making their detection unavoidable. If the goal of an eDNA metabarcoding survey is to determine the distribution of the less abundant native fish species in a region, it may be beneficial to modify the detection conditions by using blocking primers for Hypophthalmichthys species in the PCR reactions to avoid these common invasive species detections (Rojahn et al., 2021).

Historically, the lowland habitats of Black Island provided ideal habitat for the critically imperiled alligator gar (Atractosteus spatula) (Adams et al., 2024; Smith et al., 2020; Solomon et al., 2013). Wolf Bayou was identified as a priority site for reintroduction of alligator gar by the Missouri Department of Conservation, and the species was stocked five different times between 2013 and 2022 (S. Mondragon pers. comm.). It was reassuring to detect alligator gar eDNA at three sample locations in Wolf Bayou, including both spring and fall samples. Not much more can be inferred about the status of alligator gar in this bayou complex other than that the species was not detected in either of the other bayous, which are all interconnected. Other species of conservation concern detected in this survey included the imperiled bantam sunfish (Lepomis symmetricus) in both bayou and Robinson Lake habitats, and the vulnerable Mississippi silvery minnow and river darter, both of which were detected in the river slough habitat.

5. Conclusions

Though the degradation of floodplain habitats persist in modern times (Morrison et al., 2023), our study further demonstrates the wealth of information that can be extracted from the advancement of eDNA biomonitoring tools. Importantly, we were able to demonstrate significant fish assemblage distinctions among lowland habitat types using eDNA survey methods across a period when a fyke net survey identified far fewer species only months earlier. Species of conservation concern were documented, and fish assemblages associated with water permanence and river connection took shape. Information from this study and others stemming from it will contribute to the understanding of floodplain communities and how to best preserve them.