Diversity and composition of macroinvertebrate communities in a rare inland salt marsh

Abstract Inland salt marshes are rare habitats in the Great Lakes region of North America, formed on salt deposits from the Silurian period. These patchy habitats are abiotically stressful for the freshwater invertebrates that live there, and provide an opportunity to study the relationship between stress and diversity. We used morphological and COI metabarcoding data to assess changes in diversity and composition across both space (a transect from the salt seep to an adjacent freshwater area) and time (three sampling seasons). Richness was significantly lower at the seep site with both datatypes, while metabarcoding data additionally showed reduced richness at the freshwater transect end, consistent with a pattern where intermediate levels of stress show higher diversity. We found complementary, rather than redundant, patterns of community composition using the two datatypes: not all taxa were equally sequenced with the metabarcoding protocol. We identified taxa that are abundant at the salt seep of the marsh, including biting midges (Culicoides) and ostracods (Heterocypris). We conclude that (as found in other studies) molecular and morphological work should be used in tandem to identify the biodiversity in this rare habitat. Additionally, salinity may be a driver of community membership in this system, though further ecological research is needed to rule out alternate hypotheses.

The intermediate disturbance hypothesis (IDH) predicts that highest diversity will occur in areas with intermediate levels of disturbance: at high levels of disturbance, only the most disturbancetolerant species can persist, while at low levels of disturbance, strong biotic interactions dominate and reduce species diversity (Connell, 1978). The IDH model can be expanded to include abiotic stress more generally (e.g., salinity stress), whereby high levels of stress reduce the number of species that can live in an area, and more clement conditions may lead to strong biotic interactions.
Although the utility of the IDH as a theoretical framework has been questioned (Fox, 2013), as a verbal model it anticipates that diversity and stress may not be directly, linearly related.

| Inland salt marshes
Although salt marshes are generally associated with coastal environments, there are several different types of inland saline habitats (e.g., the Great Salt Lake in the United States, temporary wetlands in the Camargue region of France; Herbst, 2001;Waterkeyn et al., 2008).

Inland salt marshes surrounding the Great Lakes region of North
America are rare, patchy habitats. They are formed on salt deposits from an ocean that covered the region approximately 400 million years ago, during the Silurian period. In some places, glaciation has removed overlying sediment and groundwater now interacts with the salt (Albert, 2010;Eallonardo & Leopold, 2014;Lincoln et al., 2020). Although never an abundant habitat, inland salt marshes have been very degraded through salt extraction for commercial use, beginning during the settler/colonial period in what is now the United States. Marshes were often chosen as settlement sites due to the role of salt as a critical element in food preservation (Lincoln et al., 2020). We hereafter use "inland salt marshes" to refer exclusively to the habitats that are from this geological origin and which are located generally in the Great Lakes region of North America; we do not further address other types of inland saline habitat.
Inland salt marshes are disconnected from the ocean and exist as patchy habitats in the landscape, designated as critically endangered both on a global scale and within the state of Michigan (Lincoln et al., 2020). The salt constitutes a source of abiotic stress for the species that live there, which are either freshwater species or brackish-water specialists that have evolved from freshwater ancestors. The plant communities that have been documented at two salt marshes in Michigan and New York, for example, include halophilic species (Eallonardo & Leopod, 2014). Although the plants that have been identified in these marshes are different from those in coastal salt marshes, there is still a clear zonation of plants that can tolerate these salty conditions versus those that cannot (Eallonardo & Leopold, 2014;Lincoln et al., 2020). However, very little is known about the invertebrate species that live in these habitats (Albert, 2010), including basic species composition.
Another potential source of stress in inland salt marshes is the seasonality of the water level. During wetter seasons, the habitat is more suited for animals that require standing water, and in dry seasons, the water in the marsh can dry up altogether. The variability of the water level in ponds and other wetlands has been shown to affect the diversity and composition of plants (Riis & Hawes, 2002) and animals that live in these areas (Anderson & Smith, 2000;Strachan et al., 2014;Uzarski et al., 2004).

| Metabarcoding as a technique
One of the problems with characterizing diversity in small invertebrate species (such as those found in wetland sediments) is the difficulty in identifying the species using morphology, which requires a high degree of time and expertise to identify individuals to the species level. Additionally, morphological identification at the species level often cannot be completed for subadults in a population.
Metabarcoding uses high-throughput DNA sequencing and bioinformatics to identify all members of a community that are present in a sample. Although there are disadvantages to this method (e.g., biases in primers or other steps of the processing protocol, incompleteness of reference databases; Carugati et al., 2015;Danavaro et al., 2016;Cahill et al., 2018), it allows for rapid identification of small organisms, potentially to the species level where such a precise ID is not possible based on morphology (e.g., cryptic species or juvenile life stages; Danavaro et al., 2016;Pearman et al., 2016). The two techniques are complementary rather than redundant: metabarcoding can be inconsistently successful across taxa (Afzali et al., 2020;Cahill et al., 2018;Kelly et al., 2017;Lejzerowicz et al., 2015) while simultaneously allowing more precise and thorough identifications.
The combination of morphological and molecular tools is therefore important for verifying the conclusions of a study where possible.

| Hypotheses
In this study, we investigated the relationship between abiotic stress (salinity) and diversity in an inland salt marsh. We hypothesized that species diversity and richness would vary with changing salinity in space, though we did not have strong a priori expectations for the shape of the relationship. For instance, diversity/richness might increase with decreasing abiotic stress or might peak at intermediate stress levels. We also measured seasonal change in diversity, richness, and community composition and expected to find temporal changes in these variables. We used both morphological identification and molecular techniques (barcoding and metabarcoding) to address these questions, expecting that metabarcoding would give us the most thorough picture of community composition in the marsh.

| Study site
The study was conducted in the Maple River inland salt marsh, part of the protected Maple River State Game Area in Fowler,Michigan (USA;43.0847 N,84.7649 W). The marsh occupies 6.5 acres, with a small halozone (Lincoln et al., 2020;Figure 1), and the plants present in the area affected by the salt seep have been described by the Michigan Natural Features Inventory (Albert, 2010;Lincoln et al., 2020). Although the seep is within the larger floodplain of the Maple River, it appears to be unaffected by seasonal flooding of the river, with water levels primarily changing due to changes in precipitation and groundwater seepage (Lincoln et al., 2020).

| Field sampling
In April, July, and October of 2018, we sampled along a transect within the marsh. The first sample point on the transect was the seep site within the marsh, identifiable because of its lack of vegetation (the rest of the marsh is a dense stand of predominantly Typha spp.). We sampled at 7 points along the transect, spaced 20 m apart, for a total transect length of ~120 m. The transect ran parallel to the ridge that forms the boundary of the floodplain (Figure 1). A transect length of 120 m was sufficient to get out of the area influenced by the salt seep: it took us across a small berm, out of the Typha marsh, and into an area with different vegetation (including some trees). Three sediment samples were taken from the surface at each transect point, each approximately 450 ml. We used surface sampling because a pilot look at sediment cores showed that invertebrates were only present in the top layer of the marsh. Samples were returned to the laboratory in Albion, Michigan, where half of each sample was preserved at −80°C for molecular analysis and half was preserved in 70% ethanol for morphological identifications.
Soil salinity was measured in the laboratory following the protocol of Rhoades (1996). This involved drying the soil sample, weighing ten grams of each sample, suspending it in 100 ml of water, agitating it every ten minutes for an hour, and then measuring total salinity (in psu) with a YSI probe (model Pro2030). Since the organisms in the sample were living in the porewater of the sediment, soil salinity is an ecologically relevant measure of the environment.

| Morphological data collection
Each sediment sample was filtered at 100 µm and rinsed with DI water, then placed in a sorting tray and examined under a dissecting microscope. Animals that were larger than ~1 mm were removed and identified using a dichotomous key (University of Wisconsin Extension, 2014). Identifications were done to the lowest taxonomic level that could be easily and reliably reached; this was often the class or family level but varied among taxonomic groups. The complete list of taxa used in morphological identifications is available in Appendix S1. At least two replicate samples were sorted per transect point, totaling at least 14 samples at each of the three sampling timepoints (44 samples overall).

| Analyses: Morphology
We calculated both rarefied taxon richness and Simpson's diversity index for each sample and used two-way analyses of variance to F I G U R E 1 Sampling site in the Maple River salt marsh, after Lincoln et al. (2020). Stars indicate the salt seep (left) and freshwater (right) ends of the transect; the total length of the transect is 120 m, and samples were taken every 20 m. The transect begins at (43.0847N, 84.7649W) and ends at (43.0852N, 84.7642W). White indicates the area affected by the salt seep compare these indices across season and across the transect points (hereafter "sites"). We used permutational analysis of variance (PERMANOVA) and nonmetric multidimensional scaling (NMDS) plots to analyze and compare community composition across space and time. All analyses were conducted in R, version 3.5.3 (R Core Team, 2019), with the package vegan (version 2.5-5, Oksanen et al., 2019).

| COI barcoding of target specimens
Two species of arthropod were consistently found in high population densities at the salt seep: an ostracod and a ceratopogonid midge larva. To get a more precise species identification, we performed COI barcoding. We extracted DNA from a one midge larva and one ostracod using Qiagen DNeasy kits and amplified a ~650-bp fragment of cytochrome oxidase I using universal LCO and HCO primers (Folmer et al., 1994).

| Metabarcoding
For metabarcoding, we extracted DNA from approximately 10 g of each soil sample using Qiagen DNeasy PowerMax Soil Kit. We used all samples collected or 21 per timepoint (63 total samples). Samples from each timepoint were extracted and sequenced separately (Table 1). We purified extracted genomic DNA using Qiagen DNeasy PowerClean Cleanup Kit.
We used universal metazoan primers mlCOIint and HCO2198 from Leray et al. (2013) to amplify a fragment of the COI gene.

| Analyses: Metabarcoding
We analyzed the metabarcoding data using MOTHUR version 1.40.5 (Schloss et al., 2009), with the MiSeq standard operating procedure described in Kozich et al. (2013). Sequences from all three sampling timepoints were analyzed in a single run of the bioinformatic pipeline. The paired ends were joined with make.
contigs. We used screen.seqs (maxlength = 370, maxhomop = 8, maxambig = 0) for quality filtering and then produced unique sequences using unique.seqs. The reference sequences produced in this step were aligned using align.seqs and a reference from the Barcode of Life Database (BOLD; Ratnasingham & Herbert, 2007). Pre-clustering was undertaken with pre.cluster (diffs = 3) and chimeras removed using vsearch (Edgar, 2010). Sequences were identified using the classify.seqs command from MOTHUR v.1.39.5 implemented in Galaxy, with the public server at usegalaxy.org (Afgan et al., 2018). We compared our sequences to BOLD for identification.
Using all OTUs (not just ones that were identified), and after removing OTUs with 10 or fewer reads, we conducted a permutational multivariate analysis of variance (PERMANOVA) to analyze community composition across the seven different sites. We also conducted a separate PERMANOVA comparing the three sampling seasons. These PERMANOVAs were followed with nonmetric multidimensional scaling analyses (NMDS) based on Bray-Curtis distances. All data were fourth-root transformed prior to these analyses. We also calculated both rarefied taxon richness and To look at taxon-specific patterns of diversity and composition, we used the subset of reads that could be identified based on the BOLD database. After removing reads that were not identified or that belonged to non-animal taxa, we identified the nine taxa with the greatest number of reads both overall and within Insecta.
We then repeated the analyses (PERMANOVA, NMDS, diversity/ richness calculations) with only the OTUs that were identified as insects.

| Analyses: Comparisons
In order to compare the results from the two different datatypes, we correlated the diversity indices and rarefied richness values that were calculated in each sample with the two types of data. Since not all samples were analyzed for morphology, only those with both metabarcoding and morphological data were included here (N = 44 samples). In order to compare community composition, we compiled a matrix of Bray-Curtis distances for each site pair with each datatype and used a Mantel test to compare the distance matrices between datatypes. These analyses were conducted for each of the three sampling dates combined. For all of these comparative analyses, the morphological dataset was compared with the full metabarcoding dataset with all OTUs.

| Salinity data
The salinity in the marsh ranged from 0.025 to 3.597 psu, depending on the time of year and point on the transect. Mean salinity in April was 1.83 psu, while in July it was 0.43 psu and 0.10 psu in October.
The first point, located in the seep, was consistently the saltiest, and salinity dropped rapidly moving away from the seep (Figure 2; Appendix S5).

| Morphology
We found 25 different taxa using morphology-based identifications.
The most widely distributed taxon was gilled snails (Gastropoda: prosobranch pond snails), which appeared in nearly every sample, but the most abundant taxon in terms of overall numbers was an ostracod (possibly Heterocypris salina, see barcoding results below), which were found primarily at the seep site itself (total number found = 346; Figure 3). Community composition clearly shifted over the course of the year (Figure 3). Ostracods were most abundant in July, while the October samples were dominated by gilled snails.
Worms were more abundant in July and October than in April. The complete morphological dataset is found in Appendix S1.
Whole-community analyses (PERMANOVA) revealed significant effects of both transect point and season, as well as their interaction, on community composition ( Table 2). The NMDS analysis shows an overlap of sampling points, but clearer separation by season ( Figure 4): April and October overlap in composition but July has a larger variance. This was due to a larger number of terrestrial animals that were present in July (when there was much less water in the marsh), including ticks, spiders, collembolans, and earthworms.
There was no effect of season on the rarefied taxon richness (p = .385), but there was a significant effect of site (p = .035). The seep site had lower richness than two points in the center of the transect based on Tukey's HSD post hoc tests ( Figure 5). There was no interaction between site and season (p = .856; Table 3, Figure 5).
There was no significant effect of site, season, or their interaction on Simpson's index of taxonomic diversity (Table 3), though there was a visible trend toward lower diversity at the seep site ( Figure 5).

| COI barcoding
The COI barcoding of the midge larva resulted in only 302 bp of high-quality DNA. This relatively short fragment yielded equal identity matches in BLAST to two species of biting midge (Diptera: Ceratopogonidae): 99.7% to both Culicoides sonorensis and C. variipennis. The ostracod sample resulted in a sequence of 555 bp of high-quality DNA with a BLAST match to Heterocypris salina (Ostracoda: Cyprididae), with a 95.9% identity score. However, between a small fragment size in the midge larva and a relatively low percent identity match in the ostracod, any attempt to identify the individuals to species based on these barcodes is premature. The sequences have been submitted to GenBank, with the accession numbers MZ444699 and MZ444700. Note: For each month (= sequencing run), the number of quality reads retained after bioinformatic processing in MOTHUR, the number of OTUs these reads represent, the mean number of reads per OTU (± SD), and the number of taxonomic categories following identification using the BOLD database. The number of taxonomic categories was calculated after the removal of fungi, algae, and unidentifiable OTUs.

TA B L E 1 Summary of sequencing output
F I G U R E 2 Salinity changes in space and time. The salt seep corresponds to point A, the freshwater end of the transect to point G, and all points were sampled 20 m apart

| Metabarcoding results and community analysis
The bioinformatic pipeline recovered 16 617 unique OTUs from 1 090 805 quality-filtered reads. We conducted composition and diversity analyses on this full dataset as described above. The full set of OTUs and their abundances is available in Appendix S2; the list of identified OTUs and abundances is in Appendix S3. There was variation in the numbers of reads, OTUs, and identified OTUs recovered among months (Table 1, Table S1).
The PERMANOVAs conducted with all OTUs showed that sites were strongly different across the transect (  (Figure 4). There was also differentiation according to season (Table 2), and the NMDS plot showed that July again had a wider variation than April or October, driven by a point at the seep site in July where over 99% of sequences were from a single OTU (later identified as a mite in the order Sarcoptiformes). The interaction of site and season was also significant in the PERMANOVA (Table 2).
Based on all OTUs, there was a significant effect of sampling site on rarefied taxon richness (p < .001), but not of season (p = .095) nor of the interaction between site and season (p = .155; Table 3).
Richness showed a relationship where the seep had lower richness than all other sites (Tukey's HSD tests; Figure 6). Although not statistically significant after corrections for multiple testing, there is a hump-shaped relationship in richness (with richness highest in the center of the transect; Figure 5). Taxon diversity (Simpson's index) was not different based on site or season, nor was there an interaction between the two (Table 3, Figure 6).
The PERMANOVA of only the OTUs identified as Insecta showed a significant effect of both site and season on insect community composition, though no interaction between the two ( Table 2). As with the full OTU dataset, there was a slight difference between the seep community and other sampling points, but the difference among replicates was much less striking than with the full OTU dataset ( Figure S5). In contrast with the full OTU dataset, there was no difference in rarefied richness according to either sampling site or season ( seasons, but was lower at the seep than in the center of the transect ( Figure S6).

| Taxon identifications based on metabarcoding data
In addition to whole-community analyses, we found the animal taxa that were most commonly represented in the dataset. For these analyses, we removed OTUs that were identified as fungus (8954 reads) or algae (26,651 reads), as well as the 54% of reads which were unable to be classified using the BOLD database. This left 108 taxonomic groups from 462,308 reads. Taxonomic identifications ranged from the phylum level (e.g., "Unidentified Arthropoda") to the species level. The most abundant OTUs in all samples were from arthropods, particularly insects (Figure 7, Appendix S3). The most abundant taxonomic classes identified by the analysis pipeline, in decreasing order of abundance, were Insecta (arthropods), unclassified Arthropoda, Clitellata (annelids), Arachnida (arthropods), Hydrozoa (cnidarians),

Bdelloidea (rotifers), Gastropoda (molluscs), Monogononta (rotifers),
and Ostracoda (arthropods). We identified a few interesting patterns in the data: annelid worms were most abundant toward the center of the transect with a peak that shifted through the seasons (probably in response to changing water levels), while ostracods were patchily abundant ( Figure 6). We found many taxa with the metabarcoding that were not identified with morphological techniques; conversely, gastropods, which were commonly found in microscope samples, were relatively rare in the metabarcoding samples (Figure 7).
When we looked at community composition within the insect sequences (and after removing sequences that could not be identified lower than Insecta), the most common taxa in descending order were reads relatively evenly spaced (though rare at the seep site; Figure 8).

| Comparison of methods
Mantel tests revealed a weak relationship between the two datasets (metabarcoding and morphology) when all OTUs were considered and when data were separated based on sampling site (Mantel's r = .438, p = .07; Table 4), but no relationship between methods when data were separated based on season (Mantel's r = −1, p = 1; Table 4).
There was a weak but significant correlation between Simpson's diversity metrics calculated with morphological and molecular data (all OTUs; r = .321, p = .033; Figure 9). There was also a correlation F I G U R E 4 Nonmetric multidimensional scaling plot based on morphological data (a and b) and molecular data (c and d), with points coded according to site along the transect (panels a and c, where A corresponds to the saltwater seep and G to the freshwater end) and season (panels b and d). The stress in panels a and b is 0.244; the stress in panels c and d is 0.181 in rarefied richness between the two methods (r = .473, p = .001; Figure 9).

| D ISCUSS I ON
Our study describes the invertebrate community composition in a rare inland salt marsh over both a spatial gradient and a temporal one, corresponding with salinity stress. We found that both spatial and temporal changes are drivers of the community composition and diversity, but that the details depend on the analysis and datatype considered. As many other studies have found, the information gained by using morphological and molecular data together is complementary, though both richness and diversity were correlated between datasets.
Both community composition and species richness changed with sampling site using the metabarcoding data. The NMDS analysis showed a distinct community at the seep, a very different community at the freshwater end of the transect, and a group of sites in the middle that were not clearly differentiated from each other. Looking more closely at the OTU data revealed that the center of the transect contained both groups that were abundant in the more freshwater sites (e.g., earthworms) and those that were abundant closer to the seep (e.g., ostracods, ceratopogonid midges), and this contributed to high richness values in the transect center. Conversely, NMDS analysis showed strong overlap between the three seasons, but this was mainly due to one outlying point in the July samples. At this point, from the seep site, over 99% of reads were identified as an unclassified species of Sarcoptiformes (a mite). It seems likely that this is due to one large individual that dominated the PCR and sequencing; this in turn decreased diversity and richness metrics at that point.
The relatively low taxonomic precision of the morphological dataset precluded statistically significant patterns of composition, and NMDS did not reveal differentiated communities along the transect. However, richness increased with distance from the seep site, and diversity showed a weak trend in this direction as well. Like the metabarcoding data, there was no effect of season on diversity or richness. The NMDS analysis of the morphological data again showed that July was much more variable in terms of community composition; April and October strongly overlapped with each other.
The data recovered with traditional morphological identifications are of a much coarser scale than the molecular data and did not give clear results on how community composition changes. Additionally, in large part because of how the sediment was processed for these F I G U R E 5 Taxon richness (rarefied richness, panels a, b) and diversity (Simpson's Index, panels c, d) based on morphological data. Panels a & c present the data grouped by site (where A is the saltwater seep and G is the freshwater end of the transect), and panels b & d present the data grouped by sampling time. Letters in panel a represent significant differences according to Tukey's HSD post hoc tests morphological identifications (e.g., selecting animals >1 mm for identification), many taxa could not have been identified. For example, two of the most common taxa in the molecular data were Bdelloidea and Monogononta, classes of Rotifers. These organisms are smaller than our morphological size limit, transparent, and can survive complete desiccation (Caprioli et al., 2004), making it exceedingly unlikely that our morphological IDs would have seen them.
None of these weaknesses of morphological data are unique to our study, and all have been discussed at length elsewhere (e.g., Cahill et al., 2018;Danavaro et al., 2016;Kelly et al., 2017).
However, without the morphological dataset, we would have missed key patterns in the marsh. Most notably, some of the most common invertebrates in the morphological dataset were several groups of gastropods. However, gastropods were much rarer in the metabarcoding dataset (though still in the list of top-ten taxa with the most reads). This result is consistent with other papers, which have found that molluscan DNA is less prevalent in metabarcoding data, including with the primers that we used (Beentjes et al., 2019;Cahill et al., 2018;Kelly et al., 2017). This may relate to the difficulty of crushing the shells of microgastropods during the lysis steps of the DNA extraction process, meaning that the starting DNA pool may contain less molluscan DNA relative to the more easily digested arthropods and annelids. Additionally, the primers we used (Leray et al., 2013) are putatively universal, but seem to work better with arthropods and annelids than molluscs (Cahill et al., 2018).
Another surprising result in the metabarcoding data was the very high number of reads from lepidopterans (butterflies and moths); this was the most abundant order of insects. Lepidopteran OTUs were present in all replicates, usually with a particular OTU being present in only a single replicate. Despite this, we did not see a single lepidopteran larva or adult in the morphological analysis. One possible explanation here is puddling: a behavior that butterflies exhibit in which they dabble in puddles to increase their salt intake. Male butterflies in particular need salt for proper gamete formation (Sculley & Boggs, 1996). Through puddling, butterflies would be expected to leave eDNA that could be collected in our samples. Another possible explanation for the large number of lepidopteran reads is birds that prey on the insects and then excrete butterfly DNA; this seems like a less likely explanation because we had no identified reads from Aves in the dataset (Appendix S3). We also note that other technical sources of error (e.g., primer bias) might be better explanations, as eDNA seems unlikely to overwhelm the results in this way. It is also possible that we did not identify lepidopterans in any of the morphological samples purely by chance and that they are in fact present in the marsh sediments.
Both datasets revealed a large number of ceratopogonid midges (biting midges). DNA barcoding identified the species at the seep itself as either Culicoides sonorensis or C. variipennis, and metabarcoding identified four OTUs belonging to this family, in at least two genera. Both C. sonorensis and C. variipennis are found in brackish and polluted habitats (e.g., coastal salt marshes and liquid manure pools in agricultural facilities; Schmidtmann et al., 2000), so their presence in the Maple River marsh is unsurprising. Culicoides sonorensis is a vector of bluetongue virus, a livestock disease (Maclachlan & Mayo, 2013;Purse et al., 2005). The barcoded individual was a very close match to both species. Culicoides variipennis and C. sonorensis are part of a species complex (Holbrook et al., 2000) and were until recently described as two subspecies of C. variipennis. Although  We conducted these analyses using sites as a rough proxy for salinity, given their increasing distance from the salt seep.

TA B L E 3 Richness and diversity
However, as Figure 2 shows, salinity varied dramatically particularly at the seep site during the year, and there was also seasonal variation in the relationship between site and salinity.
When we regressed rarefied taxon richness against salinity, rather than site, there was no relationship in either the morphological dataset (p = .524) or the metabarcoding dataset (p = .370). There were also no significant relationships between salinity and taxon diversity in either the morphological dataset (p = .549) or the OTU dataset (p = .598). Given these results, and the lack of seasonal signal in community composition based on our NMDS plots, we therefore conclude that although taxa seem to respond to salinity stress in the long-term (i.e., by the development of groups of salt-and fresh-associated taxa), short-term (seasonal) variations in salinity do not strongly reshuffle the community in the marsh, even over the short spatial scale of our transect. Some taxa do show seasonal shifts in their abundance along the transect-this is particularly evident in annelid worms in both the morphological and molecular analyses, which are found in the center of the transect in April and October and move toward the seep during our July sampling. However, the worms may be responding to water level rather than salinity per se; there was no standing water at the seep in July and we saw an incursion of terrestrial fauna at that time (e.g., arachnids, beetles). In fact, a major limitation of our study is that we did not systematically examine other ecological variables, such as water depth, organic content of the sediment, or vegetation cover, and so we cannot include them into our quantitative analyses. Although water depth was spatially and seasonally variable, along the transect it was consistently shallow, never exceeding 1 m during our sampling period. Emergent vegetation is nonexistent at the seep itself, but points farther from the seep are in a Typha-dominated marsh (Lincoln et al., 2020), which should also influence the invertebrate community (though we note that the difference in vegetation across the marsh is a response to salinity; Albert, 2010;Lincoln et al., 2020).
Another possible driver of composition in the marsh is dispersal from the Maple River itself. Although we cannot rule out this idea, other sampling we have done in the floodplain near the edge of the river shows that this community is full of isopods, damselfly larvae, and other species more characteristic of flowing water (data not shown). Additionally, all of our sampling sites are roughly equidistant from the river (Figure 1), meaning that the strong spatial differences we find within the marsh are not the result of only distance to the river. In future work, we hope to explore other potential drivers of these invertebrate communities, including the relative roles of community assembly versus species sorting in this habitat.

| CON CLUS IONS
Our results demonstrate that the biodiversity and composition of invertebrate fauna in this inland salt marsh depend on distance from the salt seep, which may be related to abiotic salinity stress.
Diversity and richness were lower at the seep site, and metabarcoding data supported a hump-shaped relationship where intermediate levels of salinity allow the coexistence of taxa from both ends of the transect. Similar to other studies, we found that molecular data provide higher precision but that not all taxa were sampled equally well, making the continued use of both datatypes essential in this system. We also identified a few taxa that are common at the relatively high-salinity seep, and work is underway to investigate the salt tolerance and adaptation of these species, as well as the role of other ecological factors in determining diversity patterns in space and time in this rare habitat.

ACK N OWLED G M ENTS
We thank Chad Fedewa of the Michigan Department of Natural Activity.

F I G U R E 9
Correlation of morphological and molecular diversity (Simpson's Index, panel a) and richness (rarefied richness, panel b)

CO N FLI C T S O F I NTE R E S T S
The authors have no conflicts of interest. Bach Tran: Data curation (supporting); Investigation (supporting);

AUTH O R CO NTR I B UTI
Writing-review & editing (supporting).

Raw sequence reads have been deposited at the NCBI Sequence
Read Archive with the project number PRJNA741798; COI sequences for the ostracod and midge larva have been deposited in GenBank with accession numbers MZ444699 and MZ444700, respectively. The datasets of OTUs and taxon identification (morphological and molecular) are available as supplementary materials and as a Dryad repository, https://doi.org/10.5061/dryad.bg79c npbv.