Developing a monitoring and assessment program with diatoms, an improved metric calculation method, and causal analysis for Big Cypress National Preserve, Florida (USA)

Developing a monitoring and assessment program with diatoms, an improved metric calculation method, and causal analysis for Big Cypress National Preserve, Florida (USA) | Research Square window.SnipcartSettings = { analytics: { enabled: false } }; (function() { var accessVector = localStorage.getItem('access_vector') || ''; window.dataLayer = window.dataLayer || []; if (accessVector) { window.dataLayer.push({ user: { profile: { profileInfo: { snid: accessVector } } } }); } })(); (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0],j=d.createElement(s),dl=l!='dataLayer'?'&l='+l:'';j.async=true;j.src='https://www.googletagmanager.com/gtm.js?id='+i+dl;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-K279D39R'); Browse Preprints In Review Journals COVID-19 Preprints AJE Video Bytes Research Tools Research Promotion AJE Professional Editing AJE Rubriq About Preprint Platform In Review Editorial Policies Our Team Advisory Board Help Center Sign In Submit a Preprint Cite Share Download PDF Research Article Developing a monitoring and assessment program with diatoms, an improved metric calculation method, and causal analysis for Big Cypress National Preserve, Florida (USA) R. Jan Stevenson, Kevin R. T. Whelan, Michelle C. Prats This is a preprint; it has not been peer reviewed by a journal. https://doi.org/ 10.21203/rs.3.rs-6959033/v1 This work is licensed under a CC BY 4.0 License Status: Published Journal Publication published 15 Mar, 2026 Read the published version in Environmental Monitoring and Assessment → Version 1 posted 8 You are reading this latest preprint version Abstract We used relationships between diatom species composition and phosphorus concentration in periphyton (mat P) to develop a monitoring and assessment program for wetlands in the Big Cypress National Preserve (BCNP) in Florida, USA. Cluster analysis and regression showed limitation of taxonomic composition of assemblages to a few low P taxa in low mat P conditions, with additional species being able to invade habitats that have higher mat P concentrations. TITAN (Threshold Indicator Taxa ANalysis) and regression respectively identified 11 and 8 low P taxa and both methods identified 47 high P taxa. Congruence of our study results and results of phosphorus experiments confirmed P caused species responses. Metrics using the low and high P taxa traits were highly related to mat P. Metrics calculated with TITAN derived taxa traits and a novel log transformation of relative abundances were most highly related to mat P, differed most among regions with different human disturbance, and were unaffected by natural factors. Benchmarks for management targets were observed for minimally disturbed conditions and for greater than usual changes in assemblages along the P gradient. Diatom metrics were more highly correlated with distance from P sources than mat P, indicating species-based metrics have high value for monitoring and assessment. Our diatom metric development methods, benchmarks in ecological response and minimally disturbed condition, and causal analysis provide multiple new findings integrated for a monitoring and assessment program with effects-based management targets in BCNP. In addition, we advance important examples for applications in other ecological settings. algae ecology management metrics periphyton wetlands phosphorous Figures Figure 1 Figure 2 Figure 3 Figure 4 Figure 5 Figure 6 Figure 7 Figure 8 Figure 9 Introduction Building an assessment program depends on management goals and an understanding of the ecosystem, from which variables are selected to test their application in assessment and for establishing management targets (USEPA, 1992 ; Stevenson et al., 2004 ; Patterson et al. 2008 ). The basics of understanding the ecosystem inform which attributes are the target of management goals, what contaminants and habitat alterations (herein called stressors, sensu USEPA, 1992 ) affect those management goal attributes, and what human activities produce those contaminants and habitat alterations (Stevenson, 2011 ; Tang et al., 2020 ). Measurements of the key direct and indirect variables in the ecosystem that regulate management goal attributes should be used to characterize the condition of the ecosystem relative to management goals. Assessment involves characterizing condition. Assessment can also call for comparing condition to benchmarks in variables that can be used to determine if ecosystem condition meets management goals and whether management actions are warranted. Diagnosing which contaminants and habitat alterations are most affecting management goal attributes can also be parts of assessments, so remedial management actions can be selected (USEPA, 2000 ; Cormier et al., 2008). The goal of our research was to develop a periphyton monitoring program for Big Cypress National Preserve (BCNP). The overarching goal of the BCNP periphyton monitoring program is to support the United States National Park Service (NPS) mission to preserve the unimpaired condition of natural resources of National Parks (Patterson et al., 2008 ). Periphyton have been studied extensively in greater Everglades area, which is located southeast of the BCNP. In the greater Everglades area and BCNP, a distinctive calcareous periphyton is observed in areas with minimal human disturbance (Browder et al., 1994 ; La Hée et al., 2012). Calcareous periphyton have important ecological functions, respond directly and sensitively to phosphorus pollution, and effects of phosphorus on periphyton cascade through the rest of the ecosystem and its biota (Gleason et al., 1974; McCormick et al., 1997 ; McCormick et al., 2001 ; Gaiser et al., 2005 ; Gaiser et al., 2006 ; Hagerthey et al., 2011 ). For these reasons, periphyton were selected as one of the vital signs for BCNP (Patterson et al., 2008 ) with park-specific goals to determine status and trends in periphyton, water quality, and ecosystem function (Urgelles et al., 2019 ). Differences in soils, wetland transmissivity, and phosphorus range between the BCNP and the greater Everglades (Jarosewich & Wagner, 1985 ) were deemed sufficient to warrant development and testing of a periphyton program specifically for BCNP. Often, taxa traits for a new monitoring program are determined using past research in other regions. Lavoie et al. ( 2009 ) showed that taxa indicator traits from Van Dam et al. ( 1994 ), largely based on research and application in Europe, were related to pollution indicators in Canada. Tang et al. ( 2016 ) successfully used taxa traits developed for the western United States for assessments across the country. However, Potapova et al. (2002) found regional differences in indicator values across ecoregions in the U.S. Therefore, it seems reasonable to start developing a new program by using taxa trait characterizations from past projects in regions with similar types of human disturbance, pollution, and natural variability. Then, when sufficient data has been accumulated, as in BCNP, refine those taxa traits for the specific region that will be monitored. Many types of metrics have been proposed, ranging from diversity measurements to indicator species abundances and weighted average models that infer stressor conditions (Stevenson 2014 ). Given the goal of the BCNP program, metrics should be assessed for their ability to determine the deviation of the assessed ecosystem from minimally disturbed ecological condition. Widely accepted definitions of ‘minimally disturbed condition’ and the ‘attributes of biological condition’ have been proposed and serve as standards for many ecological assessments (Davies et al., 2006; Stoddard et al., 2006 ). In addition to measuring deviations in biological condition from a minimally disturbed condition, non-linear attributes of metrics can be used to develop consensus for management targets (Muradian, 2001 ). Tipping points in relationships showing assimilative capacity can provide justification for protecting ecosystems at pollution levels below those tipping points and thus preventing major changes in ecosystems (Soranno et al., 2008 ; Stevenson et al., 2008 ; Stevenson, 2011 ), such as threshold changes in calcareous periphyton abundance in the Everglades (Stevenson et al., 2002 ; Stevenson, 2014 ). In addition to diversity traits and weighted average models that infer specific pollution levels, measures of changes in proportion to sensitive native taxa and pollution tolerant taxa align with standard attributes of biological condition (Davies et al., 2006). Some studies show that a metric characterizing the proportion of the taxa that are sensitive or tolerant to pollution can perform as well as or better than metrics that include weights for species relative abundances in samples (Wang et al., 2005 ; Stevenson et al., 2008 ; Charles et al., 2021 ). The high performance of metrics that measure changes in the proportions of taxa with traits seems beneficial because they are more likely directly related to measuring losses in sensitive taxa or gains in tolerant taxa, even though there are problems with accuracy of these metrics measuring losses or gains in taxa when observing such small proportions of the organisms present. However, weighting relative abundances of taxa in samples should provide additional information that would make metrics perform better than simpler metrics that only account for presence and absence of taxa. Perhaps variability in relative abundance of abundant species is so great that it reduces precision of metrics that include relative abundances. If that is the case, down-weighting abundances with log-transformations could improve metric performance. The specific goals of this paper are to characterize the effect of human disturbance on biological condition of periphyton in the BCNP and to develop and test metrics for a periphyton monitoring program. First, we characterize changes in biological condition with ordination, cluster analysis, and metric changes along a human disturbance gradient. During metric evaluation we compare metrics with taxa phosphorus traits from the literature with revised traits from BCNP. We also test a novel metric calculated using log-transformed relative abundances to reduce negative effects that high variability in abundant taxa could have on metric precision. We also determine whether natural variability among habitats sampled affects metrics, and we characterized metric change along the phosphorus gradient to inform development of effects-based biological and pollution management targets. Using those analyses, we evaluate metric performance to aid selection of a subset of tested metrics for use in future monitoring. Finally, we explore relationships between metrics and the major pollutant, indicated by mat P, to provide benchmarks in the mat P gradient for potential management targets. Methods Study area Big Cypress National Preserve is located southwest of Lake Okeechobee and borders the northern edge of Everglades National Park (Fig. 1 ). It spans 295,016 ha and receives approximately 425,000 visitors annually who come to hike, canoe, camp, bird-watch, hunt, fish, and use off-road vehicle (ORV) trails (Patterson et al., 2008 ). The preserve contains a large remnant of natural wetland mosaic including cypress strands and domes, pine forests, wet prairies, marshes, sloughs, mangrove forests, and hardwood hammocks. The preserve also contains large stands of dwarf cypress, as well as rare orchids, bromeliads, and ferns (Patterson et al. 2008 ). Big Cypress is fundamentally different than its adjacent Everglades neighbor. The soil is comprised of sands (Pamlico Sands) overlying caprock (hardened limestone) of the Tamiami Formation whereas the vast majority of the greater Everglades system has an organic peat substrate overlying Miami Limestone (Jarosewich & Wagner, 1985 ). Water transmissivity is more localized and slower in Big Cypress than the Everglades (Jarosewich & Wagner, 1985 ). Concentration of surface water total phosphorus is typically higher in Big Cypress compared to the greater Everglades. This has been attributed to higher total phosphorus from natural sources such as shallow soils, rocks, and ground water (Miller et al., 2004 ). The name “Big Cypress” refers to the vast expanse of cypress ( Taxodium distichum var. imbricarium (Nuttall) Croom) rather than to the size of the trees. The larger bald cypress trees were logged during the past two centuries. Extraction of oil, gas, and minerals occurs within BCNP and surrounding areas to the north (Whelan et al., 2020 ). North of BCNP and south of Lake Okeechobee, agriculture fields have been operating since the early 1900s. To establish the farming industry, the lands were first drained, and the natural water flow was redirected. The western region of Big Cypress has headwaters that begin north of BCNP and are conveyed south via a canal system (Barron River Canal, Miller et al., 2004 ). Since BCNP establishment in 1974, the city of Immokalee’s population has grown 6 times greater (just north of the Preserve) (US Census data 1980, 2020). In addition, farming and other industrial development has expanded north of the BCNP boundary. The primary management activities have been related to restoration of the human altered regional hydrology to a more natural flow pattern, improving water quality entering the park from the north, manage invasive species, balancing recreational and extractive uses with long-term sustainability of the system, as well as protecting and preserving natural resources that include rare species (Urgelles et al., 2019 ). Sampling design Periphyton sampling has been conducted since 2008 to investigate potential and test methods for using periphyton to assess ecological condition in BCNP. The target sampling population was periphyton within mapped graminoid and broadleaf marshes in northwestern BCNP. The sampling area has been restricted to mapped and accessible marshes (which includes access by helicopter, ORV, truck, or hiking) in the northwest section of BCNP because that section has regions close to areas of human disturbance with elevated phosphorus concentrations as well as minimally disturbed conditions regions in the southern and eastern parts of the sampling area. Sample units were contiguous graminoid and broadleaf marsh habitats within a 250 m radius of selected map grid cell centroids from the Western Big Cypress National Preserve Vegetation Map (Whelan et al., 2020 ). The northwestern BCNP is divided into basins separated by artificial or natural structures that may or may not impede water flow during part or all the year (Fig. 1 ). Each basin is treated as a separate sampling strata or block. Using a restricted stratified random design, a number of potential sites are selected within the desirable habitat of each basin. Each potential site was evaluated and those that meet specific criteria were sampled annually (Londoño, 2019 ). The goal is to sample the same site year after year (Urgelles et al., 2019 ). For hydrological sample years 2013 and 2014, sampling followed very similar criteria as above except some sites were haphazardly selected to allow for: (a) maximizing spatial spread within a basin, (b) co-location with water-quality monitoring stations where available, (c) answering specific monitoring questions and (d) ensuring half the sites were accessible via helicopter and the other half of the sites accessible via ORV trails. Additionally, sampling occurred in the Kissimmee Billy Strand basin, which is an eastern basin that was only sampled in 2013 and 2014. Field collection methods Sampling goals and methods varied during early years to explore potential for periphyton monitoring. Starting in 2013, mat phosphorus (mat P) concentrations were measured as well as periphyton species composition. Mat P concentrations vary less than water column P as indicators of P pollution and have been recommended for measurement of P in the Everglades (Gaiser et al., 2004 ). Therefore, periphyton metric development was limited to hydrologic years after 2013. For hydroyears 2013, 2014, 2017, and 2019, a general overview of methods upon site arrival was to locate where water was present and collect a floating periphyton sample if possible. Attempts were made to access the same depressional marsh for revisits. Sampling followed a hydrological clock, approximately two months after the peak of the wet season, whereby sites were accessed at a similar hydroperiod stage year after year (typically November or December). Thus, water was usually present at all sites during collection. However, if the site was dry, then dry periphyton was collected. Separate grab samples were collected for diatom and mat P analyses, with each composed of a minimum of five grabs within a 5-meter radius sampling area to get two 125 mL samples. The preferred order of substrate collection was (1) floating mat, (2) algae on plants (sweaters), (3) algae on sediments (benthic mats), and (4) algae on woody debris. Grab samples are placed inside two 125-milliliter Nalgene ® opaque bottles with as little water as possible. Excess water was decanted prior to fixation. Additional site characterization information was collected including vegetation data, pH, water conductivity, and water depth. Samples were preserved as quickly as possible upon return from the field. Periphyton mat P samples were chilled in an ice slurry in the field and frozen at the lab. The periphyton composition samples were fixed with a 3% buffered formalin solution (Muxo & Shamblin, 2019 ). For hydroyears 2020, the only modification was the periphyton grab samples did not have a preferred substrate collection order but instead, periphyton was collected in relative proportion of the substrates present at the site, e.g. as floating mat, sweaters, benthic map or as debris. Otherwise, sampling was the same as above. Sample processing and analysis Periphyton samples were sent to Florida International University for the assay of mat P, which involved digesting mats and assaying total phosphorus (TP). In the laboratory, periphyton sample wet weights were recorded, thawed, and then extraneous plant material and animals were manually separated from the periphyton and placed in foil packets. These were placed in an 80°C oven for three days; once dry, the weight was recorded as “extraneous material”, and then discarded. The TP content of periphyton was expressed on a dry-weight basis because the organically incorporated P was not separable from that bound to calcite. The TP subsample was dried at 80°C and ground down to a fine powder with a mortar and pestle. Colorimetric analysis was used to estimate TP concentrations of the periphyton subsamples following the methods of Solorzano and Sharp ( 1980 ), which were then used to calculate µg g − 1 dried periphyton mat (not ash free dry mass). For hydroyears 2013, 2014, and 2017 the lab methods followed Smith ( 2018 ) and for hydrologic year 2019 and 2020 they followed Wilson ( 2020 ). The difference in methods is due to the acquisition of a newer spectrophotometer. Early, unpublished analyses of data indicated diatom species composition was related to the human disturbance gradient in BCNP as well as composition of both diatom and non-diatom algae. Afterward algal analyses were limited to diatom species composition for development of diatom metrics to reduce costs of periphyton analyses. Periphyton samples were sent to Michigan State University to be assayed for diatom species composition. There, samples were acid cleaned in nitric acid and mounted on microscope slides in NAPHRAX©. At least 600 diatom valves were identified and counted for each slide. Data analysis: data selection and descriptive analyses We created independent calibration and validation datasets with only one sample per site in each dataset to avoid problems with pseudoreplication and repeated measures. Both subsets of samples were created using samples from the 2013–2020 period when mat P was assayed and by randomly selecting one sample from each site without replacement of selected sites. This created a calibration dataset with 113 samples for characterizing changes in diatom assemblages along the P gradient and for characterizing species traits. 78 of the 113 calibration sites had all chemistry information. 105 sites had all chemistries except pH. A second set of samples from the 2013–2020 period with chemistries was selected for a validation dataset to test metrics, which resulted in 85 randomly selected samples with one sample per site that was not selected for the calibration dataset. Ordination analyses indicated that environmental factors, sampling year, and taxonomist were related to variability in diatom species composition. Follow-up analyses of the data indicated that residuals in diatom metric relationships with mat P concentration also differed among sampling years. Different taxonomists were assigned to count diatom samples from year to year with turnover of taxonomists present in the lab. Interannual residual variation could not be assigned to either taxonomist or environment with certainty, so we restricted metric calculation and testing to hydrologic years 2017, 2019, and 2020 when the same taxonomist counted all samples. That reduced the validation dataset for testing metrics to the 44 sites sampled annually in the 7 regions with one sample per site for either 2017, 2019, or 2020. However, we used the full calibration dataset from 2013–2020 for the preliminary descriptive analyses: ordination, cluster analysis, trait characterization, and characterization of minimally disturbed condition. Our rationale was despite the concern about potential effect of taxonomist or interannual environmental differences, we wanted to ensure sufficient sample size to detect species-environment relationships and characterize traits for as many taxa as possible so future results would be resilient to interannual variation caused either by environment or taxonomist. Species-environment relationships were evaluated with ordination using the calibration dataset. Non-metric multidimensional scaling (NMDS) was selected to ordinate samples in species space and then relate environmental variables to the NMDS axes with the vegan package in R (Oksanen et al., 2020 ). Cluster analyses using Bray-Curtis dissimilarity and the calibration dataset were calculated with the clust command in R to observe patterns in species composition without constraints of environmental factors. Groups of sites with low dissimilarity were identified and then differences in mat P, conductivity, and pH among these site groups were determined. To more thoroughly understand and illustrate the succession in diatom assemblages along the P gradient, we used a stacked bar chart using the most abundant low and high P indicator species. Relative rates of change in species relative abundances and the dominance of these taxa in assemblages along the P gradient were determined for the calibration dataset. The most abundant low P and high P taxa were selected for the figure. In addition, we created a heat map of all taxa in 5 samples or more using an Excel© spreadsheet. The heat map matrix had taxa abundances colored in rows and samples in columns, with rows ordered by TWINSPAN groups (R package) and columns ordered by increasing mat P concentrations. Cells of the heat map matrix were colored yellow, orange, or red to indicate increasing relative abundances of taxa. Data analysis: taxa trait characterization Taxa have been characterized by the NPS as being characteristic of oligotrophic or eutrophic habitats based on consensus and best professional interpretation of information reported in literature sources (Table S1 ). This literature review assessed 13 papers and attempted to determine the optimum TP, trophic habitat preference, pH preference, and trophic indication status of diatoms at the species level. Traits analysis was conducted to characterize taxa as low and high P taxa using the calibration dataset. Multiple trait characterization methods were used to be thorough and to evaluate consistency in characterizations. Linear and polynomial regression were used to determine whether relationships between taxa relative abundances and mat P concentration were negative or positive to identify low and high P taxa, respectively. Polynomial regression was used to complement linear regression to ensure non-linear responses were detected. TITAN (Threshold Indicator Taxa Analysis (Baker & King 2010 )) was also used to describe taxa with decreasing and increasing indicator species values along the TP gradient for low and high P taxa, respectively. Finally, weighted average TP optima were calculated for taxa. Low and high P characterizations (i.e. taxa traits) were only calculated for taxa observed in 5 or more samples from the calibration set of 78 samples from different sites and having P concentrations measured because TITAN requires 5 observations per taxon to provide reasonable statistical power for characterizing taxa traits with regression. Mat P optima and tolerances were determined for all taxa because our experience showed weighted average models (WAM) were most precise when all taxa are used in the model. The NPS does not plan to use the WAM model in future monitoring. The WAM model is used in this paper to provide a benchmark for a model that is usually the most precisely related to an environmental gradient (Reavie et al., 2008 ). Taxa traits determined with the 78-sample calibration dataset were compared to taxa traits observed in past Everglades research to compare consistency in results and evaluate causal relationships between low and high P traits and P concentrations. Pairwise relationships between taxa trait characterizations in this study and those by Hagerthey et al. ( 2012 ) and Gaiser et al. ( 2006 ) were determined. Species characterizations for BCNP were also compared to low and high P characterizations by Slate and Stevenson ( 2007 ), which were determined by experimental manipulation of P with in situ mesocosms. The experiments enabled establishing a causal relationship between P concentration and relative abundances of diatom taxa (Beyers 1998 ). Data analysis: metric evaluation We then evaluated metric design on their performance to be sensitive to environmental change, to characterize changes in biological condition along the human disturbance gradient, and to help justify management targets for P concentration. We did not expect the same metric to be best for all goals (Fig. S1 , Stevenson, 2006 , Stevenson, 2011 ). We sought a metric with a linear response to P concentration to be sensitive to environmental change, because metrics with linear response have the same sensitivity to incremental changes in P concentration at all levels of P concentration. We sought non-linear responses of metrics of biological condition to detect thresholds in biological responses that could be used to justify P management targets. We did not calculate a multimetric index because that was not the plan for BCNP application. For each sample we calculated multiple metrics with each trait, plus three diversity metrics: taxa number observed in standardized counts, Shannon diversity, and Pielou’s evenness. Counts were standardized to 600 valves, because some counts had more than 600 valves. Counts were standardized by assigning random numbers to valves in counts and picking valves with the lowest 600 numbers. For each trait, four metric types were calculated: 1) proportion of individuals with a trait (PropValves) as sum(v ijt /V j ), where v ijt is the number of valves counted for taxon i with trait t (e.g. either high or low P using a specific trait determination method) in sample j , and V j is the number of valves of all taxa with assigned traits in sample j ; 2) the proportion of taxa with a trait (PropTaxa), as t tj /T j where t tj is the number of taxa with trait t in sample j and T j is the number of taxa with traits assigned in sample j ; and 3) the number of taxa in a sample with a specific trait (noTaxa, e.g. number of low P taxa). The fourth metric type was RlogA and is described in the next paragraph. Below we also describe our use of the three trait determination methods: review of the literature, TITAN, and regression. Recent papers have shown metrics calculated as proportion of taxa with a trait are often related better to human disturbance and stressor gradients than metrics calculated as proportion of individuals with a trait (Stevenson et al., 2008 b; Carlisle et al., 2022 ). Because relative performance of species in a habitat likely varies with environmental conditions as well as their presence in a habitat, we tried down-weighting metrics by log transforming relative abundances of taxa and thereby reducing variability caused by slight changes in growth rates of abundant taxa that could produce great differences in a taxon’s abundance. We determined relative log abundance (RlogA ijt ) by calculating the natural log of valve numbers of taxon i with trait t in sample j and then the sum of natural log transformed valve numbers for all taxa with traits assigned in sample j as sum(RlogA ijt=h+l ), summed for all taxa with both high or low P value (t = h + l) assigned with a specific method (literature, TITAN, regression). The RlogA metric values for either high or low P RlogA were calculated as the sum of all RlogA ijt for taxa with either low or high P traits based on a specific trait determination divided by the sum of all RlogA ijt for taxa with both low and high P traits, i.e. RlogA j = sum(RlogA ijt )/SRlogA j . Metrics were evaluated by comparing their adjusted R 2 values determined for slopes in relationships between metrics and log transformed mat P. We log 2 transformed mat P to even the distribution of the observations along the mat P gradient, because the number of low mat P samples was much greater than the number of high mat P samples. Unless described otherwise, a metric-mat P relationship in the following text refers to the relationship between a metric and log transformed mat P. The R package lm was used for linear regression (R Core Team, 2021 ). Both actual values and ranks of R 2 values were evaluated. Three analyses of variance (ANOVA) were used to compare average adjusted R 2 for metric-mat P relationships for: 1) different metric types (noTaxa, PropTaxa, PropValves, RlogA); 2) trait source (literature, regression, or TITAN);, or 3) trait (low P versus high P). To control for other metric attributes (metric type, trait source, and trait) when evaluating performance of one of the metric attributes, Kruskal tests ( kruskal.test in R) were run to compare ranks of adjusted R 2 of metric-mat P relationships. Because these analyses are rather complicated to set up, we describe them in detail. Three Kruskal tests were run to evaluate metric attributes, and each used the same set of 24 metric-mat P relationships resulting from the 24 combinations of 4 metric types, 3 trait sources, and 2 traits. To compare performance of metric types (noTaxa, PropTaxa, PropValves, RlogA) while controlling for trait source and trait, we grouped metric-mat P relationships into the 6 different trait source-trait groups (low P-Lit, low P-Regr, low P-TITAN, high P-Lit, high P-Regr, high P-TITAN). Then we ranked (1–4) the adjusted R 2 of the metric-mat P relationships for the 4 metric types within each of the 6 trait source-trait (TS-T) groups. In the case for the Kruskal test comparing metric types, ranks ranged from 1–4 because there were 4 metric-mat P relationships having different metric types within each TS-T group. In the Kruskal test for metric type, each metric type had 6 rankings because there were 6 TS-T groups. Using the same approach we compared performance of metrics with either Lit, Regr, or TITAN trait sources while holding metric type and trait constant; we used the Kruskal test to compare ranks (1–3 trait sources) of adjusted R 2 for the 3 metric-mat P relationships in the 8 metric type-trait source (MT-TS) groups (noTaxa-Lit, PropT-Lit, PropV-Lit, RLogV-Lit, noTaxa-Regr, PropT-Regr, PropV-Regr, RLogV-Regr, noTaxa-TITAN, PropT-TITAN, PropV-TITAN, RLogV-TITAN). To compare performance of metrics with either low P or high P traits while holding metric type and trait source constant, we used the Kruskal test to compare ranks (1 or 2) of adjusted R 2 for those 2 metric-mat P relationships in the 12 trait source– metric type groups (TS-MT groups: Lit-noTaxa, Lit-PropT, Lit-PropV, Lit-RLogV, Regr-noTaxa, Regr-PropT, Regr-PropV, Regr-RLogV, TITAN-noTaxa, TITAN-PropT, TITAN-PropV, TITAN-RLogV). Mat P and metrics were related to latitude (lat) and longitude (long) to determine their relationship to distance from the main source of human disturbance in the northwestern corner of the BCNP study region. Linear regression ( lm in R) was used to relate mat P and metrics to latitude, longitude, and a lat-long interaction term. Differences in metrics among regions were determined with the 2020 data to evaluate their performance and demonstrate application by determining regional differences that are related to levels of human disturbance. The regional approach was chosen because the NPS plans to analyze and report results by region using sites as replicates as well as to determine differences among regions and changes in time by region. Single factor ANOVA were used to compare average metric values among regions. p values resulting from these analyses were used to evaluate metric performance with the average p rank for metric types (noTaxa, PropTaxa, PropValves, and RlogA), low and high P traits, and trait sources (literature, regression, and TITAN). Non-parametric Kruskal tests were again used, as used to compare metric-mat P relationships above, to determine which metric attributes (metric type, trait source, trait) were most important for metric performance when comparing regions. This set of three separate Kruskal tests for metric type, trait source, and trait used ranks of 24 ANOVA F-values for comparing differences in the 24 possible metrics among regions. Kruskal tests used ranks of F-value ranks ranging from 1–4 to compare metric types within the 6 trait source-trait groups, ranging from 1–3 to compare trait sources within the 8 metric type-trait groups, and ranging from 1–2 to compare the two traits within the 12 trait source-metric type groups. In addition, we calculated an ANOVA for mat P differences among regions to compare F-values with diatom metrics to determine if diatom metrics more precisely differed among regions than mat P. Data analysis: metric correction for natural factors Residuals in metric-mat P relationships were related to naturally varying ecological factors to determine if expectations for metrics should be adjusted for habitat type, substrate location, substrate type, and mean water depth. Habitat type included broadleaf marsh, cypress dome, cypress scrub, graminoid marsh and mixed combinations of these habitats. Substrate locations were either floating on the water surface, on benthic sediment or soils (ground), or enveloping macrophytes (sweaters). Substrate types were filamentous algae, periphyton, and soil. Water depth was measured at the location and could covary with the major human disturbance gradient originating with human alterations in the northeast corner of the sampling area. Univariate ANOVA were used to determine separately the effects of differing habitat types, substrate locations, and substrate types on metric-mat P relationships. Linear regression in R ( lm ) was used to characterize the relationship between water depth and residuals of the metric-mat P relationship. p values for statistical significance were reported in results without accounting for multiple tests, but multiple tests were accounted for by dividing reported p values by the number of tests performed when interpreting likelihood that observed results were not due to random patterns in the data. This is like a Bonferroni correction for multiple tests. In addition, we looked at the proportion of p values that were less than 0.05. Residuals in the relationships for metrics and mat P as function of latitude and longitude were related to natural environmental factors to determine if natural environmental factors affected both mat P and metrics. This second residual analysis differed from the first analysis, because it addressed whether different habitat conditions had different background P conditions and metric values versus the first analysis which evaluated whether natural environmental factors introduced bias in metric-mat P relationships. Data analysis: management benchmarks We used multiple lines of evidence to establish benchmarks for assessment and possible management targets for mat P and periphyton metrics. Here we use the term benchmark to mean a level of an environmental variable that could be used by resource managers to set management targets. To aid establishment of management targets for mat P, we evaluated non-linearities in metric-mat P relationships for tipping points (changes in the linear patterns) at the low-P end of a mat P range, thereby showing natural assimilative capacity for the biological condition attributes measured by a metric (Fig. S1 ). We used three approaches for determining benchmarks in biological response along the mat P gradient. Decision tree analysis was calculated relating the 29 trait-based, diversity, and WAM metrics to mat P using the rpart package in R (Therneau & Atkinson, 2025 ) to determine change points along the mat P gradient in trees with the lowest error variance. We used the 75th percentile of mat P concentrations for low P groups of samples identified with cluster analysis as another benchmark. We also used the point with highest TITAN z scores for low P taxa as a management benchmark (Baker et al. 2010). We compared benchmarks to changes in species composition illustrated in the heat map of taxon abundances in samples with successively higher mat P to interpret biological responses associated with these mat P benchmarks. We related mat P to latitude and longitude to determine if it was related to distance from the northwest corner of the BCNP study area where human disturbance was highest. To ensure management targets based on effects of mat P did not under or over-protect BCNP wetlands, we determined likely mat P concentrations at sites we assumed were minimally disturbed sites. We selected all sites in the Fire Prairie (FP), Monument (MN), and Little Marsh (LM) regions because they were the regions furthest south and east in the BCNP study area. However, some sites in these regions were close to canals and trails which could introduce disturbance. Results Changes in diatom assemblages along the P gradient NMDS with all samples with all chemistries in the calibration dataset (n = 78) showed diatom assemblages responded most to mat P. The NMDS had an insignificant deviation from 1 in the stress plot. The continuous variables mat P, year sampled, pH, water temperature, and conductivity as well as the categorical variables habitat, water color, and taxonomist were significantly related to NMDS axes ( p ≤ 0.041, Fig. S2). Mat P, sampling year, and habitat were the most highly correlated abiotic variables with 0.624, 0.349, and 0.273 R 2 values, respectively ( p < 0.001, Table S2). pH and water temperature were next most closely related to NMDS axes, with 0.162 and 0.128 R 2 values respectively. Mat P and pH were inversely correlated to each other (Pearson r =-0.357) along the NMDS axes. Cluster analysis of all 113 samples at unique sites in the calibration dataset isolated 5 groups of sites at a dissimilarity level of 0.75 (Fig. S3). From low to high along the mat P scale, these groups were designated 1 to 5. Group 1 had more than two thirds of the sites. Groups 2 to 5 had 5, 1, 9 and 10 sites respectively. Group 3 was dropped from the following discussion because it only had one sample. Three subgroups in Group 1 that had dissimilarity less than 0.6 among sites were identified and designated as Groups 1–1, 1–2, and 1–3 in order of dissimilarity among groups. Group 1–1 had very low dissimilarity among sites and Group 1–3 had the highest dissimilarity of the Group 1 subgroups. Dissimilarity increased significantly (Tukey HSD, adj p < 0.001, Fig. 2 ) from a mean of 0.39 among samples in Group 1–1 to 0.46 and 0.47 in Groups 1–2 and 1–3, which were themselves not statistically different. Dissimilarity also increased significantly (Tukey HSD, adj p < 0.001) from Groups 1–2 and 1–3 to Groups 2, 4, and 5, which had 0.54 average dissimilarity for those three groups of sites. Mat P, pH, and conductivity, varied significantly among site groups (ANOVA, p < 0.05, Fig. 2 ). With 75th and 90th percentiles of 142 and 201 µg P g − 1 dry mass, mat P was lower in Group 1–1 than all other site groups (Tukey HSD, adj p < 0.001). Average mat P for sites in either Group 1–2 or 1–3 was lower than for sites in either Group 2, 4, or 5. The 75th and 90th percentiles of all sites in Group 1 were 354 and 508 µg g − 1 mat P. pH of Group 1–1 sites was greater than sites in either Group 2, 4, or 5 (Tukey HSD, adj p = 0.03–0.6). Average conductivity for sites in Group 4 was less than in Groups 1–3 and 5 (Tukey HSD, adj p = 0.04). Tukey HSD comparisons of mean chemistry values for pairwise comparisons of other site groups were not significantly different. The number of diatom taxa observed in 600 valve counts increased ( p < 0.001) from about 10 in low P to 25 in high P as high P taxa invaded and low P taxa became very rare or were lost (Figs. 3 , 4 , and S4). Encyonema evergladianum , Mastogloia calcarea, Encyonopsis microcephala , and Brachysira ocalanensis , assigned low P traits in calculations described later in these results, comprised approximately 65 percent of assemblages in low mat P conditions. Their combined proportion of assemblages decreased to 45, 20 and 5 percent in samples with mat P ranging, respectively, from 350–450, 450–750, and greater than 750 µg g − 1 mat (Fig. 4 ). The five dominant high P taxa (Encyonema silesiacum, Navicula cryptotenella, Nitzschia amphibia, Gomphonema gracilis, and Gomphonema auritum ) increased from 6 percent when mat P was less than 100 to about 45 percent of samples when mat P was greater than 750 µg g − 1 (Fig. 4 ), leaving a high diversity of other diatom taxa to comprise the rest of assemblages in high P conditions (Fig. S4). Whereas other common low P taxa maintained relatively high proportions of counts at intermediate concentrations of mat P, proportions of E. evergladianum decreased more rapidly than other low P taxa (Fig. 4 ). Encyonema evergladianum was observed in 56 of 58 samples with mat P less than 300 µg g − 1 and were usually quite common, however it was only observed in 4 of 28 samples in samples with mat P greater than 550 µg g − 1 (Fig. S4). In contrast, N. amphibia was only observed in 6 of 49 samples with mat P less than 200 µg g − 1 , but it was highly abundant in 14 of 16 samples with mat P greater than 800 µg g − 1 . Taxa traits Low and high P traits were assigned to 58 taxa with TITAN and to 55 taxa with regression (Table S3). Weighted average mat P optima were assigned to 157 taxa because the criterion of observing taxa in 5 or more samples was not applied for determination of mat P optima. 62 taxa had low and high P traits assigned from the literature by the NPS. Many more taxa were assigned high P than low P traits. The number of taxa assigned low and high P traits were, respectively, 11 and 47 taxa by TITAN, 8 and 47 taxa by regression, and 23 and 39 by the literature review. TITAN and regression did not assign different traits, either low P or high P, to any taxon; but some of the less common taxa were assigned traits by only TITAN or regression. Low and high P trait assignments did differ for one taxon, E. silesiacum , which was assigned a high P trait by TITAN and regression and a low P trait using literature references. Many of the abundant low and high P taxa determined by TITAN and regression did not have trait assignments from the literature. The averages of P optima were, respectively: 179 and 1001 for low and high P TITAN taxa; 187 and 1220 for low and high P regression taxa; and 394 and 1062 for low and high P literature-defined taxa. Encyonopsis evergladianum, M. calcarea, E. microcephala , and B. ocalanensis were the four most common low P taxa (Table S3, Fig. 4 ). The other four taxa classified as low P by TITAN as well as regression were an unknown species of Nitzschia (sp. 1 ), Fragilaria synegrotesca, Adlafia bryophila , and Nitzschia serpentiraphe . The most common high P taxa were E. silesiacum, N. cryptotenella, N. amphibia, G. auritum , and G. gracile (Table S3). These most common taxa were not necessarily the most limited to high and low P conditions. The relative preference of taxa for low and high P conditions, versus the absolute preference (low or high) was important for metrics only when calculating the WAM diatom inferred mat P. Of the abundant low P taxa, Encyonema evergladianum had the lowest mat P optimum (112 µg g − 1 , Table S3). Nitzschia amphibia had one of the highest mat P optima of the common high P taxa (1198 µg g − 1 ). Encyonema silesiacum and N. cryptotenella had relatively low P optima for taxa characterized by TITAN and regression as high P taxa (439 and 358 µg g − 1 respectively). Traits assigned in this project with BCNP samples were highly correlated with traits assigned to taxa in three other studies (Fig. S5). Taxalists did vary between our BCNP list and lists in Hagerthey et al. ( 2012 ), Gaiser et al. ( 2006 ), and Slate and Stevenson ( 2007 ). We felt confident comparing BCNP traits for 35 taxa with names that were clearly the same with Hagerthey et al., 23 taxa in Gaiser et al., and 33 taxa in Slate and Stevenson. Spearman correlations among trait values in paired taxalists were positive, all highly significant ( p < 0.001), and had R 2 0.676, 0.737, and 0.783 for Hagerthey et al., Gaiser et al., and Slate and Stevenson, respectively. Metric testing and evaluation: metric-mat P relationships Almost all metrics except the diversity metrics (Shannon H and Pielou’s evenness) were related to mat P with high statistical significance and adjusted R 2 (Table S4). When comparing metric performance, metrics calculated with the relative log abundances (RlogA) were more precisely related (based on adjusted R 2 ) to the mat P gradient than all other metrics, including WAM diatom inferred TP (Figs. 5 and 6 , Table S4). According to ANOVA and Tukey HSD, average adjusted R 2 for relationships between RlogA metrics and log 2 (mat P) were significantly higher than the average adjusted R 2 for mat P relationships with noTaxa metrics (Fig. 6 A, Tukey HSD, p = 0.016), but other pairwise comparisons among metric types did not differ (ANOVA, n = 24, p = 0.017; Tukey HSD, 0.09 < p < 0.83, Fig. 6 A). However, the ranks of RlogA metric performances were greater than noTaxa, PropTaxa, and PropValves metrics (Kruskal test, n = 24, p = 0.001; Fig. 6 D; Table S4). Ranks of adjusted R 2 values for relationships between RlogA metrics and log 2 (mat P) were higher than all other 18 metric relationships with mat P (Fig. 6 D). Using the same ANOVA tests of adjusted R 2 and Kruskal tests for ranks of metric-mat P relationships, TITAN-based metrics performed better than metrics based on literature or regression traits. Averages for adjusted R 2 for metric-mat P relationships, as well as their ranks, were lower for Lit than TITAN and regression sources of traits (Figs. 6 B & 6 E, Tukey HSD p < 0.026, Kruskal test p < 0.0001). Although there was great overlap in adjusted R 2 values of metric-mat P relationships based on traits determined by TITAN and regression (Fig. 6 B), for 7 out of 8 comparisons TITAN traits had higher ranks than regression traits for adjusted R 2 ( p = 0.008, Fig. 6 E). Neither ANOVA or Kruskal tests indicated differences in performance of low and high P metrics (Figs. 6 C, 6 F). As with metric-mat P relationships, all metrics except the diversity metrics Shannon H and Pielou’s evenness were related with high statistical significance and adjusted R 2 to our geographic indicators of human disturbance, latitude, longitude, and the lat-long interaction term (Table S5). In a comparison of which metrics and mat P were most precisely related to the lat-long human disturbance gradient, all metrics using regression and TITAN traits and both WAM metrics had greater adjusted R 2 than mat P. Mat P was related to the lat-long model with a 0.44 adjusted R 2 , whereas the four metrics with highest adjusted R 2 , which included both WAM metrics, had adjusted R 2 ranging from 0.57 to 0.65 (Table S5). The metric with the high 0.65 R 2 for a lat-long model was RlogA high P PropValves. Metric testing and evaluation: metric differences among regions All metrics differed among regions with high statistical significance in 2020, with the highest P effects in the OK region and lowest in the MN region (Table S6, Figs. 7 and S6). In OK, metric values were highest for diatom inferred mat P, Shannon’s H, proportion of taxa with high P traits and relative log-transformed abundance of high P taxa valves. In OK, metric values were lowest compared to all other regions for proportion of taxa with low P traits and log-transformed valve abundances with low P traits. Overall, OK and then EH had the highest indicators of P pollution based on relatively high richness, WAM mat P metrics, and high-P trait metrics as well as relatively low low-P metrics compared to other regions (Figs. 7 and S6). MN and LM had the consistently lowest richness, WAM mat P, and high P metrics as well as the highest low P metrics. BI, EC, and FP had intermediate levels of disturbance related to the other two groups of regions (Figs. 7 and S6). Metric performance for detecting differences among regions was greater for the 2 diatom inferred mat P metrics (WAM_matP, WAM_RlogA.matP) compared to all other metrics as indicated by ANOVA F values (Table S6). A Kruskal test indicated ranks of metric types differed ( p = 0.025), with noTaxa and RlogA metrics having more highly ranked F values than PropTaxa and PropValves metrics (Fig. S7). Metric performance for distinguishing regions was better for TITAN and regression traits than literature traits (Fig. S7, Kruskal test, p = 0.009). There was no difference in low and high P metrics ability to detect differences among regions (Fig. S7). Mat P was significantly different among regions (ANOVA, p < 0.0001). Metric testing and evaluation: corrections for natural factors The following results describing the relationship between residuals from metric-mat P relationships and these habitat parameters should be considered with caution. Unequal sample sizes among the natural habitat features presented challenges for some residual analyses. Whereas substrate locations had from 11–21 sites for either floating, ground, or macrophyte sweaters, of the 6 habitat types there were 33 graminoid marshes and from 2–3 of all other habitat types. For substrate types, there was only one soil sample when filamentous algae and periphyton had bigger sample sizes, 7 and 36 samples, respectively. Analyzing differences between the two substrate types, filamentous algae and periphyton, was an issue because shifts from periphyton to filamentous algae are associated with increasing P concentrations. Similarly, in this dataset water depth can be associated with human disturbance. Little evidence indicated that natural factors affected metric-mat P relationships (Tables S7, S8). Water depth was not significantly ( p < 0.05) related to residuals in any metric-mat P relationship (Table S7). Time of year sampled (indicated by the variable named season) had a significant ( p < 0.05) negative effect on 6 of 29 metrics, which were metrics for the numbers of high P taxa and diversity (Table S7). For substrate type, location, and hydrologic year, p values for statistical significance were seldom less than 0.05 (1–4 of 29 metrics, respectively) and not less than 0.05 if accounting for the multiple statistical tests conducted (Table S8). For habitat type, 5 of 29 metrics had a p value less than 0.05, but low sample size for habitats other than graminoid marshes limited certainty in this observation. Evidence was weak for natural factors affecting metrics after effects of human disturbance (lat-long model) were accounted for. Residuals from relationships relating either metrics or mat P with the lat-long model were so seldom related to natural factors with ANOVA or regression that the few tests with p < 0.05 could have occurred by chance. Attained significance ( p ) was not greater than 0.05 for any of the mat P or metric residuals relationships with water depth, time of season when sampling occurred, or substrate type. p was < 0.05 for 1, 5, and 8 of the residual relationships with habitat type, substrate location, and hydrologic year, respectively (Table S8). Thus, substrate location and hydrologic year were the most likely natural factors affecting mat P and metrics. Interestingly, lat-long related residuals for both WAM metrics and mat P were related to hydrologic year, but residuals for regression or TITAN metrics were not related to hydrologic year (Table S8). Management benchmarks Given that metrics were evaluated along a log-transformed mat P gradient to even distribution of observations along the mat P gradient for regression analysis, that affected linearity of relationships to mat P. Metrics were linearly related to log transformed mat P according to comparisons of linear and non-linear models (e.g. Figures 5 , S8). However, most metrics were non-linearly related to the mat P gradient when mat P was not log transformed (Fig. S9). With untransformed mat P, metrics with low P taxa traits had a negative exponential relationship with mat P with rapid decreases in low ranges of mat P and little decrease in high ranges of mat P. In contrast, metrics with high P traits increased asymptotically with untransformed mat P with rapid metric increases in the low mat P range and little change in metrics in the high range of mat P. Benchmarks for relatively sudden changes in taxonomic composition were observed along the mat P gradient for 4 of 13 metrics evaluated. We evaluated mat P relationships with the 3 diversity metrics, 2 WAM metrics, and 8 TITAN-trait metrics. We excluded trait-based metrics that used literature and regression traits from the following analysis to prevent redundancy with TITAN-trait metrics that were best related to log mat P. Changepoint analysis showed that a breakpoint for low P TITAN taxa occurred at 427 ug g − 1 mat P, above which all low P taxa numbers were low (Fig. 8 , Table S9). Changepoint analysis showed numerous statistically significant ( p ≤ 0.01) breakpoints occurring at the same mat P concentrations for both the number of all taxa and high P TITAN taxa: 71, 137, 464, 627, and 1246 µg g − 1 mat P (Table S9, Figs. 3 and 8 B). WAM-mat P had a changepoint at 564 µg g − 1 mat P, above which almost all values were high (Fig. S8). Changepoints along the mat P gradient were not observed for any other diversity, WAM, or TITAN metrics (Table S9). However, review of graphs of metrics showed a strong changepoint ≤ 500 µg g − 1 mat P for proportion of low and high P literature-defined taxa that was particularly evident when plotted along a log mat P scale (Fig. S8d). From the heat map (Fig. S4), changes in assemblages around 137 µg g − 1 mat P appeared to be associated with the shift from the most likely low P co-dominants being E. evergladianum, N. serpentiraphe , and Nitzschia sp. 1 to E. microcephala and B. ocalanensis . Invasion of some high P taxa to accompany low P taxa in samples was more frequent above 137 µg g − 1 . Above 464 and especially above 627 µg g − 1 , high P taxa typically were most abundant in samples with loss of the sensitive common low P taxa E. evergladianum, Nitzschia sp. 1 , and M. calcarea . TITAN analysis supported the relatively high certainty of change in assemblages from low to high P taxa in the 300–600 µg g − 1 range of mat P (Figs. 9 and S10). The peak in sum of low P taxa z scores was approximately 250 µg g − 1 mat P. The spread of the interquartile range of changepoints for low P taxa, i.e. from the 0.25–0.75 cumulative frequency quartiles, was from 300–400 µg g − 1 mat P and from 400–600 µg g − 1 for high P taxa. The overlap in those changepoint ranges was evident in the heat map of species relative abundances in samples along the mat P gradient (Fig. S4). Between 400 and 500 µg g − 1 mat P, as mat P increased, the frequency of high abundances of most low P taxa decreased greatly and yielded to an increase in numbers and abundances of high P taxa in samples (Fig. S4). However, some low P taxa had taxon-specific changepoints in indicator species values below the 300 µg g − 1 mat P (Fig. S10). In contrast, some low P taxa were able to persist in relatively high mat P conditions, such as E. microcephala and B. ocalanensis that had likely indicator value changepoints greater than 500 µg g − 1 mat P (Fig. S10). Encyonema silesiacum and N. cryptotenella had the opposite characteristics; they were high P taxa that were commonly in high abundance in low P conditions as quantified by their median indicator species thresholds near 350–400 µg g − 1 mat P (Fig. S10). Mat P concentrations in minimally disturbed regions of the sampling area can be used to characterize the natural variability in minimally disturbed conditions. Minimum and maximum mat P were 52 and 606 µg g − 1 in the calibration data set at the 32 sites sampled in the FP, LM, and MN regions. The median and 75th percentile of mat P values in FP, LM, and MN were 140 and 367 µg g − 1 mat P. Sixty-four percent of mat P values at the 32 FP, LM, and MN sites were less than 200 µg g − 1 , and 90% of mat P values were less than 500 µg g − 1 . Discussion Changes in diatom assemblages along the P gradient Multiple lines of evidence indicated that phosphorus concentration was the determinant of changes in diatom species composition in regions of the BCNP near areas of human disturbance, as in regions of the greater Everglades. Mat P and metrics were highly related to distance from human disturbance. NMDS analyses showed mat P was the variable most strongly related to variation in diatom species composition. The mat P optima for taxa determined with the BCNP calibration dataset correlated well with P optima identified in experimental manipulations of P in the Everglades Natural Park and Everglades Water Conservation Area 2A (Pan et al., 2000 ; Gaiser et al., 2006 ; Slate & Stevenson, 2007 ). Phosphorus contamination in different regions of the Everglades, typically from canal water entering marshes, is widely recognized as a major threat to biological condition in the Everglades (McCormick et al., 1996 ; Gaiser et al., 2006 ) and many other wetland ecosystems as well (Pan et al., 1996; Lougheed et al., 2007 ; Wyatt et al., 2010 ; Pillsbury et al., 2019 ). The number of taxa observed in BCNP samples (alpha diversity) increased with mat P concentration from an assemblage restricted to very few low P taxa to a wider diversity of taxa requiring high P concentrations, as with other regions of southern Florida (Raschke, 1993 ; Pan et al., 2000 ; Gaiser et al., 2006 ). Given the widespread sampling in the greater Everglades region and low number of taxa observed in low P conditions, it seems likely that diversity of diatoms is indeed low, despite limitations of 600 valve counts. Low diversity in naturally low nutrient conditions could be due to few taxa are adapted to the stress of low nutrient supply (Worm et al., 2002 ; Cardinale et al., 2009 ). Then, as nutrient concentrations increase, habitat availability increases for a larger number of taxa that require higher P to persist. Some low P taxa were observed occasionally in moderate and even high P sites, however other low P taxa were not observed in high P sites. Although showing that taxa have indeed been extirpated from a habitat is difficult, especially for microbes for which we observe such minute proportions of populations, the lack of occurrence of taxa such as E. evergladianum in the large number of counts with sites having more than 550 µg g − 1 mat P indicates they are very rare in high P areas even though they are the most common diatom in minimally disturbed condition. Diversity of taxa among sites (beta diversity) also increased with P concentration. Cluster analysis showed high similarity among sites occurred in low P conditions, and dissimilarity in taxonomic composition was relatively high among high P sites. This differs from other studies that show human disturbance causes greater homogenization of flora among sites (reduced beta diversity, Lougheed et al., 2008 ). The mechanisms regulating apparent and absolute diversity probably vary. The low P conditions in the greater Everglades that constrain taxonomic membership to a small number of taxa may be an unusually low resource condition not observed along human disturbance gradients in other locations. Thus, in the greater Everglades and BCNP we have low nutrient supply constraining taxa numbers with release of nutrient constraint as P supply increases allowing other taxa to invade, thereby increasing apparent numbers of taxa in the habitat. In the higher range of a resource gradient that can be found in other non-Everglades locations with high human disturbance (Pan et al., 1996; Wang et al., 2006 ; Lougheed et al, 2007 ), nutrients may be sufficiently high at the high end of the human disturbance gradient so additional disturbance allows overgrowth and dominance of a few high nutrient species, thereby homogenizing biota with increasing human disturbance. Taxa traits and relative metric performance More taxa were indicative of high P than low P conditions in BCNP. This is likely related to the low mat P conditions at some BCNP sites, which are among the lowest in the greater Everglades area. All taxa identified as low P and high P with TITAN and regression were characterized the same. In contrast, there were discrepancies for low and high P characterizations from the literature evaluations by NPS for BCNP when compared to the TITAN and regression traits determined with the calibration dataset. Some of the discrepancies between NPS literature and BCNP data evaluations were due to different taxonomic treatments of taxa, but many were related to the literature dealing with higher nutrient ranges, so low P (oligotrophic) taxa in the literature were from relatively low P conditions in studies which happen to be intermediate or high P conditions for BCNP and the greater Everglades. Taxa P optima were highly correlated among different Everglades projects (Gaiser et al., 2006 ; Slate & Stevenson 2007 ; Hagerthey et al., 2012 ) despite the limited number of pairwise taxa comparisons across studies that were possible. The concordance of these multiple lines of evidence again indicate that phosphorus is a causal determinant of low and high P characterizations assigned to taxa. Almost all metrics were strongly related to mat P, whether traits were derived by regression and TITAN with BCNP data or by reviewing literature. Almost all metrics with all traits differed among BCNP regions with differing levels of human disturbance. However, metrics using traits derived using regression and TITAN consistently performed better than metrics calculated with traits derived from the NPS literature evaluation. Tests with the validation dataset and the 2020 regional comparison indicated that regression and TITAN traits performed better in metrics than literature traits. As in other studies, project specific taxonomy may affect trait assignments to species derived from the literature, but we also show environmental gradient lengths and ranges also affect trait assignment and metric performance. In a comparison of regression and TITAN, TITAN traits consistently performed better. Metrics calculated with the relative log abundance of taxa performed better than metrics calculated as the proportion of individuals or taxa. This was observed with the validation data tests of metrics relationships to mat P and differences among regions with differing mat P. Applications of simpler metrics calculated with the proportion of taxa having a trait may have good precision (low variability) when related to disturbance (Stevenson et al., 2008 b) and can be easier to relate conceptually to standard characteristics of biological condition as species loss or invasion (Davies et al., 2006), however resource managers should know that accuracy (proximity to true value) of these metrics is a concern. The actual proportions of low or high disturbance taxa could be very different if more thorough assessments than 600 valve counts were used to determine taxon presence or absence, as illustrated by the large numbers of rare taxa in log-normal, taxon-abundance distributions (Preston, 1948 ; Patrick, 1967 ). Traits based metrics for pollutants can be more precise indicators of pollution levels than measurements of pollutants that are difficult to measure in the environment because of spatial and temporal variability (Stevenson et al., 2010 ). In BCNP, WAM and metrics using either TITAN or regression-derived taxa P traits were more precisely related to the human disturbance gradient than measured mat P. Additionally, P concentrations in algae (mat P) or sediments are more precise and sensitive indicators of wetland P availability than measurements of P concentration in the water column (Pan et al., 2000 ; Gaiser et al., 2004 ). Thus, basing management targets on biological indicators, either as measures of biological condition or as indicators of P pollution, warrants consideration. Metrics were affected little by natural factors such as wetland habitat type, substrate type, substrate location, water depth, and when sampling occurred during the high water season. Analyses of both sets of residuals, for metric-mat P relationships or metric and mat P relationships with the lat-long models of human disturbance, showed little likelihood that natural factors had major effects on metrics. Metrics are often robust to pollution levels and human disturbance, but there are study scales when we need to account for how minimally disturbed condition of pollutants and metrics vary with natural factors, such as in streams and lakes assessments (Cao et al., 2007 ; Stevenson et al., 2013 ; Tang et al., 2016 ). But the range of natural factors in BCNP had little effect on metric relationships to mat P (the stressor) or the lat-long model of human disturbance. Hydrologic year did affect expected levels of mat P and metrics in BCNP according to differences in residuals for the metric and mat P relationships with the lat-long model. All years of samples were counted by the same taxonomist for data used in metric testing and evaluation, which reduces the likelihood that taxonomist caused interannual variability in residuals. Interannual differences in P availability seem to be the issue because mat P residuals also differed with hydrologic year. Perhaps weather-related events, such as long periods of stable weather or recent rainfalls, affected P availability and thereby, metric values. Management benchmarks and metric considerations Management benchmarks for protecting ecological condition call for determining how ecological systems respond to pollutants and relating those responses to different levels of human disturbance to prevent over- or under-protection of a resource. We were particularly interested in thresholds (tipping points) of ecological responses for their potential to develop consensus among stakeholders for management targets (Muradian, 2001 ), such as levels of biological condition to protect and corresponding levels of pollution that protect biological condition. Many metrics had nonlinear relationships with untransformed mat P concentrations, either increasing or decreasing asymptotically. These responses are not as helpful for identifying threshold responses as are responses having some assimilative capacity at low levels of stressors where biota respond little because they are so highly sensitive to even small changes in human disturbance (Fig. S1 , Stevenson, 2011 ). Decision tree analyses did show biological responses with assimilative capacity, such as high P taxa increasing when mat P surpassed 137 µg g − 1 mat P, and low P taxa decreasing when mat P exceeded 464 µg g − 1 . These mat P concentrations showing assimilative capacity provide breakpoints in relationships that are particularly valuable for justifying management targets in BCNP. In addition, groups of sites with similar diatom assemblage composition were observed in low P conditions, e.g. Group 1–1 was the low P subgroup within Group 1. The 75th percentile of mat P for Group 1–1 was 142 and for all of Group 1 (that included subgroups 1–1, 1–2, and 1–3) was 354 µg g − 1 . TITAN analysis showed transition from low to high P taxa occurring in the 300–600 µg g − 1 mat P range with most transition of low P taxa in the 300–400 µg g − 1 range. Breakpoints in ecological responses of calcareous periphyton in the eastern, peat-based wetlands of the Everglades could be relevant for management of BCNP because calcareous periphyton with the same species is also expected in minimally disturbed conditions of BCNP. Much of the early research in the Everglades relates periphyton responses to water column TP. Floating mat cover and species composition changed greatly with elevation of TP above 10 µg L − 1 in WCA-2A (McCormick, et al. 1996 ; Stevenson et al., 2002 ), which was important for justifying Florida’s 10 µg TP L − 1 criterion for the Everglades. According to Hagerthey et al. ( 2008 ), mat P would be 412 µg g − 1 in sloughs with a water column TP of 11 µg L − 1 , which was a threshold for slough macrophytes. Pan et al. ( 2000 ) found a characteristic low P periphyton assemblage (TWINSPAN I) at sites with epiphytic mat P averaging 290 µg g − 1 ± 50 SE, which with n = 10 would have a 75th mat P percentile of 300 µg g − 1 . Gaiser et al. ( 2005 ) found enrichment of background conditions in long-term experiments by as little as 5 µg P L − 1 can alter periphyton assemblages, and therefore, the metabolic biogeochemical processes that regulate periphyton. This sensitivity to low concentrations of P is not surprising given the asymptotic increase in periphyton growth rates and high sensitivity when P is less than 10 µg L − 1 in stream water (Bothwell, 1989 ; Rier et al., 2006). The substantial transitions from low to high P taxa in BCNP periphyton and in the Everglades occurring in the 300–400 µg g − 1 mat P range corresponded to the 75th percentile of mat P in minimally disturbed conditions of BCNP, which was 367 µg g − 1 . Differences in periphyton in ranges of mat P less than the 75th percentile of mat P in the three most southern and eastern minimally disturbed regions of BCNP indicated that habitat specific management targets could be important. For example, if the mat P management target were set at the 367 µg g − 1 mat P, then Group 1–1 assemblages found in very low P conditions would not be protected, even though the more inclusive Group 1 assemblages would be protected. Gaiser et al. ( 2006 ) characterized minimally disturbed conditions for Shark River Slough and Taylor Slough of the Everglades as having maximum TP at levels that would protect the Group 1–1 periphyton assemblages, so such habitats do occur on a large scale in southern Florida. Further research should be conducted to determine if lower mat P management targets should be used to protect some habitat types with the Group 1–1 periphyton assemblages. Natural habitat features could create different P conditions, which could justify setting different management targets for different kinds of wetland habitats, even though our analyses showed residuals in the metric mat P relationships with lat-long were not related to natural factors (e.g. local hydrologic conditions). Many factors should be considered when selecting metrics for use in a monitoring program (Jackson et al., 2000 ). All BCNP metrics were related to mat P, differed among BCNP regions in 2020, and were related to the human disturbance gradient, in most cases. Therefore, many options could be satisfactory for use in a monitoring program. Distinguishing between metrics measuring an element of biological condition (low or high P taxa indicators) versus inferring stressors (WAM mat P models) is important conceptually relative to the idea that we want to manage stressors to protect biological condition and have clear communication with stakeholders about levels of protection for measures of biological condition (Stevenson, 2011 ). The higher performance of TITAN and RlogA metrics of biological condition introduces tradeoffs between using metrics that perform best (i.e., are the most accurately and precisely related to mat P concentration, changes among regions, and human disturbance) and metrics that are simpler to explain to stakeholders, i.e. how they were calculated and relate to ecological condition. Simple metrics such as numbers or even proportions of low and high P taxa may also be more vulnerable to taxonomic error, if the number of taxa distinguished among taxonomists differs. There are relatively simple ways to explain TITAN and RlogA methods and why they are likely better than other methods, but nevertheless, they are more complicated methods. Fortunately, a library of taxa traits and a database of diatom counts allows calculation of all the metrics relatively easily, so some metrics can be used for public reports and others can be used for administrative and technical records and analysis. Such a set of tools, as we provided with this research, provides the necessary information for advancing a monitoring and assessment program and establishing management targets for biological condition and stressor levels that will protect an unimpaired condition, which is the specified goal of BCNP management (Patterson et al., 2008 ). Conclusions We have developed the information needed for a monitoring and assessment program and integrated it for BCNP with many attributes that are not included in other programs. First we explored the ecological relationships among indicators of biological condition and pollutants which are directly related to the goals of management of the resource. Then we developed multiple metrics of biological condition and tested which are the most sensitive and consistently related to the primary pollutant. We confirmed the cause-effect relationship between metrics of biological condition and P, the primary pollutant. Then we determined the natural variation in minimally disturbed condition and characterized changes in biological condition from minimally disturbed condition with increasing P concentrations which can be used to justify effects-based management targets. These are fundamental steps for developing an effects-based monitoring and assessment program, but all of them are seldom integrated and presented in one publication with state-of-the-art approaches to show clear relationships among the lines of results needed for an integrated monitoring and assessment program. We advanced the state-of-the-art of monitoring and assessment at multiple steps in our research. We showed how nutrient enrichment clearly affected biological assemblages by allowing invasion of taxa that cannot survive in the naturally low nutrient conditions and how sensitive low nutrient taxa are greatly reduced or eliminated as nutrient concentrations increase. We confirmed cause-effect relationships using multiple causal assessment methods. We developed and tested a new indicator calculation method that improved metric performance and was based on basic principles of population growth dynamics using log transformed relative abundances. And finally, we integrated the application of minimally disturbed condition and biological condition response along a pollution gradient to provide a potentially ultimate management goal of minimally disturbed condition, but also interim goals to assess incremental restoration of the ecosystem. Declarations Funding The United States National Park Service funded this research. Author Contribution K.R.T.W. was responsible for N.P.S. administration of the project. K.R.T.W. and M.C.P. collected samples. R.J.S. analyzed samples and data, and was lead writer of the manuscript. All authors contributed to writing and reviewing the manuscript. Acknowledgement The National Park Service (NPS) funded this research. Raul Urgelles and Andrea Atkinson at NPS were founding members of this monitoring project. Mario Londono at the NPS provided database assistance. The NPS Big Cypress Aviation unit supported sampling missions; and specially we would like to thank Fred Goodwin (pilot) and Michael Roof (Helicopter Program Manager). Data Availability The data that support the findings of this study are not openly available because they are being prepared for deposit in an appropriate archive. They will be published before this manuscript is accepted for publication. References Baker, M. E., & King, R. S. (2010). A new method for detecting and interpreting biodiversity and ecological community thresholds. Methods in Ecology and Evolution, 1 (1), 25-37. doi:https://doi.org/10.1111/j.2041-210X.2009.00007 Beyers, D. W. (1998). Causal inference in environmental impact studies. Journal of the North American Benthological Society, 17 , 367-373. Bothwell, M. L. (1989). Phosphorus-limited growth dynamics of lotic periphytic diatom communities: areal biomass and cellular growth rate responses. Canadian Journal of Fisheries and Aquatic Sciences, 46 , 1293-1301. Browder, J. A., Gleason, P. J., & Swift, D. R. (1994). Periphyton in the Everglades: spatial variation, environmental correlates,and ecological implications. In S. M. Davis & J. C. Ogden (Eds.), Everglades, the ecosystem and its restoration (pp. pp.379-418). Delray Beach, Florida: St. Lucie Press. Cao, Y., Hawkins, C. P., Olson, J., & Kosterman, M. A. (2007). Modeling natural environmental gradients improves the accuracy and precision of diatom-based indicators. Journal of the North American Benthological Society, 26 (3), 566-585. Cardinale, B. J., Hillebrand, H., Harpole, W. S., Gross, K., & Ptacnik, R. (2009). Separating the influence of resource ‘availability’ from resource ‘imbalance’ on productivity–diversity relationships. Ecology Letters, 12 (6), 475-487. doi:https://doi.org/10.1111/j.1461-0248.2009.01317 Carlisle, D. M., Spaulding, S. A., Tyree, M. A., Schulte, N. O., Lee, S. S., Mitchell, R. M., & Pollard, A. A. (2022). A web-based tool for assessing the condition of benthic diatom assemblages in streams and rivers of the conterminous United States. Ecological Indicators, 135 , 108513. doi:https://doi.org/10.1016/j.ecolind.2021.108513 Charles, D. F., Kelly, M. G., Stevenson, R. J., Poikane, S., Theroux, S., Zgrundo, A., & Cantonati, M. (2021). Benthic algae assessments in the EU and the US: Striving for consistency in the face of great ecological diversity. Ecological Indicators, 121 , 107082. doi:https://doi.org/10.1016/j.ecolind.2020.107082 Cormier, S. M., & Suter II, G. W. (2008). A framework for fully integrating environmental assessment. Environmental Management, 42 , 543-556. Davies, S. P., & Jackson, S. K. (2006). The biological condition gradient: a descriptive model for interpreting change in aquatic ecosystems. Ecological Applications, 16 , 1251-1266. Gaiser, E. E., Childers, D. L., Jones, R. D., Richards, J. H., & Scinto, L. J. (2006). Periphyton responses to eutrophication in the Florida Everglades: cross-system patterns of strucural and composition changes. Limnology & Oceanography, 51 , 617-630. Gaiser, E. E., L. J. Scinto, J. H. Richards, K. Jayachandran, D. L. Childers, J. D. Trexler, & Jones, R. D. (2004). Phosphorus in periphyton mats provides best metric for detecting low-level P enrichment in an oligotrophic wetland. Water Research, 38 , 507-516. Gaiser, E. E., Trexler, J. C., Richards, J. H., Childers, D. L., Lee, D., Edwards, A. L., . . . Jones, R. D. (2005). Cascading Ecological Effects of Low-Level Phosphorus Enrichment in the Florida Everglades. Journal of Environmental Quality, 34 (2), 717-723. doi:https://doi.org/10.2134/jeq2005.0717 Gleason, P. J., & W. Spackman, J. (1974). Calcareous periphyton and water chemistry in the Everglades. In P. J. Gleason (Ed.), Environments of South Florida: Present and Past, Memoir No. 2. (pp. pp. 146-181). Caral Gables, Fla.: Miami Geological Society. Hagerthey, S. E., Bellinger, B. J., Wheeler, K., Gantar, M., & Gaiser, E. (2011). Everglades Periphyton: A Biogeochemical Perspective. Critical Reviews in Environmental Science and Technology, 41 (sup1), 309-343. doi:10.1080/10643389.2010.531218 Hagerthey, S. E., Newman, S., Rutchey, K., Smith, E. P., & Godin, J. (2008). Multiple regime shifts in a subtropical peatland: community-specific thresholds in eutrophication. Ecological Monographs, 78 (4), 547-565. doi:https://doi.org/10.1890/07-0538. Hagerthey, S. E., Newman, S., & Xue, S. (2012). Periphyton-based transfer functions to assess ecological imbalance and management of a subtropical ombrotrophic peatland. Freshwater Biology, 57 , 1947-1965. Jackson, L. E., Kutz, J. C., & Fisher, W. S. (Eds.). (2000). Evaluation guidelines for ecological indicators . Research Triangle Park, North Carolina: United States Environmental Protection Agency. Jarosewich, M. & Wagner, J. (1985) Geologic Structure of the Surficial Aquifer System Underlying Everglades National Park and the Big Cypress National Preserve. South Florida Research Center Report SFRC-85/02. 43p. La Hée, J. M., & Gaiser, E. E. (2012). Benthic diatom assemblages as indicators of water quality in the Everglades and three tropical karstic wetlands. Freshwater Science, 31 , 205-221. Lavoie, I., Hamilton, P. B., Wang, Y. K., Dillon, P. J., & Campeau, S. (2009). A comparison of stream bioassessment in Québec (Canada) using six European and North American diatom-based indices. Nova Hedwigia, 135 , 37-56. Londoño, M. 2019. Selecting Sites in GIS—Version 1.0. South Florida/Caribbean Network Standard Operating Procedure. South Florida/Caribbean Network, Palmetto Bay, Florida. Available at: https://irma.nps.gov/DataStore/Reference/Edit/2260145 (last accessed April 2019). Lougheed, V. L., McIntosh, M. D., Parker, C. A., & Stevenson, R. J. (2008). Wetland degradation leads to homogenization of biota at local and landscape scales. Freshwater Biology, 52 , 2402-2413. Lougheed, V. L., Parker, C. A., & Stevenson, R. J. (2007). Using non-linear responses of multiple taxonomic groups to establish criteria indicative of wetland biological condition. Wetlands, 27 , 96-109. McCormick, P. V., Chimney, M. J., & Swift, D. R. (1997). Diel oxygen profiles and water column community metabolism in the Florida Everglades, USA. Archiv Fur Hydrobiologie, 140 (1), 117-129. McCormick, P. V., O'Dell, M. B., Shuford, R. B. E., Backus, J. G., & Kennedy, W. C. (2001). Periphyton responses to experimental phosphorus enrichment in a subtropical wetland. Aquatic Botany, 71 (2), 119-139. McCormick, P. V., Rawlick, P. S., Lurding, K., Smith, E. P., & Sklar, F. H. (1996). Periphyton-water quality relationships along a nutrient gradient in the northern Florida Everglades. Journal of the North American Benthological Society, 15 , 433-449. Miller, R., McPherson, B., Sobczak, R., & Clark, C. (2004) Water Quality in Big Cypress National Preserve and Everglades National Park — Trends and Spatial Characteristics of Selected Constituents. U.S. Geological Survey (USGS) Water Resources Investigations Report. 03-4249. U.S. Geological Survey; U.S. Department of the Interior. Denver, Colorado. Muradian, R. (2001). Ecological thresholds: a survey. Ecological Economics, 38 , 7-24. Muxo, R. & Shamblin, R. B. (2019) Standard operating procedure 4: Periphyton Sample Collection Methods. South Florida/Caribbean Network, Palmetto Bay, Florida. Available at: https://irma.nps.gov/DataStore/Reference/Edit/2260145 (last accessed April 2019). Oksanen, J., F. Guillaume Blanchet, M. Friendly, R. Kindt, P. Legendre, D. McGlinn, . . . Wagner, H. (2020). vegan: Community Ecology Package (Version R package version 2.5-7). Retrieved from https://CRAN.R-project.org/package=vegan Pan, Y., & Stevenson, R. J. (1996). Gradient analysis of diatom assemblages in Western Kentucky wetlands. Journal of Phycology, 32 , 222-232. Pan, Y., Stevenson, R. J., Vaithiyanathan, P., Slate, J., & Richardson, C. J. (2000). Changes in algal assemblages along observed and experimental phosphorus gradients in a subtropical wetland, U.S.A. Freshwater Biology, 44 , 339-353. Patrick, R. (1967). The effect of invasion rate, species pool, and size of area on the structure of the diatom community. Proceedings of the Academy of Natural Sciences, 58 , 1335-1342. Patterson, M. E., A.J. Atkinson, B.D. Witcher, K.R.T. Whelan, W.J. Miller, R.J. Waara, . . . Shamblin, R. B. (2008). South Florida / Caribbean Network vital signs monitoring plan. . Retrieved from Pillsbury, R., Stevenson, R. J., Munn, M. D., & Waite, I. (2019). Relationships between diatom metrics based on species nutrient traits and agricultural land use. Environmental Monitoring and Assessment, 2019 (191). doi:https://doi.org/10.1007/s10661-019-7357-8 Potapova, M., & Charles, D. (2002). Benthic diatoms in USA rivers: distributions along spatial and environmental gradients. Journal of Biogeography, 29 , 167-187. Preston, F. W. (1948). The commonness, and rarity, of species. Ecology, 29 , 254-283. R Core Team (2021). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL https://www.R-project.org/. Raschke, R. L. (1993). Diatom (Bacillariophyta) community response to phosphorus in the Everglades National Park USA. Phycologia, 32 , 48-58. Reavie, E. D., Kireta, A. R., Kingston, J. C., Sgro, G. V., Danz, N. P., Axler, R. P., & Hollenhorst, T. P. (2008). Comparison of simple and multimetric diatom-based indicators for Great Lakes coastline disturbance. Journal of Phycology, 44 , 787-802. Rier, S. T., & Stevenson, R. J. (2006). Response of periphytic algae to gradients in nitrogen and phosphorus in streamside mesocosms. Hydrobiologia, 561 , 131-147. Slate, J. E., & Stevenson, R. J. (2007). The diatom flora of phosphorus-enriched and unenriched sites in an Everglades marsh. Diatom Research, 22 , 355-386. doi:10.1080/0269249X.2007.9705721 Smith, S. (2018) Total Phosphorus (TP) Protocol, Seagrass Ecosystems Research Laboratory (SERL) Coastlines and Oceans Division, Institute of Environment, Florida International University. Solorzano, L. & Sharp, J. H. (1980) Determination of total dissolved phosphorus and particulate phosphorus in natural waters. Limnology and Oceanography 25: 754–758. Soranno, P. A., Cheruvelil, K. S., Stevenson, R. J., Rollins, S. L., Holden, S. W., Heaton, S., & Torng, E. (2008). A framework for developing ecosystem-specific nutrient criteria: Integrating biological thresholds with predictive modeling. Limnology and Oceanography, 53 , 773-787. Stevenson, R. J. (2006). Refining diatom indicators for valued ecological attributes and development of water quality criteria. In N. Ognjanova-Rumenova & K. Manoylov (Eds.), Advances in Phycological Studies (pp. 365-383). Moscow, Russia: Pensoft Publishers. Stevenson, R. J. (2011). A revised framework for coupled human and natural systems, propagating thresholds, and managing environmental problems. Physics and Chemistry of the Earth, 36 , 342-351. Stevenson, R. J. (2014). Ecological assessment with algae: a review and synthesis. Journal of Phycology, 50 , 437-461. Stevenson, R. J., Bailey, R. C., Harass, M. C., Hawkins, C. P., Alba-Tercedor, J., Couch, C., . . . Johnson, R. K. (2004). Designing data collection for ecological assessments. In M. T. Barbour, S. B. Norton, H. R. Preston, & K. W. Thornton (Eds.), Ecological Assessment of Aquatic Resources: Linking Science to Decision-Making (pp. 55-84). Pensacola, Florida: Society of Environmental Toxicology and Contamination Publication. Stevenson, R. J., Hill, B. E., A.T. Herlihy, Yuan, L. L., & Norton, S. B. (2008). Algal-P relationships, thresholds, and frequency distributions guide nutrient criterion development. Journal of the North American Benthological Society, 27 , 783-799. Stevenson, R. J., Pan, Y., Manoylov, K., Parker, C., Larsen, D. P., & Herlihy, A. T. (2008). Development of diatom indicators of ecological conditions for streams of the western United States. Journal of the North American Benthological Society, 27 , 1000-1016. Stevenson, R. J., Pan, Y., & Vaithiyanathan., P. (2002). Ecological assessment and indicator development in wetlands: the case of algae in the Everglades, USA. Verhandlungen Internationale Vereinigung für Theoretische und Andgewandte Limnologie, 28 , 1248-1252. Stevenson, R. J., Pan, Y., & van Dam, H. (2010). Assessing ecological conditions in rivers and streams with diatoms. In J. P. Smol & E. F. Stoermer (Eds.), The Diatoms: Applications to the Environmental and Earth Sciences (Second ed., pp. 57-85). Cambridge, UK: Cambridge University Press. Stevenson, R. J., Zalack, J., & Wolin, J. (2013). A multimetric index of lake diatom condition using surface sediment assemblages. Freshwater Science, 32 , 1005–1025. Stoddard, J. L., Larsen, D. P., Hawkins, C. P., Johnson, R. K., & Norris, R. H. (2006). Setting expectations for the ecological condition of streams: the concept of reference condition. Ecological Applications, 16 , 1267-1276. Tang, T., Stevenson, R. J., & Grace, J. B. (2020). The importance of natural versus human factors for ecological conditions of streams and rivers. Science of the Total Environment, 704 , 135268. doi:https://doi.org/10.1016/j.scitotenv.2019.135268 Tang, T., Stevenson, R. J., & Infante, D. M. (2016). Accounting for regional variation in both natural environment and human disturbance to improve performance of multimetric indices of lotic benthic diatoms. Science of the Total Environment 268 , 1124-1134. Therneau, T, & Atkinson, B. (2025). rpart: Recursive Partioning and Regression Trees (Version R package (4.1.24). DOI 10.32614/CRAN.package.rpart Urgelles, R., Whelan, K. R. T., Muxo, R., Shamblin, R. B., Patterson, J. M., & Atkinson, A. J. (2019) Periphyton monitoring in Big Cypress National Preserve: Protocol narrative. Natural Resource Report. NPS/SFCN/NRR—2019/1911. National Park Service. Fort Collins, Colorado USEPA. (1992). Framework for Ecological Risk Assessment (EPA/630/R-92/001) United States Environmental Protection Agency, Washington, DC: USEPA. (2000). Stressor Identification Guidance Document (EPA 600/R/96/059) United States Environmental Protection Agency, Washington, DC: Van Dam, H., Mertenes, A., & Sinkeldam, J. (1994). A coded checklist and ecological indicator values of freshwater diatoms from the Netherlands. Netherlands Journal of Aquatic Ecology, 28 , 117-133. Wang, Y. K., Stevenson, R., Sweets, P., & DiFranco, J. (2006). Developing and Testing Diatom Indicators for Wetlands in the Casco Bay Watershed, Maine, USA. Hydrobiologia, 561 (1), 191-206. Retrieved from http://springerlink.metapress.com/openurl.asp?genre=article&id=doi:10.1007/s10750-005-1614-2 Wang, Y. K., Stevenson, R. J., & Metzmeier, L. (2005). Development and evaluation of a diatom-based index of biotic integrity for the Interior Plateau Ecoregion. Journal of the North American Benthological Society, 24 , 990-1008. Whelan, K. R. T., Prats, M., & Atkinson, A. J. (2020) Western Big Cypress National Preserve vegetation map. Natural Resource Report NPS/SFCN/NRR—2020/2170. National Park Service, Fort Collins, Colorado. https://doi.org/10.36967/nrr-2278240. Wilson, S. (2020) BCAL phosphorus sample preparation protocol., Blue Carbon Analysis Laboratory and Seagrass Ecosystems Research Laboratory, Coastlines and Oceans Worm, B., Lotze, H. K., Hillebrand, H., & Sommer, U. (2002). Consumer versus resource control of species diversity and ecosystem functioning. Nature, 417 (6891), 848-851. Wyatt, K. H., Stevenson, R. J., & Turetsky, M. R. (2010). The importance of nutrient co-limitation in regulating algal community composition, productivity, and algal-derived DOC in an oligotrophic marsh in interior Alaska Freshwater Biology 55 , 1845-1860. Additional Declarations No competing interests reported. Supplementary Files TablesandFiguresSupplemental20250206.pdf Cite Share Download PDF Status: Published Journal Publication published 15 Mar, 2026 Read the published version in Environmental Monitoring and Assessment → Version 1 posted Editorial decision: Revision requested 02 Aug, 2025 Reviews received at journal 13 Jul, 2025 Reviewers agreed at journal 08 Jul, 2025 Reviewers agreed at journal 08 Jul, 2025 Reviewers invited by journal 07 Jul, 2025 Editor assigned by journal 01 Jul, 2025 Submission checks completed at journal 01 Jul, 2025 First submitted to journal 23 Jun, 2025 You are reading this latest preprint version Research Square lets you share your work early, gain feedback from the community, and start making changes to your manuscript prior to peer review in a journal. As a division of Research Square Company, we’re committed to making research communication faster, fairer, and more useful. We do this by developing innovative software and high quality services for the global research community. Our growing team is made up of researchers and industry professionals working together to solve the most critical problems facing scientific publishing. Also discoverable on Platform About Our Team In Review Editorial Policies Advisory Board Help Center Resources Author Services Accessibility API Access RSS feed Manage Cookie Preferences © Research Square 2026 | ISSN 2693-5015 (online) Privacy Policy Terms of Service Do Not Sell My Personal Information {"props":{"pageProps":{"initialData":{"identity":"rs-6959033","acceptedTermsAndConditions":true,"allowDirectSubmit":false,"archivedVersions":[],"articleType":"Research Article","associatedPublications":[],"authors":[{"id":483219934,"identity":"d69df580-6744-4c02-bf8f-aaac4c05251c","order_by":0,"name":"R. Jan Stevenson","email":"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAZAAAAAyAQMAAABI0h/eAAAABlBMVEX///8AAABVwtN+AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAAs0lEQVRIiWNgGAWjYFCCAyDCBsLmIVILYwMDQxpJWhhAWg6ToMW88fjzBx/bzkfLz0hgfPC2jQgtMgfOGDbObLudu+FGArPhXGK0SDCcYWzmOQPUIpHAJs1LnJbjD5v/nDmXO39GAvtvIrUcMGxmqDiQ23AjgY2ZSC1nDGf2VCTnbjjzsFlyzjlitEgcf/Dhh4Fd7vz25IMf3pQRoYVB4gCMBYofogA/sQpHwSgYBaNg5AIAPas7vKy7ux8AAAAASUVORK5CYII=","orcid":"","institution":"Michigan State University","correspondingAuthor":true,"prefix":"","firstName":"R.","middleName":"Jan","lastName":"Stevenson","suffix":""},{"id":483219935,"identity":"dd541f6b-39d3-4de8-8af0-619d1fa10712","order_by":1,"name":"Kevin R. T. Whelan","email":"","orcid":"","institution":"National Park Service","correspondingAuthor":false,"prefix":"","firstName":"Kevin","middleName":"R. T.","lastName":"Whelan","suffix":""},{"id":483219936,"identity":"a6d597fd-7446-414c-9a25-8777ce9e08d8","order_by":2,"name":"Michelle C. Prats","email":"","orcid":"","institution":"National Park Service","correspondingAuthor":false,"prefix":"","firstName":"Michelle","middleName":"C.","lastName":"Prats","suffix":""}],"badges":[],"createdAt":"2025-06-23 17:38:18","currentVersionCode":1,"declarations":"","doi":"10.21203/rs.3.rs-6959033/v1","doiUrl":"https://doi.org/10.21203/rs.3.rs-6959033/v1","draftVersion":[],"editorialEvents":[{"content":"https://doi.org/10.1007/s10661-026-15144-0","type":"published","date":"2026-03-15T15:58:41+00:00"}],"editorialNote":"","failedWorkflow":false,"files":[{"id":86521324,"identity":"f02490f9-78a1-4851-8e11-00d9f7822234","added_by":"auto","created_at":"2025-07-11 15:07:24","extension":"png","order_by":1,"title":"Figure 1","display":"","copyAsset":false,"role":"figure","size":1458253,"visible":true,"origin":"","legend":"\u003cp\u003eLeft panel: Map of BCNP and the study area in northwest corner. Right panel: Names of regions within the study area.\u003c/p\u003e","description":"","filename":"1.png","url":"https://assets-eu.researchsquare.com/files/rs-6959033/v1/ab8ff50fc932fafa3243728e.png"},{"id":86521325,"identity":"887b35c8-b95a-485d-b60c-2b54e735b322","added_by":"auto","created_at":"2025-07-11 15:07:24","extension":"png","order_by":2,"title":"Figure 2","display":"","copyAsset":false,"role":"figure","size":110205,"visible":true,"origin":"","legend":"\u003cp\u003eBoxplots of dissimilarity, mat P, pH, and conductivity for the groups of sites delineated in cluster analysis of the calibration data set.\u003c/p\u003e","description":"","filename":"2.png","url":"https://assets-eu.researchsquare.com/files/rs-6959033/v1/9b9f2333caa48499f954b00a.png"},{"id":86521326,"identity":"a842d5b2-8e64-4d77-b7d6-367dc2fd1c4d","added_by":"auto","created_at":"2025-07-11 15:07:24","extension":"png","order_by":3,"title":"Figure 3","display":"","copyAsset":false,"role":"figure","size":55590,"visible":true,"origin":"","legend":"\u003cp\u003eRelationship between species richness and log mat P. Levels of mat P associated with changepoints having the lowest probability of occurring by chance are shown by vertical arrows (Table S9).\u003c/p\u003e","description":"","filename":"3.png","url":"https://assets-eu.researchsquare.com/files/rs-6959033/v1/868f97200cfdd7203ea3edcd.png"},{"id":86521327,"identity":"9cc7e854-1fe4-4c98-a030-a06f0e7584da","added_by":"auto","created_at":"2025-07-11 15:07:24","extension":"png","order_by":4,"title":"Figure 4","display":"","copyAsset":false,"role":"figure","size":53323,"visible":true,"origin":"","legend":"\u003cp\u003eA stacked bar graph of the dominant species in samples in groups of samples with increase mat P concentrations. Braocala = \u003cem\u003eBrachysira ocalanensis\u003c/em\u003e, Enceverg = \u003cem\u003eEncyonema evergladianum\u003c/em\u003e, Encmicro = \u003cem\u003eEncyonema microcephala\u003c/em\u003e, Encsiles = \u003cem\u003eEncyonema silesiacum\u003c/em\u003e, Gomaurit = \u003cem\u003eGomphonema auritum\u003c/em\u003e, Mascalca = \u003cem\u003eMastogloia calcarea\u003c/em\u003e, Navcrytn = \u003cem\u003eNavicula cryptotenella\u003c/em\u003e, and Nitamphi = \u003cem\u003eNitzschia amphibia\u003c/em\u003e.\u003c/p\u003e","description":"","filename":"4.png","url":"https://assets-eu.researchsquare.com/files/rs-6959033/v1/ded87002908278ffb70d2732.png"},{"id":86521333,"identity":"cb21fc02-bac6-4017-9970-9d655cacb0b7","added_by":"auto","created_at":"2025-07-11 15:07:24","extension":"png","order_by":5,"title":"Figure 5","display":"","copyAsset":false,"role":"figure","size":148541,"visible":true,"origin":"","legend":"\u003cp\u003eRelationship between TITAN trait metrics and log mat P. The metric-mat P relationships for all other metrics are in Figs. S8 and S9.\u003c/p\u003e","description":"","filename":"5.png","url":"https://assets-eu.researchsquare.com/files/rs-6959033/v1/d0c812cc590d884b3db92262.png"},{"id":86521331,"identity":"c82b0da1-65fb-47de-af1a-8aff0e58fad3","added_by":"auto","created_at":"2025-07-11 15:07:24","extension":"png","order_by":6,"title":"Figure 6","display":"","copyAsset":false,"role":"figure","size":295933,"visible":true,"origin":"","legend":"\u003cp\u003eAdjusted R\u003csup\u003e2\u003c/sup\u003e and ranks of metric performance for metrics with different traits. Panels A-C show effects of metric type, trait source, and trait on adj. R\u003csup\u003e2\u003c/sup\u003e of metric-mat P relationships. Panels D-F plot the ranks of each adj. R\u003csup\u003e2\u003c/sup\u003e of metric-mat P relationships to show effects of metric type, trait source, and trait on metric-mat P relationships, using rank of metric-mat P relationships for metric type (ranks 1-4) within trait-source-trait groups, for trait source (ranks 1-3) within metric type-trait groups, and trait (ranks 1-2) within trait source-metric type groups.\u0026nbsp;\u003c/p\u003e","description":"","filename":"6.png","url":"https://assets-eu.researchsquare.com/files/rs-6959033/v1/294664cb0eb5445bcd3854cc.png"},{"id":86522267,"identity":"8a16c4e2-a4c3-4031-96e3-72a53bbf710d","added_by":"auto","created_at":"2025-07-11 15:15:24","extension":"png","order_by":7,"title":"Figure 7","display":"","copyAsset":false,"role":"figure","size":277805,"visible":true,"origin":"","legend":"\u003cp\u003eBoxplots of selected metrics for regions in BCNP sampling area for 2020 samples. Region codes are BI = Bear Island, EC = East Crossing, EH =East Hinson Marsh, FP = Fire Prairie, LM = Little Marsh, MN = Monument, and OK = Okaloacoochee Slough.\u003c/p\u003e","description":"","filename":"7.png","url":"https://assets-eu.researchsquare.com/files/rs-6959033/v1/c62699cce1a5ecacd91c995d.png"},{"id":86521329,"identity":"7cf726a9-3530-4ce1-a2e2-130b3a79b37a","added_by":"auto","created_at":"2025-07-11 15:07:24","extension":"png","order_by":8,"title":"Figure 8","display":"","copyAsset":false,"role":"figure","size":131111,"visible":true,"origin":"","legend":"\u003cp\u003eRelationship between number of low and high P TITAN taxa and log mat P with delineations for changepoints determined with decision tree analysis. Levels of mat P associated with changepoints having the lowest probability of occurring by chance are shown by vertical arrows.\u003c/p\u003e","description":"","filename":"8.png","url":"https://assets-eu.researchsquare.com/files/rs-6959033/v1/9fc66aed25ddab97893a091a.png"},{"id":86523018,"identity":"80dbc246-4efa-4e7e-80c5-9cd102131368","added_by":"auto","created_at":"2025-07-11 15:23:24","extension":"png","order_by":9,"title":"Figure 9","display":"","copyAsset":false,"role":"figure","size":85236,"visible":true,"origin":"","legend":"\u003cp\u003eTITAN sum(z-) for low P taxa and sum(z+) for high P taxa for all candidate change points along the mat P gradient (labeled environmental gradient). The solid and dashed black lines show the cumulative frequency distribution of change points along the mat P gradient for low P (z-) and high P (z+) taxa, respectively.\u003c/p\u003e","description":"","filename":"9.png","url":"https://assets-eu.researchsquare.com/files/rs-6959033/v1/43165e119703ec22d63ef95f.png"},{"id":104739558,"identity":"dc8d56f9-768b-4afb-b79e-61352f0493dc","added_by":"auto","created_at":"2026-03-16 16:09:13","extension":"pdf","order_by":0,"title":"","display":"","copyAsset":false,"role":"manuscript-pdf","size":3032821,"visible":true,"origin":"","legend":"","description":"","filename":"manuscript.pdf","url":"https://assets-eu.researchsquare.com/files/rs-6959033/v1/68152803-5f4e-4f60-94ca-5231442b786a.pdf"},{"id":86522265,"identity":"6c96f32b-9ae2-4b15-910f-61f8882ef4d3","added_by":"auto","created_at":"2025-07-11 15:15:24","extension":"pdf","order_by":0,"title":"","display":"","copyAsset":false,"role":"supplement","size":1140116,"visible":true,"origin":"","legend":"","description":"","filename":"TablesandFiguresSupplemental20250206.pdf","url":"https://assets-eu.researchsquare.com/files/rs-6959033/v1/bc14b918b51ebdc61c4a9441.pdf"}],"financialInterests":"No competing interests reported.","formattedTitle":"Developing a monitoring and assessment program with diatoms, an improved metric calculation method, and causal analysis for Big Cypress National Preserve, Florida (USA)","fulltext":[{"header":"Introduction","content":"\u003cp\u003eBuilding an assessment program depends on management goals and an understanding of the ecosystem, from which variables are selected to test their application in assessment and for establishing management targets (USEPA, \u003cspan citationid=\"CR61\" class=\"CitationRef\"\u003e1992\u003c/span\u003e; Stevenson et al., \u003cspan citationid=\"CR50\" class=\"CitationRef\"\u003e2004\u003c/span\u003e; Patterson et al. \u003cspan citationid=\"CR35\" class=\"CitationRef\"\u003e2008\u003c/span\u003e). The basics of understanding the ecosystem inform which attributes are the target of management goals, what contaminants and habitat alterations (herein called stressors, \u003cem\u003esensu\u003c/em\u003e USEPA, \u003cspan citationid=\"CR61\" class=\"CitationRef\"\u003e1992\u003c/span\u003e) affect those management goal attributes, and what human activities produce those contaminants and habitat alterations (Stevenson, \u003cspan citationid=\"CR48\" class=\"CitationRef\"\u003e2011\u003c/span\u003e; Tang et al., \u003cspan citationid=\"CR57\" class=\"CitationRef\"\u003e2020\u003c/span\u003e). Measurements of the key direct and indirect variables in the ecosystem that regulate management goal attributes should be used to characterize the condition of the ecosystem relative to management goals. Assessment involves characterizing condition. Assessment can also call for comparing condition to benchmarks in variables that can be used to determine if ecosystem condition meets management goals and whether management actions are warranted. Diagnosing which contaminants and habitat alterations are most affecting management goal attributes can also be parts of assessments, so remedial management actions can be selected (USEPA, \u003cspan citationid=\"CR62\" class=\"CitationRef\"\u003e2000\u003c/span\u003e; Cormier et al., 2008).\u003c/p\u003e\u003cp\u003eThe goal of our research was to develop a periphyton monitoring program for Big Cypress National Preserve (BCNP). The overarching goal of the BCNP periphyton monitoring program is to support the United States National Park Service (NPS) mission to preserve the unimpaired condition of natural resources of National Parks (Patterson et al., \u003cspan citationid=\"CR35\" class=\"CitationRef\"\u003e2008\u003c/span\u003e). Periphyton have been studied extensively in greater Everglades area, which is located southeast of the BCNP. In the greater Everglades area and BCNP, a distinctive calcareous periphyton is observed in areas with minimal human disturbance (Browder et al., \u003cspan citationid=\"CR4\" class=\"CitationRef\"\u003e1994\u003c/span\u003e; La H\u0026eacute;e et al., 2012). Calcareous periphyton have important ecological functions, respond directly and sensitively to phosphorus pollution, and effects of phosphorus on periphyton cascade through the rest of the ecosystem and its biota (Gleason et al., 1974; McCormick et al., \u003cspan citationid=\"CR25\" class=\"CitationRef\"\u003e1997\u003c/span\u003e; McCormick et al., \u003cspan citationid=\"CR26\" class=\"CitationRef\"\u003e2001\u003c/span\u003e; Gaiser et al., \u003cspan citationid=\"CR13\" class=\"CitationRef\"\u003e2005\u003c/span\u003e; Gaiser et al., \u003cspan citationid=\"CR11\" class=\"CitationRef\"\u003e2006\u003c/span\u003e; Hagerthey et al., \u003cspan citationid=\"CR15\" class=\"CitationRef\"\u003e2011\u003c/span\u003e). For these reasons, periphyton were selected as one of the vital signs for BCNP (Patterson et al., \u003cspan citationid=\"CR35\" class=\"CitationRef\"\u003e2008\u003c/span\u003e) with park-specific goals to determine status and trends in periphyton, water quality, and ecosystem function (Urgelles et al., \u003cspan citationid=\"CR60\" class=\"CitationRef\"\u003e2019\u003c/span\u003e). Differences in soils, wetland transmissivity, and phosphorus range between the BCNP and the greater Everglades (Jarosewich \u0026amp; Wagner, \u003cspan citationid=\"CR19\" class=\"CitationRef\"\u003e1985\u003c/span\u003e) were deemed sufficient to warrant development and testing of a periphyton program specifically for BCNP.\u003c/p\u003e\u003cp\u003eOften, taxa traits for a new monitoring program are determined using past research in other regions. Lavoie et al. (\u003cspan citationid=\"CR21\" class=\"CitationRef\"\u003e2009\u003c/span\u003e) showed that taxa indicator traits from Van Dam et al. (\u003cspan citationid=\"CR63\" class=\"CitationRef\"\u003e1994\u003c/span\u003e), largely based on research and application in Europe, were related to pollution indicators in Canada. Tang et al. (\u003cspan citationid=\"CR58\" class=\"CitationRef\"\u003e2016\u003c/span\u003e) successfully used taxa traits developed for the western United States for assessments across the country. However, Potapova et al. (2002) found regional differences in indicator values across ecoregions in the U.S. Therefore, it seems reasonable to start developing a new program by using taxa trait characterizations from past projects in regions with similar types of human disturbance, pollution, and natural variability. Then, when sufficient data has been accumulated, as in BCNP, refine those taxa traits for the specific region that will be monitored.\u003c/p\u003e\u003cp\u003eMany types of metrics have been proposed, ranging from diversity measurements to indicator species abundances and weighted average models that infer stressor conditions (Stevenson \u003cspan citationid=\"CR49\" class=\"CitationRef\"\u003e2014\u003c/span\u003e). Given the goal of the BCNP program, metrics should be assessed for their ability to determine the deviation of the assessed ecosystem from minimally disturbed ecological condition. Widely accepted definitions of \u0026lsquo;minimally disturbed condition\u0026rsquo; and the \u0026lsquo;attributes of biological condition\u0026rsquo; have been proposed and serve as standards for many ecological assessments (Davies et al., 2006; Stoddard et al., \u003cspan citationid=\"CR56\" class=\"CitationRef\"\u003e2006\u003c/span\u003e). In addition to measuring deviations in biological condition from a minimally disturbed condition, non-linear attributes of metrics can be used to develop consensus for management targets (Muradian, \u003cspan citationid=\"CR29\" class=\"CitationRef\"\u003e2001\u003c/span\u003e). Tipping points in relationships showing assimilative capacity can provide justification for protecting ecosystems at pollution levels below those tipping points and thus preventing major changes in ecosystems (Soranno et al., \u003cspan citationid=\"CR46\" class=\"CitationRef\"\u003e2008\u003c/span\u003e; Stevenson et al., \u003cspan citationid=\"CR52\" class=\"CitationRef\"\u003e2008\u003c/span\u003e; Stevenson, \u003cspan citationid=\"CR48\" class=\"CitationRef\"\u003e2011\u003c/span\u003e), such as threshold changes in calcareous periphyton abundance in the Everglades (Stevenson et al., \u003cspan citationid=\"CR53\" class=\"CitationRef\"\u003e2002\u003c/span\u003e; Stevenson, \u003cspan citationid=\"CR49\" class=\"CitationRef\"\u003e2014\u003c/span\u003e).\u003c/p\u003e\u003cp\u003eIn addition to diversity traits and weighted average models that infer specific pollution levels, measures of changes in proportion to sensitive native taxa and pollution tolerant taxa align with standard attributes of biological condition (Davies et al., 2006). Some studies show that a metric characterizing the proportion of the taxa that are sensitive or tolerant to pollution can perform as well as or better than metrics that include weights for species relative abundances in samples (Wang et al., \u003cspan citationid=\"CR65\" class=\"CitationRef\"\u003e2005\u003c/span\u003e; Stevenson et al., \u003cspan citationid=\"CR52\" class=\"CitationRef\"\u003e2008\u003c/span\u003e; Charles et al., \u003cspan citationid=\"CR8\" class=\"CitationRef\"\u003e2021\u003c/span\u003e). The high performance of metrics that measure changes in the proportions of taxa with traits seems beneficial because they are more likely directly related to measuring losses in sensitive taxa or gains in tolerant taxa, even though there are problems with accuracy of these metrics measuring losses or gains in taxa when observing such small proportions of the organisms present. However, weighting relative abundances of taxa in samples should provide additional information that would make metrics perform better than simpler metrics that only account for presence and absence of taxa. Perhaps variability in relative abundance of abundant species is so great that it reduces precision of metrics that include relative abundances. If that is the case, down-weighting abundances with log-transformations could improve metric performance.\u003c/p\u003e\u003cp\u003eThe specific goals of this paper are to characterize the effect of human disturbance on biological condition of periphyton in the BCNP and to develop and test metrics for a periphyton monitoring program. First, we characterize changes in biological condition with ordination, cluster analysis, and metric changes along a human disturbance gradient. During metric evaluation we compare metrics with taxa phosphorus traits from the literature with revised traits from BCNP. We also test a novel metric calculated using log-transformed relative abundances to reduce negative effects that high variability in abundant taxa could have on metric precision. We also determine whether natural variability among habitats sampled affects metrics, and we characterized metric change along the phosphorus gradient to inform development of effects-based biological and pollution management targets. Using those analyses, we evaluate metric performance to aid selection of a subset of tested metrics for use in future monitoring. Finally, we explore relationships between metrics and the major pollutant, indicated by mat P, to provide benchmarks in the mat P gradient for potential management targets.\u003c/p\u003e"},{"header":"Methods","content":"\u003cp\u003e\u003c/p\u003e\u003cp\u003eStudy area\u003c/p\u003e\u003cp\u003eBig Cypress National Preserve is located southwest of Lake Okeechobee and borders the northern edge of Everglades National Park (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003e). It spans 295,016 ha and receives approximately 425,000 visitors annually who come to hike, canoe, camp, bird-watch, hunt, fish, and use off-road vehicle (ORV) trails (Patterson et al., \u003cspan citationid=\"CR35\" class=\"CitationRef\"\u003e2008\u003c/span\u003e). The preserve contains a large remnant of natural wetland mosaic including cypress strands and domes, pine forests, wet prairies, marshes, sloughs, mangrove forests, and hardwood hammocks. The preserve also contains large stands of dwarf cypress, as well as rare orchids, bromeliads, and ferns (Patterson et al. \u003cspan citationid=\"CR35\" class=\"CitationRef\"\u003e2008\u003c/span\u003e).\u003c/p\u003e\u003cp\u003eBig Cypress is fundamentally different than its adjacent Everglades neighbor. The soil is comprised of sands (Pamlico Sands) overlying caprock (hardened limestone) of the Tamiami Formation whereas the vast majority of the greater Everglades system has an organic peat substrate overlying Miami Limestone (Jarosewich \u0026amp; Wagner, \u003cspan citationid=\"CR19\" class=\"CitationRef\"\u003e1985\u003c/span\u003e). Water transmissivity is more localized and slower in Big Cypress than the Everglades (Jarosewich \u0026amp; Wagner, \u003cspan citationid=\"CR19\" class=\"CitationRef\"\u003e1985\u003c/span\u003e). Concentration of surface water total phosphorus is typically higher in Big Cypress compared to the greater Everglades. This has been attributed to higher total phosphorus from natural sources such as shallow soils, rocks, and ground water (Miller et al., \u003cspan citationid=\"CR28\" class=\"CitationRef\"\u003e2004\u003c/span\u003e).\u003c/p\u003e\u003cp\u003eThe name \u0026ldquo;Big Cypress\u0026rdquo; refers to the vast expanse of cypress (\u003cem\u003eTaxodium distichum\u003c/em\u003e var. \u003cem\u003eimbricarium\u003c/em\u003e (Nuttall) Croom) rather than to the size of the trees. The larger bald cypress trees were logged during the past two centuries. Extraction of oil, gas, and minerals occurs within BCNP and surrounding areas to the north (Whelan et al., \u003cspan citationid=\"CR66\" class=\"CitationRef\"\u003e2020\u003c/span\u003e). North of BCNP and south of Lake Okeechobee, agriculture fields have been operating since the early 1900s. To establish the farming industry, the lands were first drained, and the natural water flow was redirected. The western region of Big Cypress has headwaters that begin north of BCNP and are conveyed south via a canal system (Barron River Canal, Miller et al., \u003cspan citationid=\"CR28\" class=\"CitationRef\"\u003e2004\u003c/span\u003e). Since BCNP establishment in 1974, the city of Immokalee\u0026rsquo;s population has grown 6 times greater (just north of the Preserve) (US Census data 1980, 2020). In addition, farming and other industrial development has expanded north of the BCNP boundary. The primary management activities have been related to restoration of the human altered regional hydrology to a more natural flow pattern, improving water quality entering the park from the north, manage invasive species, balancing recreational and extractive uses with long-term sustainability of the system, as well as protecting and preserving natural resources that include rare species (Urgelles et al., \u003cspan citationid=\"CR60\" class=\"CitationRef\"\u003e2019\u003c/span\u003e).\u003c/p\u003e\u003cp\u003eSampling design\u003c/p\u003e\u003cp\u003ePeriphyton sampling has been conducted since 2008 to investigate potential and test methods for using periphyton to assess ecological condition in BCNP. The target sampling population was periphyton within mapped graminoid and broadleaf marshes in northwestern BCNP. The sampling area has been restricted to mapped and accessible marshes (which includes access by helicopter, ORV, truck, or hiking) in the northwest section of BCNP because that section has regions close to areas of human disturbance with elevated phosphorus concentrations as well as minimally disturbed conditions regions in the southern and eastern parts of the sampling area. Sample units were contiguous graminoid and broadleaf marsh habitats within a 250 m radius of selected map grid cell centroids from the Western Big Cypress National Preserve Vegetation Map (Whelan et al., \u003cspan citationid=\"CR66\" class=\"CitationRef\"\u003e2020\u003c/span\u003e).\u003c/p\u003e\u003cp\u003eThe northwestern BCNP is divided into basins separated by artificial or natural structures that may or may not impede water flow during part or all the year (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003e). Each basin is treated as a separate sampling strata or block. Using a restricted stratified random design, a number of potential sites are selected within the desirable habitat of each basin. Each potential site was evaluated and those that meet specific criteria were sampled annually (Londo\u0026ntilde;o, \u003cspan citationid=\"CR22\" class=\"CitationRef\"\u003e2019\u003c/span\u003e). The goal is to sample the same site year after year (Urgelles et al., \u003cspan citationid=\"CR60\" class=\"CitationRef\"\u003e2019\u003c/span\u003e).\u003c/p\u003e\u003cp\u003eFor hydrological sample years 2013 and 2014, sampling followed very similar criteria as above except some sites were haphazardly selected to allow for: (a) maximizing spatial spread within a basin, (b) co-location with water-quality monitoring stations where available, (c) answering specific monitoring questions and (d) ensuring half the sites were accessible via helicopter and the other half of the sites accessible via ORV trails. Additionally, sampling occurred in the Kissimmee Billy Strand basin, which is an eastern basin that was only sampled in 2013 and 2014.\u003c/p\u003e\u003cp\u003eField collection methods\u003c/p\u003e\u003cp\u003eSampling goals and methods varied during early years to explore potential for periphyton monitoring. Starting in 2013, mat phosphorus (mat P) concentrations were measured as well as periphyton species composition. Mat P concentrations vary less than water column P as indicators of P pollution and have been recommended for measurement of P in the Everglades (Gaiser et al., \u003cspan citationid=\"CR12\" class=\"CitationRef\"\u003e2004\u003c/span\u003e). Therefore, periphyton metric development was limited to hydrologic years after 2013.\u003c/p\u003e\u003cp\u003eFor hydroyears 2013, 2014, 2017, and 2019, a general overview of methods upon site arrival was to locate where water was present and collect a floating periphyton sample if possible. Attempts were made to access the same depressional marsh for revisits. Sampling followed a hydrological clock, approximately two months after the peak of the wet season, whereby sites were accessed at a similar hydroperiod stage year after year (typically November or December). Thus, water was usually present at all sites during collection. However, if the site was dry, then dry periphyton was collected.\u003c/p\u003e\u003cp\u003eSeparate grab samples were collected for diatom and mat P analyses, with each composed of a minimum of five grabs within a 5-meter radius sampling area to get two 125 mL samples. The preferred order of substrate collection was (1) floating mat, (2) algae on plants (sweaters), (3) algae on sediments (benthic mats), and (4) algae on woody debris. Grab samples are placed inside two 125-milliliter Nalgene\u003csup\u003e\u0026reg;\u003c/sup\u003e opaque bottles with as little water as possible. Excess water was decanted prior to fixation. Additional site characterization information was collected including vegetation data, pH, water conductivity, and water depth. Samples were preserved as quickly as possible upon return from the field. Periphyton mat P samples were chilled in an ice slurry in the field and frozen at the lab. The periphyton composition samples were fixed with a 3% buffered formalin solution (Muxo \u0026amp; Shamblin, \u003cspan citationid=\"CR30\" class=\"CitationRef\"\u003e2019\u003c/span\u003e).\u003c/p\u003e\u003cp\u003eFor hydroyears 2020, the only modification was the periphyton grab samples did not have a preferred substrate collection order but instead, periphyton was collected in relative proportion of the substrates present at the site, e.g. as floating mat, sweaters, benthic map or as debris. Otherwise, sampling was the same as above.\u003c/p\u003e\u003cp\u003eSample processing and analysis\u003c/p\u003e\u003cp\u003ePeriphyton samples were sent to Florida International University for the assay of mat P, which involved digesting mats and assaying total phosphorus (TP). In the laboratory, periphyton sample wet weights were recorded, thawed, and then extraneous plant material and animals were manually separated from the periphyton and placed in foil packets. These were placed in an 80\u0026deg;C oven for three days; once dry, the weight was recorded as \u0026ldquo;extraneous material\u0026rdquo;, and then discarded. The TP content of periphyton was expressed on a dry-weight basis because the organically incorporated P was not separable from that bound to calcite. The TP subsample was dried at 80\u0026deg;C and ground down to a fine powder with a mortar and pestle. Colorimetric analysis was used to estimate TP concentrations of the periphyton subsamples following the methods of Solorzano and Sharp (\u003cspan citationid=\"CR45\" class=\"CitationRef\"\u003e1980\u003c/span\u003e), which were then used to calculate \u0026micro;g g\u003csup\u003e\u0026minus;\u0026thinsp;1\u003c/sup\u003e dried periphyton mat (not ash free dry mass). For hydroyears 2013, 2014, and 2017 the lab methods followed Smith (\u003cspan citationid=\"CR44\" class=\"CitationRef\"\u003e2018\u003c/span\u003e) and for hydrologic year 2019 and 2020 they followed Wilson (\u003cspan citationid=\"CR67\" class=\"CitationRef\"\u003e2020\u003c/span\u003e). The difference in methods is due to the acquisition of a newer spectrophotometer.\u003c/p\u003e\u003cp\u003eEarly, unpublished analyses of data indicated diatom species composition was related to the human disturbance gradient in BCNP as well as composition of both diatom and non-diatom algae. Afterward algal analyses were limited to diatom species composition for development of diatom metrics to reduce costs of periphyton analyses. Periphyton samples were sent to Michigan State University to be assayed for diatom species composition. There, samples were acid cleaned in nitric acid and mounted on microscope slides in NAPHRAX\u0026copy;. At least 600 diatom valves were identified and counted for each slide.\u003c/p\u003e\u003cp\u003eData analysis: data selection and descriptive analyses\u003c/p\u003e\u003cp\u003eWe created independent calibration and validation datasets with only one sample per site in each dataset to avoid problems with pseudoreplication and repeated measures. Both subsets of samples were created using samples from the 2013\u0026ndash;2020 period when mat P was assayed and by randomly selecting one sample from each site without replacement of selected sites. This created a calibration dataset with 113 samples for characterizing changes in diatom assemblages along the P gradient and for characterizing species traits. 78 of the 113 calibration sites had all chemistry information. 105 sites had all chemistries except pH. A second set of samples from the 2013\u0026ndash;2020 period with chemistries was selected for a validation dataset to test metrics, which resulted in 85 randomly selected samples with one sample per site that was not selected for the calibration dataset.\u003c/p\u003e\u003cp\u003eOrdination analyses indicated that environmental factors, sampling year, and taxonomist were related to variability in diatom species composition. Follow-up analyses of the data indicated that residuals in diatom metric relationships with mat P concentration also differed among sampling years. Different taxonomists were assigned to count diatom samples from year to year with turnover of taxonomists present in the lab. Interannual residual variation could not be assigned to either taxonomist or environment with certainty, so we restricted metric calculation and testing to hydrologic years 2017, 2019, and 2020 when the same taxonomist counted all samples. That reduced the validation dataset for testing metrics to the 44 sites sampled annually in the 7 regions with one sample per site for either 2017, 2019, or 2020. However, we used the full calibration dataset from 2013\u0026ndash;2020 for the preliminary descriptive analyses: ordination, cluster analysis, trait characterization, and characterization of minimally disturbed condition. Our rationale was despite the concern about potential effect of taxonomist or interannual environmental differences, we wanted to ensure sufficient sample size to detect species-environment relationships and characterize traits for as many taxa as possible so future results would be resilient to interannual variation caused either by environment or taxonomist.\u003c/p\u003e\u003cp\u003eSpecies-environment relationships were evaluated with ordination using the calibration dataset. Non-metric multidimensional scaling (NMDS) was selected to ordinate samples in species space and then relate environmental variables to the NMDS axes with the \u003cem\u003evegan\u003c/em\u003e package in R (Oksanen et al., \u003cspan citationid=\"CR31\" class=\"CitationRef\"\u003e2020\u003c/span\u003e). Cluster analyses using Bray-Curtis dissimilarity and the calibration dataset were calculated with the \u003cem\u003eclust\u003c/em\u003e command in R to observe patterns in species composition without constraints of environmental factors. Groups of sites with low dissimilarity were identified and then differences in mat P, conductivity, and pH among these site groups were determined.\u003c/p\u003e\u003cp\u003eTo more thoroughly understand and illustrate the succession in diatom assemblages along the P gradient, we used a stacked bar chart using the most abundant low and high P indicator species. Relative rates of change in species relative abundances and the dominance of these taxa in assemblages along the P gradient were determined for the calibration dataset. The most abundant low P and high P taxa were selected for the figure. In addition, we created a heat map of all taxa in 5 samples or more using an Excel\u0026copy; spreadsheet. The heat map matrix had taxa abundances colored in rows and samples in columns, with rows ordered by TWINSPAN groups (R package) and columns ordered by increasing mat P concentrations. Cells of the heat map matrix were colored yellow, orange, or red to indicate increasing relative abundances of taxa.\u003c/p\u003e\u003cp\u003eData analysis: taxa trait characterization\u003c/p\u003e\u003cp\u003eTaxa have been characterized by the NPS as being characteristic of oligotrophic or eutrophic habitats based on consensus and best professional interpretation of information reported in literature sources (Table \u003cspan refid=\"MOESM1\" class=\"InternalRef\"\u003eS1\u003c/span\u003e). This literature review assessed 13 papers and attempted to determine the optimum TP, trophic habitat preference, pH preference, and trophic indication status of diatoms at the species level.\u003c/p\u003e\u003cp\u003eTraits analysis was conducted to characterize taxa as low and high P taxa using the calibration dataset. Multiple trait characterization methods were used to be thorough and to evaluate consistency in characterizations. Linear and polynomial regression were used to determine whether relationships between taxa relative abundances and mat P concentration were negative or positive to identify low and high P taxa, respectively. Polynomial regression was used to complement linear regression to ensure non-linear responses were detected. TITAN (Threshold Indicator Taxa Analysis (Baker \u0026amp; King \u003cspan citationid=\"CR1\" class=\"CitationRef\"\u003e2010\u003c/span\u003e)) was also used to describe taxa with decreasing and increasing indicator species values along the TP gradient for low and high P taxa, respectively. Finally, weighted average TP optima were calculated for taxa. Low and high P characterizations (i.e. taxa traits) were only calculated for taxa observed in 5 or more samples from the calibration set of 78 samples from different sites and having P concentrations measured because TITAN requires 5 observations per taxon to provide reasonable statistical power for characterizing taxa traits with regression. Mat P optima and tolerances were determined for all taxa because our experience showed weighted average models (WAM) were most precise when all taxa are used in the model. The NPS does not plan to use the WAM model in future monitoring. The WAM model is used in this paper to provide a benchmark for a model that is usually the most precisely related to an environmental gradient (Reavie et al., \u003cspan citationid=\"CR41\" class=\"CitationRef\"\u003e2008\u003c/span\u003e).\u003c/p\u003e\u003cp\u003eTaxa traits determined with the 78-sample calibration dataset were compared to taxa traits observed in past Everglades research to compare consistency in results and evaluate causal relationships between low and high P traits and P concentrations. Pairwise relationships between taxa trait characterizations in this study and those by Hagerthey et al. (\u003cspan citationid=\"CR17\" class=\"CitationRef\"\u003e2012\u003c/span\u003e) and Gaiser et al. (\u003cspan citationid=\"CR11\" class=\"CitationRef\"\u003e2006\u003c/span\u003e) were determined. Species characterizations for BCNP were also compared to low and high P characterizations by Slate and Stevenson (\u003cspan citationid=\"CR43\" class=\"CitationRef\"\u003e2007\u003c/span\u003e), which were determined by experimental manipulation of P with \u003cem\u003ein situ\u003c/em\u003e mesocosms. The experiments enabled establishing a causal relationship between P concentration and relative abundances of diatom taxa (Beyers \u003cspan citationid=\"CR2\" class=\"CitationRef\"\u003e1998\u003c/span\u003e).\u003c/p\u003e\u003cp\u003eData analysis: metric evaluation\u003c/p\u003e\u003cp\u003eWe then evaluated metric design on their performance to be sensitive to environmental change, to characterize changes in biological condition along the human disturbance gradient, and to help justify management targets for P concentration. We did not expect the same metric to be best for all goals (Fig. \u003cspan refid=\"MOESM1\" class=\"InternalRef\"\u003eS1\u003c/span\u003e, Stevenson, \u003cspan citationid=\"CR47\" class=\"CitationRef\"\u003e2006\u003c/span\u003e, Stevenson, \u003cspan citationid=\"CR48\" class=\"CitationRef\"\u003e2011\u003c/span\u003e). We sought a metric with a linear response to P concentration to be sensitive to environmental change, because metrics with linear response have the same sensitivity to incremental changes in P concentration at all levels of P concentration. We sought non-linear responses of metrics of biological condition to detect thresholds in biological responses that could be used to justify P management targets. We did not calculate a multimetric index because that was not the plan for BCNP application.\u003c/p\u003e\u003cp\u003eFor each sample we calculated multiple metrics with each trait, plus three diversity metrics: taxa number observed in standardized counts, Shannon diversity, and Pielou\u0026rsquo;s evenness. Counts were standardized to 600 valves, because some counts had more than 600 valves. Counts were standardized by assigning random numbers to valves in counts and picking valves with the lowest 600 numbers. For each trait, four metric types were calculated: 1) proportion of individuals with a trait (PropValves) as sum(v\u003csub\u003eijt\u003c/sub\u003e/V\u003csub\u003ej\u003c/sub\u003e), where v\u003csub\u003eijt\u003c/sub\u003e is the number of valves counted for taxon \u003cem\u003ei\u003c/em\u003e with trait \u003cem\u003et\u003c/em\u003e (e.g. either high or low P using a specific trait determination method) in sample \u003cem\u003ej\u003c/em\u003e, and V\u003csub\u003ej\u003c/sub\u003e is the number of valves of all taxa with assigned traits in sample \u003cem\u003ej\u003c/em\u003e; 2) the proportion of taxa with a trait (PropTaxa), as t\u003csub\u003etj\u003c/sub\u003e/T\u003csub\u003ej\u003c/sub\u003e where t\u003csub\u003etj\u003c/sub\u003e is the number of taxa with trait \u003cem\u003et\u003c/em\u003e in sample \u003cem\u003ej\u003c/em\u003e and T\u003csub\u003ej\u003c/sub\u003e is the number of taxa with traits assigned in sample \u003cem\u003ej\u003c/em\u003e; and 3) the number of taxa in a sample with a specific trait (noTaxa, e.g. number of low P taxa). The fourth metric type was RlogA and is described in the next paragraph. Below we also describe our use of the three trait determination methods: review of the literature, TITAN, and regression.\u003c/p\u003e\u003cp\u003eRecent papers have shown metrics calculated as proportion of taxa with a trait are often related better to human disturbance and stressor gradients than metrics calculated as proportion of individuals with a trait (Stevenson et al., \u003cspan citationid=\"CR52\" class=\"CitationRef\"\u003e2008\u003c/span\u003eb; Carlisle et al., \u003cspan citationid=\"CR7\" class=\"CitationRef\"\u003e2022\u003c/span\u003e). Because relative performance of species in a habitat likely varies with environmental conditions as well as their presence in a habitat, we tried down-weighting metrics by log transforming relative abundances of taxa and thereby reducing variability caused by slight changes in growth rates of abundant taxa that could produce great differences in a taxon\u0026rsquo;s abundance. We determined relative log abundance (RlogA\u003csub\u003eijt\u003c/sub\u003e) by calculating the natural log of valve numbers of taxon \u003cem\u003ei\u003c/em\u003e with trait \u003cem\u003et\u003c/em\u003e in sample \u003cem\u003ej\u003c/em\u003e and then the sum of natural log transformed valve numbers for all taxa with traits assigned in sample \u003cem\u003ej\u003c/em\u003e as sum(RlogA\u003csub\u003eijt=h+l\u003c/sub\u003e), summed for all taxa with both high or low P value (t\u0026thinsp;=\u0026thinsp;h\u0026thinsp;+\u0026thinsp;l) assigned with a specific method (literature, TITAN, regression). The RlogA metric values for either high or low P RlogA were calculated as the sum of all RlogA\u003csub\u003eijt\u003c/sub\u003e for taxa with either low or high P traits based on a specific trait determination divided by the sum of all RlogA\u003csub\u003eijt\u003c/sub\u003e for taxa with both low and high P traits, i.e. RlogA\u003csub\u003ej\u003c/sub\u003e= sum(RlogA\u003csub\u003eijt\u003c/sub\u003e)/SRlogA\u003csub\u003ej\u003c/sub\u003e.\u003c/p\u003e\u003cp\u003eMetrics were evaluated by comparing their adjusted \u003cem\u003eR\u003c/em\u003e\u003csup\u003e\u003cem\u003e2\u003c/em\u003e\u003c/sup\u003e values determined for slopes in relationships between metrics and log transformed mat P. We log\u003csub\u003e2\u003c/sub\u003e transformed mat P to even the distribution of the observations along the mat P gradient, because the number of low mat P samples was much greater than the number of high mat P samples. Unless described otherwise, a metric-mat P relationship in the following text refers to the relationship between a metric and log transformed mat P. The R package \u003cem\u003elm\u003c/em\u003e was used for linear regression (R Core Team, \u003cspan citationid=\"CR39\" class=\"CitationRef\"\u003e2021\u003c/span\u003e). Both actual values and ranks of \u003cem\u003eR\u003c/em\u003e\u003csup\u003e\u003cem\u003e2\u003c/em\u003e\u003c/sup\u003e values were evaluated. Three analyses of variance (ANOVA) were used to compare average adjusted \u003cem\u003eR\u003c/em\u003e\u003csup\u003e\u003cem\u003e2\u003c/em\u003e\u003c/sup\u003e for metric-mat P relationships for: 1) different metric types (noTaxa, PropTaxa, PropValves, RlogA); 2) trait source (literature, regression, or TITAN);, or 3) trait (low P versus high P).\u003c/p\u003e\u003cp\u003eTo control for other metric attributes (metric type, trait source, and trait) when evaluating performance of one of the metric attributes, Kruskal tests (\u003cem\u003ekruskal.test\u003c/em\u003e in R) were run to compare ranks of adjusted \u003cem\u003eR\u003c/em\u003e\u003csup\u003e\u003cem\u003e2\u003c/em\u003e\u003c/sup\u003e of metric-mat P relationships. Because these analyses are rather complicated to set up, we describe them in detail. Three Kruskal tests were run to evaluate metric attributes, and each used the same set of 24 metric-mat P relationships resulting from the 24 combinations of 4 metric types, 3 trait sources, and 2 traits. To compare performance of metric types (noTaxa, PropTaxa, PropValves, RlogA) while controlling for trait source and trait, we grouped metric-mat P relationships into the 6 different trait source-trait groups (low P-Lit, low P-Regr, low P-TITAN, high P-Lit, high P-Regr, high P-TITAN). Then we ranked (1\u0026ndash;4) the adjusted \u003cem\u003eR\u003c/em\u003e\u003csup\u003e\u003cem\u003e2\u003c/em\u003e\u003c/sup\u003e of the metric-mat P relationships for the 4 metric types within each of the 6 trait source-trait (TS-T) groups. In the case for the Kruskal test comparing metric types, ranks ranged from 1\u0026ndash;4 because there were 4 metric-mat P relationships having different metric types within each TS-T group. In the Kruskal test for metric type, each metric type had 6 rankings because there were 6 TS-T groups. Using the same approach we compared performance of metrics with either Lit, Regr, or TITAN trait sources while holding metric type and trait constant; we used the Kruskal test to compare ranks (1\u0026ndash;3 trait sources) of adjusted \u003cem\u003eR\u003c/em\u003e\u003csup\u003e\u003cem\u003e2\u003c/em\u003e\u003c/sup\u003e for the 3 metric-mat P relationships in the 8 metric type-trait source (MT-TS) groups (noTaxa-Lit, PropT-Lit, PropV-Lit, RLogV-Lit, noTaxa-Regr, PropT-Regr, PropV-Regr, RLogV-Regr, noTaxa-TITAN, PropT-TITAN, PropV-TITAN, RLogV-TITAN). To compare performance of metrics with either low P or high P traits while holding metric type and trait source constant, we used the Kruskal test to compare ranks (1 or 2) of adjusted \u003cem\u003eR\u003c/em\u003e\u003csup\u003e\u003cem\u003e2\u003c/em\u003e\u003c/sup\u003e for those 2 metric-mat P relationships in the 12 trait source\u0026ndash; metric type groups (TS-MT groups: Lit-noTaxa, Lit-PropT, Lit-PropV, Lit-RLogV, Regr-noTaxa, Regr-PropT, Regr-PropV, Regr-RLogV, TITAN-noTaxa, TITAN-PropT, TITAN-PropV, TITAN-RLogV).\u003c/p\u003e\u003cp\u003eMat P and metrics were related to latitude (lat) and longitude (long) to determine their relationship to distance from the main source of human disturbance in the northwestern corner of the BCNP study region. Linear regression (\u003cem\u003elm\u003c/em\u003e in R) was used to relate mat P and metrics to latitude, longitude, and a lat-long interaction term.\u003c/p\u003e\u003cp\u003eDifferences in metrics among regions were determined with the 2020 data to evaluate their performance and demonstrate application by determining regional differences that are related to levels of human disturbance. The regional approach was chosen because the NPS plans to analyze and report results by region using sites as replicates as well as to determine differences among regions and changes in time by region. Single factor ANOVA were used to compare average metric values among regions. \u003cem\u003ep\u003c/em\u003e values resulting from these analyses were used to evaluate metric performance with the average \u003cem\u003ep\u003c/em\u003e rank for metric types (noTaxa, PropTaxa, PropValves, and RlogA), low and high P traits, and trait sources (literature, regression, and TITAN).\u003c/p\u003e\u003cp\u003eNon-parametric Kruskal tests were again used, as used to compare metric-mat P relationships above, to determine which metric attributes (metric type, trait source, trait) were most important for metric performance when comparing regions. This set of three separate Kruskal tests for metric type, trait source, and trait used ranks of 24 ANOVA F-values for comparing differences in the 24 possible metrics among regions. Kruskal tests used ranks of F-value ranks ranging from 1\u0026ndash;4 to compare metric types within the 6 trait source-trait groups, ranging from 1\u0026ndash;3 to compare trait sources within the 8 metric type-trait groups, and ranging from 1\u0026ndash;2 to compare the two traits within the 12 trait source-metric type groups. In addition, we calculated an ANOVA for mat P differences among regions to compare F-values with diatom metrics to determine if diatom metrics more precisely differed among regions than mat P.\u003c/p\u003e\u003cp\u003eData analysis: metric correction for natural factors\u003c/p\u003e\u003cp\u003eResiduals in metric-mat P relationships were related to naturally varying ecological factors to determine if expectations for metrics should be adjusted for habitat type, substrate location, substrate type, and mean water depth. Habitat type included broadleaf marsh, cypress dome, cypress scrub, graminoid marsh and mixed combinations of these habitats. Substrate locations were either floating on the water surface, on benthic sediment or soils (ground), or enveloping macrophytes (sweaters). Substrate types were filamentous algae, periphyton, and soil. Water depth was measured at the location and could covary with the major human disturbance gradient originating with human alterations in the northeast corner of the sampling area. Univariate ANOVA were used to determine separately the effects of differing habitat types, substrate locations, and substrate types on metric-mat P relationships. Linear regression in R (\u003cem\u003elm\u003c/em\u003e) was used to characterize the relationship between water depth and residuals of the metric-mat P relationship. \u003cem\u003ep\u003c/em\u003e values for statistical significance were reported in results without accounting for multiple tests, but multiple tests were accounted for by dividing reported \u003cem\u003ep\u003c/em\u003e values by the number of tests performed when interpreting likelihood that observed results were not due to random patterns in the data. This is like a Bonferroni correction for multiple tests. In addition, we looked at the proportion of \u003cem\u003ep\u003c/em\u003e values that were less than 0.05.\u003c/p\u003e\u003cp\u003eResiduals in the relationships for metrics and mat P as function of latitude and longitude were related to natural environmental factors to determine if natural environmental factors affected both mat P and metrics. This second residual analysis differed from the first analysis, because it addressed whether different habitat conditions had different background P conditions and metric values versus the first analysis which evaluated whether natural environmental factors introduced bias in metric-mat P relationships.\u003c/p\u003e\u003cp\u003eData analysis: management benchmarks\u003c/p\u003e\u003cp\u003eWe used multiple lines of evidence to establish benchmarks for assessment and possible management targets for mat P and periphyton metrics. Here we use the term benchmark to mean a level of an environmental variable that could be used by resource managers to set management targets. To aid establishment of management targets for mat P, we evaluated non-linearities in metric-mat P relationships for tipping points (changes in the linear patterns) at the low-P end of a mat P range, thereby showing natural assimilative capacity for the biological condition attributes measured by a metric (Fig. \u003cspan refid=\"MOESM1\" class=\"InternalRef\"\u003eS1\u003c/span\u003e).\u003c/p\u003e\u003cp\u003eWe used three approaches for determining benchmarks in biological response along the mat P gradient. Decision tree analysis was calculated relating the 29 trait-based, diversity, and WAM metrics to mat P using the \u003cem\u003erpart\u003c/em\u003e package in R (Therneau \u0026amp; Atkinson, \u003cspan citationid=\"CR59\" class=\"CitationRef\"\u003e2025\u003c/span\u003e) to determine change points along the mat P gradient in trees with the lowest error variance. We used the 75th percentile of mat P concentrations for low P groups of samples identified with cluster analysis as another benchmark. We also used the point with highest TITAN z scores for low P taxa as a management benchmark (Baker et al. 2010). We compared benchmarks to changes in species composition illustrated in the heat map of taxon abundances in samples with successively higher mat P to interpret biological responses associated with these mat P benchmarks.\u003c/p\u003e\u003cp\u003eWe related mat P to latitude and longitude to determine if it was related to distance from the northwest corner of the BCNP study area where human disturbance was highest. To ensure management targets based on effects of mat P did not under or over-protect BCNP wetlands, we determined likely mat P concentrations at sites we assumed were minimally disturbed sites. We selected all sites in the Fire Prairie (FP), Monument (MN), and Little Marsh (LM) regions because they were the regions furthest south and east in the BCNP study area. However, some sites in these regions were close to canals and trails which could introduce disturbance.\u003c/p\u003e"},{"header":"Results","content":"\u003cp\u003eChanges in diatom assemblages along the P gradient\u003c/p\u003e\u003cp\u003eNMDS with all samples with all chemistries in the calibration dataset (n\u0026thinsp;=\u0026thinsp;78) showed diatom assemblages responded most to mat P. The NMDS had an insignificant deviation from 1 in the stress plot. The continuous variables mat P, year sampled, pH, water temperature, and conductivity as well as the categorical variables habitat, water color, and taxonomist were significantly related to NMDS axes (\u003cem\u003ep\u003c/em\u003e\u0026thinsp;\u0026le;\u0026thinsp;0.041, Fig. S2). Mat P, sampling year, and habitat were the most highly correlated abiotic variables with 0.624, 0.349, and 0.273 \u003cem\u003eR\u003c/em\u003e\u003csup\u003e2\u003c/sup\u003e values, respectively (\u003cem\u003ep\u003c/em\u003e\u0026thinsp;\u0026lt;\u0026thinsp;0.001, Table S2). pH and water temperature were next most closely related to NMDS axes, with 0.162 and 0.128 \u003cem\u003eR\u003c/em\u003e\u003csup\u003e\u003cem\u003e2\u003c/em\u003e\u003c/sup\u003e values respectively. Mat P and pH were inversely correlated to each other (Pearson \u003cem\u003er\u003c/em\u003e=-0.357) along the NMDS axes.\u003c/p\u003e\u003cp\u003eCluster analysis of all 113 samples at unique sites in the calibration dataset isolated 5 groups of sites at a dissimilarity level of 0.75 (Fig. S3). From low to high along the mat P scale, these groups were designated 1 to 5. Group 1 had more than two thirds of the sites. Groups 2 to 5 had 5, 1, 9 and 10 sites respectively. Group 3 was dropped from the following discussion because it only had one sample. Three subgroups in Group 1 that had dissimilarity less than 0.6 among sites were identified and designated as Groups 1\u0026ndash;1, 1\u0026ndash;2, and 1\u0026ndash;3 in order of dissimilarity among groups. Group 1\u0026ndash;1 had very low dissimilarity among sites and Group 1\u0026ndash;3 had the highest dissimilarity of the Group 1 subgroups. Dissimilarity increased significantly (Tukey HSD, adj \u003cem\u003ep\u003c/em\u003e\u0026thinsp;\u0026lt;\u0026thinsp;0.001, Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003e) from a mean of 0.39 among samples in Group 1\u0026ndash;1 to 0.46 and 0.47 in Groups 1\u0026ndash;2 and 1\u0026ndash;3, which were themselves not statistically different. Dissimilarity also increased significantly (Tukey HSD, adj \u003cem\u003ep\u003c/em\u003e\u0026thinsp;\u0026lt;\u0026thinsp;0.001) from Groups 1\u0026ndash;2 and 1\u0026ndash;3 to Groups 2, 4, and 5, which had 0.54 average dissimilarity for those three groups of sites.\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003cp\u003eMat P, pH, and conductivity, varied significantly among site groups (ANOVA, \u003cem\u003ep\u003c/em\u003e\u0026thinsp;\u0026lt;\u0026thinsp;0.05, Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003e). With 75th and 90th percentiles of 142 and 201 \u0026micro;g P g\u003csup\u003e\u0026minus;\u0026thinsp;1\u003c/sup\u003e dry mass, mat P was lower in Group 1\u0026ndash;1 than all other site groups (Tukey HSD, adj \u003cem\u003ep\u003c/em\u003e\u0026thinsp;\u0026lt;\u0026thinsp;0.001). Average mat P for sites in either Group 1\u0026ndash;2 or 1\u0026ndash;3 was lower than for sites in either Group 2, 4, or 5. The 75th and 90th percentiles of all sites in Group 1 were 354 and 508 \u0026micro;g g\u003csup\u003e\u0026minus;\u0026thinsp;1\u003c/sup\u003e mat P. pH of Group 1\u0026ndash;1 sites was greater than sites in either Group 2, 4, or 5 (Tukey HSD, adj \u003cem\u003ep\u003c/em\u003e\u0026thinsp;=\u0026thinsp;0.03\u0026ndash;0.6). Average conductivity for sites in Group 4 was less than in Groups 1\u0026ndash;3 and 5 (Tukey HSD, adj \u003cem\u003ep\u003c/em\u003e\u0026thinsp;=\u0026thinsp;0.04). Tukey HSD comparisons of mean chemistry values for pairwise comparisons of other site groups were not significantly different.\u003c/p\u003e\u003cp\u003eThe number of diatom taxa observed in 600 valve counts increased (\u003cem\u003ep\u003c/em\u003e\u0026thinsp;\u0026lt;\u0026thinsp;0.001) from about 10 in low P to 25 in high P as high P taxa invaded and low P taxa became very rare or were lost (Figs.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003e, \u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003e, and S4). \u003cem\u003eEncyonema evergladianum\u003c/em\u003e, \u003cem\u003eMastogloia calcarea, Encyonopsis microcephala\u003c/em\u003e, and \u003cem\u003eBrachysira ocalanensis\u003c/em\u003e, assigned low P traits in calculations described later in these results, comprised approximately 65 percent of assemblages in low mat P conditions. Their combined proportion of assemblages decreased to 45, 20 and 5 percent in samples with mat P ranging, respectively, from 350\u0026ndash;450, 450\u0026ndash;750, and greater than 750 \u0026micro;g g\u003csup\u003e\u0026minus;\u0026thinsp;1\u003c/sup\u003e mat (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003e). The five dominant high P taxa \u003cem\u003e(Encyonema silesiacum, Navicula cryptotenella, Nitzschia amphibia, Gomphonema gracilis, and Gomphonema auritum\u003c/em\u003e) increased from 6 percent when mat P was less than 100 to about 45 percent of samples when mat P was greater than 750 \u0026micro;g g\u003csup\u003e\u0026minus;\u0026thinsp;1\u003c/sup\u003e (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003e), leaving a high diversity of other diatom taxa to comprise the rest of assemblages in high P conditions (Fig. S4). Whereas other common low P taxa maintained relatively high proportions of counts at intermediate concentrations of mat P, proportions of \u003cem\u003eE. evergladianum\u003c/em\u003e decreased more rapidly than other low P taxa (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003e). \u003cem\u003eEncyonema evergladianum\u003c/em\u003e was observed in 56 of 58 samples with mat P less than 300 \u0026micro;g g\u003csup\u003e\u0026minus;\u0026thinsp;1\u003c/sup\u003e and were usually quite common, however it was only observed in 4 of 28 samples in samples with mat P greater than 550 \u0026micro;g g\u003csup\u003e\u0026minus;\u0026thinsp;1\u003c/sup\u003e (Fig. S4). In contrast, \u003cem\u003eN. amphibia\u003c/em\u003e was only observed in 6 of 49 samples with mat P less than 200 \u0026micro;g g\u003csup\u003e\u0026minus;\u0026thinsp;1\u003c/sup\u003e, but it was highly abundant in 14 of 16 samples with mat P greater than 800 \u0026micro;g g\u003csup\u003e\u0026minus;\u0026thinsp;1\u003c/sup\u003e.\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003cp\u003eTaxa traits\u003c/p\u003e\u003cp\u003eLow and high P traits were assigned to 58 taxa with TITAN and to 55 taxa with regression (Table S3). Weighted average mat P optima were assigned to 157 taxa because the criterion of observing taxa in 5 or more samples was not applied for determination of mat P optima. 62 taxa had low and high P traits assigned from the literature by the NPS. Many more taxa were assigned high P than low P traits. The number of taxa assigned low and high P traits were, respectively, 11 and 47 taxa by TITAN, 8 and 47 taxa by regression, and 23 and 39 by the literature review. TITAN and regression did not assign different traits, either low P or high P, to any taxon; but some of the less common taxa were assigned traits by only TITAN or regression. Low and high P trait assignments did differ for one taxon, \u003cem\u003eE. silesiacum\u003c/em\u003e, which was assigned a high P trait by TITAN and regression and a low P trait using literature references. Many of the abundant low and high P taxa determined by TITAN and regression did not have trait assignments from the literature. The averages of P optima were, respectively: 179 and 1001 for low and high P TITAN taxa; 187 and 1220 for low and high P regression taxa; and 394 and 1062 for low and high P literature-defined taxa.\u003c/p\u003e\u003cp\u003e\u003cem\u003eEncyonopsis evergladianum, M. calcarea, E. microcephala\u003c/em\u003e, and \u003cem\u003eB. ocalanensis\u003c/em\u003e were the four most common low P taxa (Table S3, Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003e). The other four taxa classified as low P by TITAN as well as regression were an unknown species of \u003cem\u003eNitzschia\u003c/em\u003e (sp. 1\u003cem\u003e), Fragilaria synegrotesca, Adlafia bryophila\u003c/em\u003e, and \u003cem\u003eNitzschia serpentiraphe\u003c/em\u003e. The most common high P taxa were \u003cem\u003eE. silesiacum, N. cryptotenella, N. amphibia, G. auritum\u003c/em\u003e, and \u003cem\u003eG. gracile\u003c/em\u003e (Table S3). These most common taxa were not necessarily the most limited to high and low P conditions. The relative preference of taxa for low and high P conditions, versus the absolute preference (low or high) was important for metrics only when calculating the WAM diatom inferred mat P. Of the abundant low P taxa, \u003cem\u003eEncyonema evergladianum\u003c/em\u003e had the lowest mat P optimum (112 \u0026micro;g g\u003csup\u003e\u0026minus;\u0026thinsp;1\u003c/sup\u003e, Table S3). \u003cem\u003eNitzschia amphibia\u003c/em\u003e had one of the highest mat P optima of the common high P taxa (1198 \u0026micro;g g\u003csup\u003e\u0026minus;\u0026thinsp;1\u003c/sup\u003e). \u003cem\u003eEncyonema silesiacum\u003c/em\u003e and \u003cem\u003eN. cryptotenella\u003c/em\u003e had relatively low P optima for taxa characterized by TITAN and regression as high P taxa (439 and 358 \u0026micro;g g\u003csup\u003e\u0026minus;\u0026thinsp;1\u003c/sup\u003e respectively).\u003c/p\u003e\u003cp\u003eTraits assigned in this project with BCNP samples were highly correlated with traits assigned to taxa in three other studies (Fig. S5). Taxalists did vary between our BCNP list and lists in Hagerthey et al. (\u003cspan citationid=\"CR17\" class=\"CitationRef\"\u003e2012\u003c/span\u003e), Gaiser et al. (\u003cspan citationid=\"CR11\" class=\"CitationRef\"\u003e2006\u003c/span\u003e), and Slate and Stevenson (\u003cspan citationid=\"CR43\" class=\"CitationRef\"\u003e2007\u003c/span\u003e). We felt confident comparing BCNP traits for 35 taxa with names that were clearly the same with Hagerthey et al., 23 taxa in Gaiser et al., and 33 taxa in Slate and Stevenson. Spearman correlations among trait values in paired taxalists were positive, all highly significant (\u003cem\u003ep\u003c/em\u003e\u0026thinsp;\u0026lt;\u0026thinsp;0.001), and had \u003cem\u003eR\u003c/em\u003e\u003csup\u003e2\u003c/sup\u003e 0.676, 0.737, and 0.783 for Hagerthey et al., Gaiser et al., and Slate and Stevenson, respectively.\u003c/p\u003e\u003cp\u003eMetric testing and evaluation: metric-mat P relationships\u003c/p\u003e\u003cp\u003eAlmost all metrics except the diversity metrics (Shannon H and Pielou\u0026rsquo;s evenness) were related to mat P with high statistical significance and adjusted \u003cem\u003eR\u003c/em\u003e\u003csup\u003e\u003cem\u003e2\u003c/em\u003e\u003c/sup\u003e (Table S4). When comparing metric performance, metrics calculated with the relative log abundances (RlogA) were more precisely related (based on adjusted \u003cem\u003eR\u003c/em\u003e\u003csup\u003e\u003cem\u003e2\u003c/em\u003e\u003c/sup\u003e) to the mat P gradient than all other metrics, including WAM diatom inferred TP (Figs.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e5\u003c/span\u003e and \u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e6\u003c/span\u003e, Table S4). According to ANOVA and Tukey HSD, average adjusted \u003cem\u003eR\u003c/em\u003e\u003csup\u003e\u003cem\u003e2\u003c/em\u003e\u003c/sup\u003e for relationships between RlogA metrics and log\u003csub\u003e2\u003c/sub\u003e(mat P) were significantly higher than the average adjusted \u003cem\u003eR\u003c/em\u003e\u003csup\u003e\u003cem\u003e2\u003c/em\u003e\u003c/sup\u003e for mat P relationships with noTaxa metrics (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e6\u003c/span\u003eA, Tukey HSD, \u003cem\u003ep\u003c/em\u003e\u0026thinsp;=\u0026thinsp;0.016), but other pairwise comparisons among metric types did not differ (ANOVA, n\u0026thinsp;=\u0026thinsp;24, \u003cem\u003ep\u003c/em\u003e\u0026thinsp;=\u0026thinsp;0.017; Tukey HSD, 0.09\u0026thinsp;\u0026lt;\u0026thinsp;\u003cem\u003ep\u003c/em\u003e\u0026thinsp;\u0026lt;\u0026thinsp;0.83, Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e6\u003c/span\u003eA). However, the ranks of RlogA metric performances were greater than noTaxa, PropTaxa, and PropValves metrics (Kruskal test, n\u0026thinsp;=\u0026thinsp;24, \u003cem\u003ep\u003c/em\u003e\u0026thinsp;=\u0026thinsp;0.001; Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e6\u003c/span\u003eD; Table S4). Ranks of adjusted \u003cem\u003eR\u003c/em\u003e\u003csup\u003e\u003cem\u003e2\u003c/em\u003e\u003c/sup\u003e values for relationships between RlogA metrics and log\u003csub\u003e2\u003c/sub\u003e(mat P) were higher than all other 18 metric relationships with mat P (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e6\u003c/span\u003eD).\u003c/p\u003e\u003cp\u003eUsing the same ANOVA tests of adjusted \u003cem\u003eR\u003c/em\u003e\u003csup\u003e\u003cem\u003e2\u003c/em\u003e\u003c/sup\u003e and Kruskal tests for ranks of metric-mat P relationships, TITAN-based metrics performed better than metrics based on literature or regression traits. Averages for adjusted \u003cem\u003eR\u003c/em\u003e\u003csup\u003e\u003cem\u003e2\u003c/em\u003e\u003c/sup\u003e for metric-mat P relationships, as well as their ranks, were lower for Lit than TITAN and regression sources of traits (Figs.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e6\u003c/span\u003eB \u0026amp; \u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e6\u003c/span\u003eE, Tukey HSD \u003cem\u003ep\u003c/em\u003e\u0026thinsp;\u0026lt;\u0026thinsp;0.026, Kruskal test \u003cem\u003ep\u003c/em\u003e\u0026thinsp;\u0026lt;\u0026thinsp;0.0001). Although there was great overlap in adjusted \u003cem\u003eR\u003c/em\u003e\u003csup\u003e\u003cem\u003e2\u003c/em\u003e\u003c/sup\u003e values of metric-mat P relationships based on traits determined by TITAN and regression (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e6\u003c/span\u003eB), for 7 out of 8 comparisons TITAN traits had higher ranks than regression traits for adjusted R\u003csup\u003e2\u003c/sup\u003e (\u003cem\u003ep\u003c/em\u003e\u0026thinsp;=\u0026thinsp;0.008, Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e6\u003c/span\u003eE). Neither ANOVA or Kruskal tests indicated differences in performance of low and high P metrics (Figs.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e6\u003c/span\u003eC, \u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e6\u003c/span\u003eF).\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003cp\u003eAs with metric-mat P relationships, all metrics except the diversity metrics Shannon H and Pielou\u0026rsquo;s evenness were related with high statistical significance and adjusted \u003cem\u003eR\u003c/em\u003e\u003csup\u003e\u003cem\u003e2\u003c/em\u003e\u003c/sup\u003e to our geographic indicators of human disturbance, latitude, longitude, and the lat-long interaction term (Table S5). In a comparison of which metrics and mat P were most precisely related to the lat-long human disturbance gradient, all metrics using regression and TITAN traits and both WAM metrics had greater adjusted \u003cem\u003eR\u003c/em\u003e\u003csup\u003e\u003cem\u003e2\u003c/em\u003e\u003c/sup\u003e than mat P. Mat P was related to the lat-long model with a 0.44 adjusted \u003cem\u003eR\u003c/em\u003e\u003csup\u003e\u003cem\u003e2\u003c/em\u003e\u003c/sup\u003e, whereas the four metrics with highest adjusted \u003cem\u003eR\u003c/em\u003e\u003csup\u003e\u003cem\u003e2\u003c/em\u003e\u003c/sup\u003e, which included both WAM metrics, had adjusted \u003cem\u003eR\u003c/em\u003e\u003csup\u003e\u003cem\u003e2\u003c/em\u003e\u003c/sup\u003e ranging from 0.57 to 0.65 (Table S5). The metric with the high 0.65 \u003cem\u003eR\u003c/em\u003e\u003csup\u003e\u003cem\u003e2\u003c/em\u003e\u003c/sup\u003e for a lat-long model was RlogA high P PropValves.\u003c/p\u003e\u003cp\u003eMetric testing and evaluation: metric differences among regions\u003c/p\u003e\u003cp\u003eAll metrics differed among regions with high statistical significance in 2020, with the highest P effects in the OK region and lowest in the MN region (Table S6, Figs.\u0026nbsp;\u003cspan refid=\"Fig7\" class=\"InternalRef\"\u003e7\u003c/span\u003e and S6). In OK, metric values were highest for diatom inferred mat P, Shannon\u0026rsquo;s H, proportion of taxa with high P traits and relative log-transformed abundance of high P taxa valves. In OK, metric values were lowest compared to all other regions for proportion of taxa with low P traits and log-transformed valve abundances with low P traits.\u003c/p\u003e\u003cp\u003eOverall, OK and then EH had the highest indicators of P pollution based on relatively high richness, WAM mat P metrics, and high-P trait metrics as well as relatively low low-P metrics compared to other regions (Figs.\u0026nbsp;\u003cspan refid=\"Fig7\" class=\"InternalRef\"\u003e7\u003c/span\u003e and S6). MN and LM had the consistently lowest richness, WAM mat P, and high P metrics as well as the highest low P metrics. BI, EC, and FP had intermediate levels of disturbance related to the other two groups of regions (Figs.\u0026nbsp;\u003cspan refid=\"Fig7\" class=\"InternalRef\"\u003e7\u003c/span\u003e and S6).\u003c/p\u003e\u003cp\u003eMetric performance for detecting differences among regions was greater for the 2 diatom inferred mat P metrics (WAM_matP, WAM_RlogA.matP) compared to all other metrics as indicated by ANOVA F values (Table S6). A Kruskal test indicated ranks of metric types differed (\u003cem\u003ep\u003c/em\u003e\u0026thinsp;=\u0026thinsp;0.025), with noTaxa and RlogA metrics having more highly ranked F values than PropTaxa and PropValves metrics (Fig. S7). Metric performance for distinguishing regions was better for TITAN and regression traits than literature traits (Fig. S7, Kruskal test, \u003cem\u003ep\u003c/em\u003e\u0026thinsp;=\u0026thinsp;0.009). There was no difference in low and high P metrics ability to detect differences among regions (Fig. S7). Mat P was significantly different among regions (ANOVA, \u003cem\u003ep\u0026thinsp;\u0026lt;\u003c/em\u003e\u0026thinsp;0.0001).\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003cp\u003eMetric testing and evaluation: corrections for natural factors\u003c/p\u003e\u003cp\u003eThe following results describing the relationship between residuals from metric-mat P relationships and these habitat parameters should be considered with caution. Unequal sample sizes among the natural habitat features presented challenges for some residual analyses. Whereas substrate locations had from 11\u0026ndash;21 sites for either floating, ground, or macrophyte sweaters, of the 6 habitat types there were 33 graminoid marshes and from 2\u0026ndash;3 of all other habitat types. For substrate types, there was only one soil sample when filamentous algae and periphyton had bigger sample sizes, 7 and 36 samples, respectively. Analyzing differences between the two substrate types, filamentous algae and periphyton, was an issue because shifts from periphyton to filamentous algae are associated with increasing P concentrations. Similarly, in this dataset water depth can be associated with human disturbance.\u003c/p\u003e\u003cp\u003eLittle evidence indicated that natural factors affected metric-mat P relationships (Tables S7, S8). Water depth was not significantly (\u003cem\u003ep\u003c/em\u003e\u0026thinsp;\u0026lt;\u0026thinsp;0.05) related to residuals in any metric-mat P relationship (Table S7). Time of year sampled (indicated by the variable named season) had a significant (\u003cem\u003ep\u003c/em\u003e\u0026thinsp;\u0026lt;\u0026thinsp;0.05) negative effect on 6 of 29 metrics, which were metrics for the numbers of high P taxa and diversity (Table S7). For substrate type, location, and hydrologic year, \u003cem\u003ep\u003c/em\u003e values for statistical significance were seldom less than 0.05 (1\u0026ndash;4 of 29 metrics, respectively) and not less than 0.05 if accounting for the multiple statistical tests conducted (Table S8). For habitat type, 5 of 29 metrics had a \u003cem\u003ep\u003c/em\u003e value less than 0.05, but low sample size for habitats other than graminoid marshes limited certainty in this observation.\u003c/p\u003e\u003cp\u003eEvidence was weak for natural factors affecting metrics after effects of human disturbance (lat-long model) were accounted for. Residuals from relationships relating either metrics or mat P with the lat-long model were so seldom related to natural factors with ANOVA or regression that the few tests with p\u0026thinsp;\u0026lt;\u0026thinsp;0.05 could have occurred by chance. Attained significance (\u003cem\u003ep\u003c/em\u003e) was not greater than 0.05 for any of the mat P or metric residuals relationships with water depth, time of season when sampling occurred, or substrate type. \u003cem\u003ep\u003c/em\u003e was \u0026lt;\u0026thinsp;0.05 for 1, 5, and 8 of the residual relationships with habitat type, substrate location, and hydrologic year, respectively (Table S8). Thus, substrate location and hydrologic year were the most likely natural factors affecting mat P and metrics. Interestingly, lat-long related residuals for both WAM metrics and mat P were related to hydrologic year, but residuals for regression or TITAN metrics were not related to hydrologic year (Table S8).\u003c/p\u003e\u003cp\u003eManagement benchmarks\u003c/p\u003e\u003cp\u003eGiven that metrics were evaluated along a log-transformed mat P gradient to even distribution of observations along the mat P gradient for regression analysis, that affected linearity of relationships to mat P. Metrics were linearly related to log transformed mat P according to comparisons of linear and non-linear models (e.g. Figures\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e5\u003c/span\u003e, S8). However, most metrics were non-linearly related to the mat P gradient when mat P was not log transformed (Fig. S9). With untransformed mat P, metrics with low P taxa traits had a negative exponential relationship with mat P with rapid decreases in low ranges of mat P and little decrease in high ranges of mat P. In contrast, metrics with high P traits increased asymptotically with untransformed mat P with rapid metric increases in the low mat P range and little change in metrics in the high range of mat P.\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003cp\u003eBenchmarks for relatively sudden changes in taxonomic composition were observed along the mat P gradient for 4 of 13 metrics evaluated. We evaluated mat P relationships with the 3 diversity metrics, 2 WAM metrics, and 8 TITAN-trait metrics. We excluded trait-based metrics that used literature and regression traits from the following analysis to prevent redundancy with TITAN-trait metrics that were best related to log mat P. Changepoint analysis showed that a breakpoint for low P TITAN taxa occurred at 427 ug g\u003csup\u003e\u0026minus;\u0026thinsp;1\u003c/sup\u003e mat P, above which all low P taxa numbers were low (Fig.\u0026nbsp;\u003cspan refid=\"Fig8\" class=\"InternalRef\"\u003e8\u003c/span\u003e, Table S9). Changepoint analysis showed numerous statistically significant (\u003cem\u003ep\u003c/em\u003e\u0026thinsp;\u0026le;\u0026thinsp;0.01) breakpoints occurring at the same mat P concentrations for both the number of all taxa and high P TITAN taxa: 71, 137, 464, 627, and 1246 \u0026micro;g g\u003csup\u003e\u0026minus;\u0026thinsp;1\u003c/sup\u003e mat P (Table S9, Figs.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003e and \u003cspan refid=\"Fig8\" class=\"InternalRef\"\u003e8\u003c/span\u003eB). WAM-mat P had a changepoint at 564 \u0026micro;g g\u003csup\u003e\u0026minus;\u0026thinsp;1\u003c/sup\u003e mat P, above which almost all values were high (Fig. S8). Changepoints along the mat P gradient were not observed for any other diversity, WAM, or TITAN metrics (Table S9). However, review of graphs of metrics showed a strong changepoint\u0026thinsp;\u0026le;\u0026thinsp;500 \u0026micro;g g\u003csup\u003e\u0026minus;\u0026thinsp;1\u003c/sup\u003e mat P for proportion of low and high P literature-defined taxa that was particularly evident when plotted along a log mat P scale (Fig. S8d).\u003c/p\u003e\u003cp\u003eFrom the heat map (Fig. S4), changes in assemblages around 137 \u0026micro;g g\u003csup\u003e\u0026minus;\u0026thinsp;1\u003c/sup\u003e mat P appeared to be associated with the shift from the most likely low P co-dominants being \u003cem\u003eE. evergladianum, N. serpentiraphe\u003c/em\u003e, and \u003cem\u003eNitzschia\u003c/em\u003e sp. 1 to \u003cem\u003eE. microcephala\u003c/em\u003e and \u003cem\u003eB. ocalanensis\u003c/em\u003e. Invasion of some high P taxa to accompany low P taxa in samples was more frequent above 137 \u0026micro;g g\u003csup\u003e\u0026minus;\u0026thinsp;1\u003c/sup\u003e. Above 464 and especially above 627 \u0026micro;g g\u003csup\u003e\u0026minus;\u0026thinsp;1\u003c/sup\u003e, high P taxa typically were most abundant in samples with loss of the sensitive common low P taxa \u003cem\u003eE. evergladianum, Nitzschia sp. 1\u003c/em\u003e, and \u003cem\u003eM. calcarea\u003c/em\u003e.\u003c/p\u003e\u003cp\u003eTITAN analysis supported the relatively high certainty of change in assemblages from low to high P taxa in the 300\u0026ndash;600 \u0026micro;g g\u003csup\u003e\u0026minus;\u0026thinsp;1\u003c/sup\u003e range of mat P (Figs.\u0026nbsp;\u003cspan refid=\"Fig9\" class=\"InternalRef\"\u003e9\u003c/span\u003e and S10). The peak in sum of low P taxa z scores was approximately 250 \u0026micro;g g\u003csup\u003e\u0026minus;\u0026thinsp;1\u003c/sup\u003e mat P. The spread of the interquartile range of changepoints for low P taxa, i.e. from the 0.25\u0026ndash;0.75 cumulative frequency quartiles, was from 300\u0026ndash;400 \u0026micro;g g\u003csup\u003e\u0026minus;\u0026thinsp;1\u003c/sup\u003e mat P and from 400\u0026ndash;600 \u0026micro;g g\u003csup\u003e\u0026minus;\u0026thinsp;1\u003c/sup\u003e for high P taxa. The overlap in those changepoint ranges was evident in the heat map of species relative abundances in samples along the mat P gradient (Fig. S4). Between 400 and 500 \u0026micro;g g\u003csup\u003e\u0026minus;\u0026thinsp;1\u003c/sup\u003e mat P, as mat P increased, the frequency of high abundances of most low P taxa decreased greatly and yielded to an increase in numbers and abundances of high P taxa in samples (Fig. S4). However, some low P taxa had taxon-specific changepoints in indicator species values below the 300 \u0026micro;g g\u003csup\u003e\u0026minus;\u0026thinsp;1\u003c/sup\u003e mat P (Fig. S10). In contrast, some low P taxa were able to persist in relatively high mat P conditions, such as \u003cem\u003eE. microcephala\u003c/em\u003e and \u003cem\u003eB. ocalanensis\u003c/em\u003e that had likely indicator value changepoints greater than 500 \u0026micro;g g\u003csup\u003e\u0026minus;\u0026thinsp;1\u003c/sup\u003e mat P (Fig. S10). \u003cem\u003eEncyonema silesiacum\u003c/em\u003e and \u003cem\u003eN. cryptotenella\u003c/em\u003e had the opposite characteristics; they were high P taxa that were commonly in high abundance in low P conditions as quantified by their median indicator species thresholds near 350\u0026ndash;400 \u0026micro;g g\u003csup\u003e\u0026minus;\u0026thinsp;1\u003c/sup\u003e mat P (Fig. S10).\u003c/p\u003e\u003cp\u003e\u003c/p\u003e\u003cp\u003eMat P concentrations in minimally disturbed regions of the sampling area can be used to characterize the natural variability in minimally disturbed conditions. Minimum and maximum mat P were 52 and 606 \u0026micro;g g\u003csup\u003e\u0026minus;\u0026thinsp;1\u003c/sup\u003e in the calibration data set at the 32 sites sampled in the FP, LM, and MN regions. The median and 75th percentile of mat P values in FP, LM, and MN were 140 and 367 \u0026micro;g g\u003csup\u003e\u0026minus;\u0026thinsp;1\u003c/sup\u003e mat P. Sixty-four percent of mat P values at the 32 FP, LM, and MN sites were less than 200 \u0026micro;g g\u003csup\u003e\u0026minus;\u0026thinsp;1\u003c/sup\u003e, and 90% of mat P values were less than 500 \u0026micro;g g\u003csup\u003e\u0026minus;\u0026thinsp;1\u003c/sup\u003e.\u003c/p\u003e"},{"header":"Discussion","content":"\u003cp\u003eChanges in diatom assemblages along the P gradient\u003c/p\u003e\u003cp\u003eMultiple lines of evidence indicated that phosphorus concentration was the determinant of changes in diatom species composition in regions of the BCNP near areas of human disturbance, as in regions of the greater Everglades. Mat P and metrics were highly related to distance from human disturbance. NMDS analyses showed mat P was the variable most strongly related to variation in diatom species composition. The mat P optima for taxa determined with the BCNP calibration dataset correlated well with P optima identified in experimental manipulations of P in the Everglades Natural Park and Everglades Water Conservation Area 2A (Pan et al., \u003cspan citationid=\"CR33\" class=\"CitationRef\"\u003e2000\u003c/span\u003e; Gaiser et al., \u003cspan citationid=\"CR11\" class=\"CitationRef\"\u003e2006\u003c/span\u003e; Slate \u0026amp; Stevenson, \u003cspan citationid=\"CR43\" class=\"CitationRef\"\u003e2007\u003c/span\u003e). Phosphorus contamination in different regions of the Everglades, typically from canal water entering marshes, is widely recognized as a major threat to biological condition in the Everglades (McCormick et al., \u003cspan citationid=\"CR27\" class=\"CitationRef\"\u003e1996\u003c/span\u003e; Gaiser et al., \u003cspan citationid=\"CR11\" class=\"CitationRef\"\u003e2006\u003c/span\u003e) and many other wetland ecosystems as well (Pan et al., 1996; Lougheed et al., \u003cspan citationid=\"CR24\" class=\"CitationRef\"\u003e2007\u003c/span\u003e; Wyatt et al., \u003cspan citationid=\"CR69\" class=\"CitationRef\"\u003e2010\u003c/span\u003e; Pillsbury et al., \u003cspan citationid=\"CR36\" class=\"CitationRef\"\u003e2019\u003c/span\u003e).\u003c/p\u003e\u003cp\u003eThe number of taxa observed in BCNP samples (alpha diversity) increased with mat P concentration from an assemblage restricted to very few low P taxa to a wider diversity of taxa requiring high P concentrations, as with other regions of southern Florida (Raschke, \u003cspan citationid=\"CR40\" class=\"CitationRef\"\u003e1993\u003c/span\u003e; Pan et al., \u003cspan citationid=\"CR33\" class=\"CitationRef\"\u003e2000\u003c/span\u003e; Gaiser et al., \u003cspan citationid=\"CR11\" class=\"CitationRef\"\u003e2006\u003c/span\u003e). Given the widespread sampling in the greater Everglades region and low number of taxa observed in low P conditions, it seems likely that diversity of diatoms is indeed low, despite limitations of 600 valve counts. Low diversity in naturally low nutrient conditions could be due to few taxa are adapted to the stress of low nutrient supply (Worm et al., \u003cspan citationid=\"CR68\" class=\"CitationRef\"\u003e2002\u003c/span\u003e; Cardinale et al., \u003cspan citationid=\"CR6\" class=\"CitationRef\"\u003e2009\u003c/span\u003e). Then, as nutrient concentrations increase, habitat availability increases for a larger number of taxa that require higher P to persist. Some low P taxa were observed occasionally in moderate and even high P sites, however other low P taxa were not observed in high P sites. Although showing that taxa have indeed been extirpated from a habitat is difficult, especially for microbes for which we observe such minute proportions of populations, the lack of occurrence of taxa such as \u003cem\u003eE. evergladianum\u003c/em\u003e in the large number of counts with sites having more than 550 \u0026micro;g g\u003csup\u003e\u0026minus;\u0026thinsp;1\u003c/sup\u003e mat P indicates they are very rare in high P areas even though they are the most common diatom in minimally disturbed condition.\u003c/p\u003e\u003cp\u003eDiversity of taxa among sites (beta diversity) also increased with P concentration. Cluster analysis showed high similarity among sites occurred in low P conditions, and dissimilarity in taxonomic composition was relatively high among high P sites. This differs from other studies that show human disturbance causes greater homogenization of flora among sites (reduced beta diversity, Lougheed et al., \u003cspan citationid=\"CR23\" class=\"CitationRef\"\u003e2008\u003c/span\u003e). The mechanisms regulating apparent and absolute diversity probably vary. The low P conditions in the greater Everglades that constrain taxonomic membership to a small number of taxa may be an unusually low resource condition not observed along human disturbance gradients in other locations. Thus, in the greater Everglades and BCNP we have low nutrient supply constraining taxa numbers with release of nutrient constraint as P supply increases allowing other taxa to invade, thereby increasing apparent numbers of taxa in the habitat. In the higher range of a resource gradient that can be found in other non-Everglades locations with high human disturbance (Pan et al., 1996; Wang et al., \u003cspan citationid=\"CR64\" class=\"CitationRef\"\u003e2006\u003c/span\u003e; Lougheed et al, \u003cspan citationid=\"CR24\" class=\"CitationRef\"\u003e2007\u003c/span\u003e), nutrients may be sufficiently high at the high end of the human disturbance gradient so additional disturbance allows overgrowth and dominance of a few high nutrient species, thereby homogenizing biota with increasing human disturbance.\u003c/p\u003e\u003cp\u003eTaxa traits and relative metric performance\u003c/p\u003e\u003cp\u003eMore taxa were indicative of high P than low P conditions in BCNP. This is likely related to the low mat P conditions at some BCNP sites, which are among the lowest in the greater Everglades area. All taxa identified as low P and high P with TITAN and regression were characterized the same. In contrast, there were discrepancies for low and high P characterizations from the literature evaluations by NPS for BCNP when compared to the TITAN and regression traits determined with the calibration dataset. Some of the discrepancies between NPS literature and BCNP data evaluations were due to different taxonomic treatments of taxa, but many were related to the literature dealing with higher nutrient ranges, so low P (oligotrophic) taxa in the literature were from relatively low P conditions in studies which happen to be intermediate or high P conditions for BCNP and the greater Everglades. Taxa P optima were highly correlated among different Everglades projects (Gaiser et al., \u003cspan citationid=\"CR11\" class=\"CitationRef\"\u003e2006\u003c/span\u003e; Slate \u0026amp; Stevenson \u003cspan citationid=\"CR43\" class=\"CitationRef\"\u003e2007\u003c/span\u003e; Hagerthey et al., \u003cspan citationid=\"CR17\" class=\"CitationRef\"\u003e2012\u003c/span\u003e) despite the limited number of pairwise taxa comparisons across studies that were possible. The concordance of these multiple lines of evidence again indicate that phosphorus is a causal determinant of low and high P characterizations assigned to taxa.\u003c/p\u003e\u003cp\u003eAlmost all metrics were strongly related to mat P, whether traits were derived by regression and TITAN with BCNP data or by reviewing literature. Almost all metrics with all traits differed among BCNP regions with differing levels of human disturbance. However, metrics using traits derived using regression and TITAN consistently performed better than metrics calculated with traits derived from the NPS literature evaluation. Tests with the validation dataset and the 2020 regional comparison indicated that regression and TITAN traits performed better in metrics than literature traits. As in other studies, project specific taxonomy may affect trait assignments to species derived from the literature, but we also show environmental gradient lengths and ranges also affect trait assignment and metric performance. In a comparison of regression and TITAN, TITAN traits consistently performed better.\u003c/p\u003e\u003cp\u003eMetrics calculated with the relative log abundance of taxa performed better than metrics calculated as the proportion of individuals or taxa. This was observed with the validation data tests of metrics relationships to mat P and differences among regions with differing mat P. Applications of simpler metrics calculated with the proportion of taxa having a trait may have good precision (low variability) when related to disturbance (Stevenson et al., \u003cspan citationid=\"CR52\" class=\"CitationRef\"\u003e2008\u003c/span\u003eb) and can be easier to relate conceptually to standard characteristics of biological condition as species loss or invasion (Davies et al., 2006), however resource managers should know that accuracy (proximity to true value) of these metrics is a concern. The actual proportions of low or high disturbance taxa could be very different if more thorough assessments than 600 valve counts were used to determine taxon presence or absence, as illustrated by the large numbers of rare taxa in log-normal, taxon-abundance distributions (Preston, \u003cspan citationid=\"CR38\" class=\"CitationRef\"\u003e1948\u003c/span\u003e; Patrick, \u003cspan citationid=\"CR34\" class=\"CitationRef\"\u003e1967\u003c/span\u003e).\u003c/p\u003e\u003cp\u003eTraits based metrics for pollutants can be more precise indicators of pollution levels than measurements of pollutants that are difficult to measure in the environment because of spatial and temporal variability (Stevenson et al., \u003cspan citationid=\"CR54\" class=\"CitationRef\"\u003e2010\u003c/span\u003e). In BCNP, WAM and metrics using either TITAN or regression-derived taxa P traits were more precisely related to the human disturbance gradient than measured mat P. Additionally, P concentrations in algae (mat P) or sediments are more precise and sensitive indicators of wetland P availability than measurements of P concentration in the water column (Pan et al., \u003cspan citationid=\"CR33\" class=\"CitationRef\"\u003e2000\u003c/span\u003e; Gaiser et al., \u003cspan citationid=\"CR12\" class=\"CitationRef\"\u003e2004\u003c/span\u003e). Thus, basing management targets on biological indicators, either as measures of biological condition or as indicators of P pollution, warrants consideration.\u003c/p\u003e\u003cp\u003eMetrics were affected little by natural factors such as wetland habitat type, substrate type, substrate location, water depth, and when sampling occurred during the high water season. Analyses of both sets of residuals, for metric-mat P relationships or metric and mat P relationships with the lat-long models of human disturbance, showed little likelihood that natural factors had major effects on metrics. Metrics are often robust to pollution levels and human disturbance, but there are study scales when we need to account for how minimally disturbed condition of pollutants and metrics vary with natural factors, such as in streams and lakes assessments (Cao et al., \u003cspan citationid=\"CR5\" class=\"CitationRef\"\u003e2007\u003c/span\u003e; Stevenson et al., \u003cspan citationid=\"CR55\" class=\"CitationRef\"\u003e2013\u003c/span\u003e; Tang et al., \u003cspan citationid=\"CR58\" class=\"CitationRef\"\u003e2016\u003c/span\u003e). But the range of natural factors in BCNP had little effect on metric relationships to mat P (the stressor) or the lat-long model of human disturbance.\u003c/p\u003e\u003cp\u003eHydrologic year did affect expected levels of mat P and metrics in BCNP according to differences in residuals for the metric and mat P relationships with the lat-long model. All years of samples were counted by the same taxonomist for data used in metric testing and evaluation, which reduces the likelihood that taxonomist caused interannual variability in residuals. Interannual differences in P availability seem to be the issue because mat P residuals also differed with hydrologic year. Perhaps weather-related events, such as long periods of stable weather or recent rainfalls, affected P availability and thereby, metric values.\u003c/p\u003e\u003cp\u003eManagement benchmarks and metric considerations\u003c/p\u003e\u003cp\u003eManagement benchmarks for protecting ecological condition call for determining how ecological systems respond to pollutants and relating those responses to different levels of human disturbance to prevent over- or under-protection of a resource. We were particularly interested in thresholds (tipping points) of ecological responses for their potential to develop consensus among stakeholders for management targets (Muradian, \u003cspan citationid=\"CR29\" class=\"CitationRef\"\u003e2001\u003c/span\u003e), such as levels of biological condition to protect and corresponding levels of pollution that protect biological condition. Many metrics had nonlinear relationships with untransformed mat P concentrations, either increasing or decreasing asymptotically. These responses are not as helpful for identifying threshold responses as are responses having some assimilative capacity at low levels of stressors where biota respond little because they are so highly sensitive to even small changes in human disturbance (Fig. \u003cspan refid=\"MOESM1\" class=\"InternalRef\"\u003eS1\u003c/span\u003e, Stevenson, \u003cspan citationid=\"CR48\" class=\"CitationRef\"\u003e2011\u003c/span\u003e). Decision tree analyses did show biological responses with assimilative capacity, such as high P taxa increasing when mat P surpassed 137 \u0026micro;g g\u003csup\u003e\u0026minus;\u0026thinsp;1\u003c/sup\u003e mat P, and low P taxa decreasing when mat P exceeded 464 \u0026micro;g g\u003csup\u003e\u0026minus;\u0026thinsp;1\u003c/sup\u003e. These mat P concentrations showing assimilative capacity provide breakpoints in relationships that are particularly valuable for justifying management targets in BCNP. In addition, groups of sites with similar diatom assemblage composition were observed in low P conditions, e.g. Group 1\u0026ndash;1 was the low P subgroup within Group 1. The 75th percentile of mat P for Group 1\u0026ndash;1 was 142 and for all of Group 1 (that included subgroups 1\u0026ndash;1, 1\u0026ndash;2, and 1\u0026ndash;3) was 354 \u0026micro;g g\u003csup\u003e\u0026minus;\u0026thinsp;1\u003c/sup\u003e. TITAN analysis showed transition from low to high P taxa occurring in the 300\u0026ndash;600 \u0026micro;g g\u003csup\u003e\u0026minus;\u0026thinsp;1\u003c/sup\u003e mat P range with most transition of low P taxa in the 300\u0026ndash;400 \u0026micro;g g\u003csup\u003e\u0026minus;\u0026thinsp;1\u003c/sup\u003e range.\u003c/p\u003e\u003cp\u003eBreakpoints in ecological responses of calcareous periphyton in the eastern, peat-based wetlands of the Everglades could be relevant for management of BCNP because calcareous periphyton with the same species is also expected in minimally disturbed conditions of BCNP. Much of the early research in the Everglades relates periphyton responses to water column TP. Floating mat cover and species composition changed greatly with elevation of TP above 10 \u0026micro;g L\u003csup\u003e\u0026minus;\u0026thinsp;1\u003c/sup\u003e in WCA-2A (McCormick, et al. \u003cspan citationid=\"CR27\" class=\"CitationRef\"\u003e1996\u003c/span\u003e; Stevenson et al., \u003cspan citationid=\"CR53\" class=\"CitationRef\"\u003e2002\u003c/span\u003e), which was important for justifying Florida\u0026rsquo;s 10 \u0026micro;g TP L\u003csup\u003e\u0026minus;\u0026thinsp;1\u003c/sup\u003e criterion for the Everglades. According to Hagerthey et al. (\u003cspan citationid=\"CR16\" class=\"CitationRef\"\u003e2008\u003c/span\u003e), mat P would be 412 \u0026micro;g g\u003csup\u003e\u0026minus;\u0026thinsp;1\u003c/sup\u003e in sloughs with a water column TP of 11 \u0026micro;g L\u003csup\u003e\u0026minus;\u0026thinsp;1\u003c/sup\u003e, which was a threshold for slough macrophytes. Pan et al. (\u003cspan citationid=\"CR33\" class=\"CitationRef\"\u003e2000\u003c/span\u003e) found a characteristic low P periphyton assemblage (TWINSPAN I) at sites with epiphytic mat P averaging 290 \u0026micro;g g\u003csup\u003e\u0026minus;\u0026thinsp;1\u003c/sup\u003e \u0026plusmn; 50 SE, which with n\u0026thinsp;=\u0026thinsp;10 would have a 75th mat P percentile of 300 \u0026micro;g g\u003csup\u003e\u0026minus;\u0026thinsp;1\u003c/sup\u003e. Gaiser et al. (\u003cspan citationid=\"CR13\" class=\"CitationRef\"\u003e2005\u003c/span\u003e) found enrichment of background conditions in long-term experiments by as little as 5 \u0026micro;g P L\u003csup\u003e\u0026minus;\u0026thinsp;1\u003c/sup\u003e can alter periphyton assemblages, and therefore, the metabolic biogeochemical processes that regulate periphyton. This sensitivity to low concentrations of P is not surprising given the asymptotic increase in periphyton growth rates and high sensitivity when P is less than 10 \u0026micro;g L\u003csup\u003e\u0026minus;\u0026thinsp;1\u003c/sup\u003e in stream water (Bothwell, \u003cspan citationid=\"CR3\" class=\"CitationRef\"\u003e1989\u003c/span\u003e; Rier et al., 2006).\u003c/p\u003e\u003cp\u003eThe substantial transitions from low to high P taxa in BCNP periphyton and in the Everglades occurring in the 300\u0026ndash;400 \u0026micro;g g\u003csup\u003e\u0026minus;\u0026thinsp;1\u003c/sup\u003e mat P range corresponded to the 75th percentile of mat P in minimally disturbed conditions of BCNP, which was 367 \u0026micro;g g\u003csup\u003e\u0026minus;\u0026thinsp;1\u003c/sup\u003e. Differences in periphyton in ranges of mat P less than the 75th percentile of mat P in the three most southern and eastern minimally disturbed regions of BCNP indicated that habitat specific management targets could be important. For example, if the mat P management target were set at the 367 \u0026micro;g g\u003csup\u003e\u0026minus;\u0026thinsp;1\u003c/sup\u003e mat P, then Group 1\u0026ndash;1 assemblages found in very low P conditions would not be protected, even though the more inclusive Group 1 assemblages would be protected. Gaiser et al. (\u003cspan citationid=\"CR11\" class=\"CitationRef\"\u003e2006\u003c/span\u003e) characterized minimally disturbed conditions for Shark River Slough and Taylor Slough of the Everglades as having maximum TP at levels that would protect the Group 1\u0026ndash;1 periphyton assemblages, so such habitats do occur on a large scale in southern Florida. Further research should be conducted to determine if lower mat P management targets should be used to protect some habitat types with the Group 1\u0026ndash;1 periphyton assemblages. Natural habitat features could create different P conditions, which could justify setting different management targets for different kinds of wetland habitats, even though our analyses showed residuals in the metric mat P relationships with lat-long were not related to natural factors (e.g. local hydrologic conditions).\u003c/p\u003e\u003cp\u003eMany factors should be considered when selecting metrics for use in a monitoring program (Jackson et al., \u003cspan citationid=\"CR18\" class=\"CitationRef\"\u003e2000\u003c/span\u003e). All BCNP metrics were related to mat P, differed among BCNP regions in 2020, and were related to the human disturbance gradient, in most cases. Therefore, many options could be satisfactory for use in a monitoring program. Distinguishing between metrics measuring an element of biological condition (low or high P taxa indicators) versus inferring stressors (WAM mat P models) is important conceptually relative to the idea that we want to manage stressors to protect biological condition and have clear communication with stakeholders about levels of protection for measures of biological condition (Stevenson, \u003cspan citationid=\"CR48\" class=\"CitationRef\"\u003e2011\u003c/span\u003e). The higher performance of TITAN and RlogA metrics of biological condition introduces tradeoffs between using metrics that perform best (i.e., are the most accurately and precisely related to mat P concentration, changes among regions, and human disturbance) and metrics that are simpler to explain to stakeholders, i.e. how they were calculated and relate to ecological condition. Simple metrics such as numbers or even proportions of low and high P taxa may also be more vulnerable to taxonomic error, if the number of taxa distinguished among taxonomists differs. There are relatively simple ways to explain TITAN and RlogA methods and why they are likely better than other methods, but nevertheless, they are more complicated methods.\u003c/p\u003e\u003cp\u003eFortunately, a library of taxa traits and a database of diatom counts allows calculation of all the metrics relatively easily, so some metrics can be used for public reports and others can be used for administrative and technical records and analysis. Such a set of tools, as we provided with this research, provides the necessary information for advancing a monitoring and assessment program and establishing management targets for biological condition and stressor levels that will protect an unimpaired condition, which is the specified goal of BCNP management (Patterson et al., \u003cspan citationid=\"CR35\" class=\"CitationRef\"\u003e2008\u003c/span\u003e).\u003c/p\u003e"},{"header":"Conclusions","content":"\u003cp\u003eWe have developed the information needed for a monitoring and assessment program and integrated it for BCNP with many attributes that are not included in other programs. First we explored the ecological relationships among indicators of biological condition and pollutants which are directly related to the goals of management of the resource. Then we developed multiple metrics of biological condition and tested which are the most sensitive and consistently related to the primary pollutant. We confirmed the cause-effect relationship between metrics of biological condition and P, the primary pollutant. Then we determined the natural variation in minimally disturbed condition and characterized changes in biological condition from minimally disturbed condition with increasing P concentrations which can be used to justify effects-based management targets. These are fundamental steps for developing an effects-based monitoring and assessment program, but all of them are seldom integrated and presented in one publication with state-of-the-art approaches to show clear relationships among the lines of results needed for an integrated monitoring and assessment program.\u003c/p\u003e\u003cp\u003eWe advanced the state-of-the-art of monitoring and assessment at multiple steps in our research. We showed how nutrient enrichment clearly affected biological assemblages by allowing invasion of taxa that cannot survive in the naturally low nutrient conditions and how sensitive low nutrient taxa are greatly reduced or eliminated as nutrient concentrations increase. We confirmed cause-effect relationships using multiple causal assessment methods. We developed and tested a new indicator calculation method that improved metric performance and was based on basic principles of population growth dynamics using log transformed relative abundances. And finally, we integrated the application of minimally disturbed condition and biological condition response along a pollution gradient to provide a potentially ultimate management goal of minimally disturbed condition, but also interim goals to assess incremental restoration of the ecosystem.\u003c/p\u003e"},{"header":"Declarations","content":"\u003ch2\u003eFunding\u003c/h2\u003e\u003cp\u003eThe United States National Park Service funded this research.\u003c/p\u003e\u003ch2\u003eAuthor Contribution\u003c/h2\u003e\u003cp\u003eK.R.T.W. was responsible for N.P.S. administration of the project. K.R.T.W. and M.C.P. collected samples. R.J.S. analyzed samples and data, and was lead writer of the manuscript. All authors contributed to writing and reviewing the manuscript.\u003c/p\u003e\u003ch2\u003eAcknowledgement\u003c/h2\u003e\u003cp\u003eThe National Park Service (NPS) funded this research. Raul Urgelles and Andrea Atkinson at NPS were founding members of this monitoring project. Mario Londono at the NPS provided database assistance. The NPS Big Cypress Aviation unit supported sampling missions; and specially we would like to thank Fred Goodwin (pilot) and Michael Roof (Helicopter Program Manager).\u003c/p\u003e\u003ch2\u003eData Availability\u003c/h2\u003e\u003cp\u003eThe data that support the findings of this study are not openly available because they are being prepared for deposit in an appropriate archive. They will be published before this manuscript is accepted for publication.\u003c/p\u003e"},{"header":"References","content":"\u003col\u003e\n\u003cli\u003eBaker, M. E., \u0026amp; King, R. S. (2010). A new method for detecting and interpreting biodiversity and ecological community thresholds. \u003cem\u003eMethods in Ecology and Evolution, 1\u003c/em\u003e(1), 25-37. doi:https://doi.org/10.1111/j.2041-210X.2009.00007\u003c/li\u003e\n\u003cli\u003eBeyers, D. W. (1998). Causal inference in environmental impact studies. \u003cem\u003eJournal of the North American Benthological Society, 17\u003c/em\u003e, 367-373. \u003c/li\u003e\n\u003cli\u003eBothwell, M. L. (1989). Phosphorus-limited growth dynamics of lotic periphytic diatom communities: areal biomass and cellular growth rate responses. \u003cem\u003eCanadian Journal of Fisheries and Aquatic Sciences, 46\u003c/em\u003e, 1293-1301. \u003c/li\u003e\n\u003cli\u003eBrowder, J. A., Gleason, P. J., \u0026amp; Swift, D. R. (1994). Periphyton in the Everglades: spatial variation, environmental correlates,and ecological implications. In S. M. Davis \u0026amp; J. C. Ogden (Eds.), \u003cem\u003eEverglades, the ecosystem and its restoration\u003c/em\u003e (pp. pp.379-418). Delray Beach, Florida: St. Lucie Press.\u003c/li\u003e\n\u003cli\u003eCao, Y., Hawkins, C. P., Olson, J., \u0026amp; Kosterman, M. A. (2007). Modeling natural environmental gradients improves the accuracy and precision of diatom-based indicators. \u003cem\u003eJournal of the North American Benthological Society, 26\u003c/em\u003e(3), 566-585. \u003c/li\u003e\n\u003cli\u003eCardinale, B. J., Hillebrand, H., Harpole, W. S., Gross, K., \u0026amp; Ptacnik, R. (2009). Separating the influence of resource \u0026lsquo;availability\u0026rsquo; from resource \u0026lsquo;imbalance\u0026rsquo; on productivity\u0026ndash;diversity relationships. \u003cem\u003eEcology Letters, 12\u003c/em\u003e(6), 475-487. doi:https://doi.org/10.1111/j.1461-0248.2009.01317\u003c/li\u003e\n\u003cli\u003eCarlisle, D. M., Spaulding, S. A., Tyree, M. A., Schulte, N. O., Lee, S. S., Mitchell, R. M., \u0026amp; Pollard, A. A. (2022). A web-based tool for assessing the condition of benthic diatom assemblages in streams and rivers of the conterminous United States. \u003cem\u003eEcological Indicators, 135\u003c/em\u003e, 108513. doi:https://doi.org/10.1016/j.ecolind.2021.108513\u003c/li\u003e\n\u003cli\u003eCharles, D. F., Kelly, M. G., Stevenson, R. J., Poikane, S., Theroux, S., Zgrundo, A., \u0026amp; Cantonati, M. (2021). Benthic algae assessments in the EU and the US: Striving for consistency in the face of great ecological diversity. \u003cem\u003eEcological Indicators, 121\u003c/em\u003e, 107082. doi:https://doi.org/10.1016/j.ecolind.2020.107082\u003c/li\u003e\n\u003cli\u003eCormier, S. M., \u0026amp; Suter II, G. W. (2008). A framework for fully integrating environmental assessment. \u003cem\u003eEnvironmental Management, 42\u003c/em\u003e, 543-556. \u003c/li\u003e\n\u003cli\u003eDavies, S. P., \u0026amp; Jackson, S. K. (2006). The biological condition gradient: a descriptive model for interpreting change in aquatic ecosystems. \u003cem\u003eEcological Applications, 16\u003c/em\u003e, 1251-1266. \u003c/li\u003e\n\u003cli\u003eGaiser, E. E., Childers, D. L., Jones, R. D., Richards, J. H., \u0026amp; Scinto, L. J. (2006). Periphyton responses to eutrophication in the Florida Everglades: cross-system patterns of strucural and composition changes. \u003cem\u003eLimnology \u0026amp; Oceanography, 51\u003c/em\u003e, 617-630. \u003c/li\u003e\n\u003cli\u003eGaiser, E. E., L. J. Scinto, J. H. Richards, K. Jayachandran, D. L. Childers, J. D. Trexler, \u0026amp; Jones, R. D. (2004). Phosphorus in periphyton mats provides best metric for detecting low-level P enrichment in an oligotrophic wetland. \u003cem\u003eWater Research, 38\u003c/em\u003e, 507-516. \u003c/li\u003e\n\u003cli\u003eGaiser, E. E., Trexler, J. C., Richards, J. H., Childers, D. L., Lee, D., Edwards, A. L., . . . Jones, R. D. (2005). Cascading Ecological Effects of Low-Level Phosphorus Enrichment in the Florida Everglades. \u003cem\u003eJournal of Environmental Quality, 34\u003c/em\u003e(2), 717-723. doi:https://doi.org/10.2134/jeq2005.0717\u003c/li\u003e\n\u003cli\u003eGleason, P. J., \u0026amp; W. Spackman, J. (1974). Calcareous periphyton and water chemistry in the Everglades. In P. J. Gleason (Ed.), \u003cem\u003eEnvironments of South Florida: Present and Past, Memoir No. 2.\u003c/em\u003e (pp. pp. 146-181). Caral Gables, Fla.: Miami Geological Society.\u003c/li\u003e\n\u003cli\u003eHagerthey, S. E., Bellinger, B. J., Wheeler, K., Gantar, M., \u0026amp; Gaiser, E. (2011). Everglades Periphyton: A Biogeochemical Perspective. \u003cem\u003eCritical Reviews in Environmental Science and Technology, 41\u003c/em\u003e(sup1), 309-343. doi:10.1080/10643389.2010.531218\u003c/li\u003e\n\u003cli\u003eHagerthey, S. E., Newman, S., Rutchey, K., Smith, E. P., \u0026amp; Godin, J. (2008). Multiple regime shifts in a subtropical peatland: community-specific thresholds in eutrophication. \u003cem\u003eEcological Monographs, 78\u003c/em\u003e(4), 547-565. doi:https://doi.org/10.1890/07-0538.\u003c/li\u003e\n\u003cli\u003eHagerthey, S. E., Newman, S., \u0026amp; Xue, S. (2012). Periphyton-based transfer functions to assess ecological imbalance and management of a subtropical ombrotrophic peatland. \u003cem\u003eFreshwater Biology, 57\u003c/em\u003e, 1947-1965. \u003c/li\u003e\n\u003cli\u003eJackson, L. E., Kutz, J. C., \u0026amp; Fisher, W. S. (Eds.). (2000). \u003cem\u003eEvaluation guidelines for ecological indicators\u003c/em\u003e. Research Triangle Park, North Carolina: United States Environmental Protection Agency.\u003c/li\u003e\n\u003cli\u003eJarosewich, M. \u0026amp; Wagner, J. (1985) Geologic Structure of the Surficial Aquifer System Underlying Everglades National Park and the Big Cypress National Preserve. South Florida Research Center Report SFRC-85/02. 43p.\u003c/li\u003e\n\u003cli\u003eLa H\u0026eacute;e, J. M., \u0026amp; Gaiser, E. E. (2012). Benthic diatom assemblages as indicators of water quality in the Everglades and three tropical karstic wetlands. \u003cem\u003eFreshwater Science, 31\u003c/em\u003e, 205-221. \u003c/li\u003e\n\u003cli\u003eLavoie, I., Hamilton, P. B., Wang, Y. K., Dillon, P. J., \u0026amp; Campeau, S. (2009). A comparison of stream bioassessment in Qu\u0026eacute;bec (Canada) using six European and North American diatom-based indices. \u003cem\u003eNova Hedwigia, 135\u003c/em\u003e, 37-56. \u003c/li\u003e\n\u003cli\u003eLondo\u0026ntilde;o, M. 2019. Selecting Sites in GIS\u0026mdash;Version 1.0. South Florida/Caribbean Network Standard Operating Procedure. South Florida/Caribbean Network, Palmetto Bay, Florida. Available at: https://irma.nps.gov/DataStore/Reference/Edit/2260145 (last accessed April 2019).\u003c/li\u003e\n\u003cli\u003eLougheed, V. L., McIntosh, M. D., Parker, C. A., \u0026amp; Stevenson, R. J. (2008). Wetland degradation leads to homogenization of biota at local and landscape scales. \u003cem\u003eFreshwater Biology, 52\u003c/em\u003e, 2402-2413. \u003c/li\u003e\n\u003cli\u003eLougheed, V. L., Parker, C. A., \u0026amp; Stevenson, R. J. (2007). Using non-linear responses of multiple taxonomic groups to establish criteria indicative of wetland biological condition. \u003cem\u003eWetlands, 27\u003c/em\u003e, 96-109. \u003c/li\u003e\n\u003cli\u003eMcCormick, P. V., Chimney, M. J., \u0026amp; Swift, D. R. (1997). Diel oxygen profiles and water column community metabolism in the Florida Everglades, USA. \u003cem\u003eArchiv Fur Hydrobiologie, 140\u003c/em\u003e(1), 117-129. \u003c/li\u003e\n\u003cli\u003eMcCormick, P. V., O\u0026apos;Dell, M. B., Shuford, R. B. E., Backus, J. G., \u0026amp; Kennedy, W. C. (2001). Periphyton responses to experimental phosphorus enrichment in a subtropical wetland. \u003cem\u003eAquatic Botany, 71\u003c/em\u003e(2), 119-139. \u003c/li\u003e\n\u003cli\u003eMcCormick, P. V., Rawlick, P. S., Lurding, K., Smith, E. P., \u0026amp; Sklar, F. H. (1996). Periphyton-water quality relationships along a nutrient gradient in the northern Florida Everglades. \u003cem\u003eJournal of the North American Benthological Society, 15\u003c/em\u003e, 433-449. \u003c/li\u003e\n\u003cli\u003eMiller, R., McPherson, B., Sobczak, R., \u0026amp; Clark, C. (2004) Water Quality in Big Cypress National Preserve and Everglades National Park \u0026mdash; Trends and Spatial Characteristics of Selected Constituents. U.S. Geological Survey (USGS) Water Resources Investigations Report. 03-4249. U.S. Geological Survey; U.S. Department of the Interior. Denver, Colorado.\u003c/li\u003e\n\u003cli\u003eMuradian, R. (2001). Ecological thresholds: a survey. \u003cem\u003eEcological Economics, 38\u003c/em\u003e, 7-24. \u003c/li\u003e\n\u003cli\u003eMuxo, R. \u0026amp; Shamblin, R. B. (2019) Standard operating procedure 4: Periphyton Sample Collection Methods. South Florida/Caribbean Network, Palmetto Bay, Florida. Available at: https://irma.nps.gov/DataStore/Reference/Edit/2260145 (last accessed April 2019).\u003c/li\u003e\n\u003cli\u003eOksanen, J., F. Guillaume Blanchet, M. Friendly, R. Kindt, P. Legendre, D. McGlinn, . . . Wagner, H. (2020). vegan: Community Ecology Package (Version R package version 2.5-7). Retrieved from https://CRAN.R-project.org/package=vegan\u003c/li\u003e\n\u003cli\u003e\u003c/li\u003e\n\u003cli\u003ePan, Y., \u0026amp; Stevenson, R. J. (1996). Gradient analysis of diatom assemblages in Western Kentucky wetlands. \u003cem\u003eJournal of Phycology, 32\u003c/em\u003e, 222-232. \u003c/li\u003e\n\u003cli\u003ePan, Y., Stevenson, R. J., Vaithiyanathan, P., Slate, J., \u0026amp; Richardson, C. J. (2000). Changes in algal assemblages along observed and experimental phosphorus gradients in a subtropical wetland, U.S.A. \u003cem\u003eFreshwater Biology, 44\u003c/em\u003e, 339-353. \u003c/li\u003e\n\u003cli\u003ePatrick, R. (1967). The effect of invasion rate, species pool, and size of area on the structure of the diatom community. \u003cem\u003eProceedings of the Academy of Natural Sciences, 58\u003c/em\u003e, 1335-1342. \u003c/li\u003e\n\u003cli\u003ePatterson, M. E., A.J. Atkinson, B.D. Witcher, K.R.T. Whelan, W.J. Miller, R.J. Waara, . . . Shamblin, R. B. (2008). \u003cem\u003eSouth Florida / Caribbean Network vital signs monitoring plan. \u003c/em\u003e. Retrieved from \u003c/li\u003e\n\u003cli\u003ePillsbury, R., Stevenson, R. J., Munn, M. D., \u0026amp; Waite, I. (2019). Relationships between diatom metrics based on species nutrient traits and agricultural land use. \u003cem\u003eEnvironmental Monitoring and Assessment, 2019\u003c/em\u003e(191). doi:https://doi.org/10.1007/s10661-019-7357-8\u003c/li\u003e\n\u003cli\u003ePotapova, M., \u0026amp; Charles, D. (2002). Benthic diatoms in USA rivers: distributions along spatial and environmental gradients. \u003cem\u003eJournal of Biogeography, 29\u003c/em\u003e, 167-187. \u003c/li\u003e\n\u003cli\u003ePreston, F. W. (1948). The commonness, and rarity, of species. \u003cem\u003eEcology, 29\u003c/em\u003e, 254-283. \u003c/li\u003e\n\u003cli\u003eR Core Team (2021). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL https://www.R-project.org/.\u003c/li\u003e\n\u003cli\u003eRaschke, R. L. (1993). Diatom (Bacillariophyta) community response to phosphorus in the Everglades National Park USA. \u003cem\u003ePhycologia, 32\u003c/em\u003e, 48-58. \u003c/li\u003e\n\u003cli\u003eReavie, E. D., Kireta, A. R., Kingston, J. C., Sgro, G. V., Danz, N. P., Axler, R. P., \u0026amp; Hollenhorst, T. P. (2008). Comparison of simple and multimetric diatom-based indicators for Great Lakes coastline disturbance. \u003cem\u003eJournal of Phycology, 44\u003c/em\u003e, 787-802. \u003c/li\u003e\n\u003cli\u003eRier, S. T., \u0026amp; Stevenson, R. J. (2006). Response of periphytic algae to gradients in nitrogen and phosphorus in streamside mesocosms. \u003cem\u003eHydrobiologia, 561\u003c/em\u003e, 131-147. \u003c/li\u003e\n\u003cli\u003eSlate, J. E., \u0026amp; Stevenson, R. J. (2007). The diatom flora of phosphorus-enriched and unenriched sites in an Everglades marsh. \u003cem\u003eDiatom Research, 22\u003c/em\u003e, 355-386. doi:10.1080/0269249X.2007.9705721\u003c/li\u003e\n\u003cli\u003eSmith, S. (2018) Total Phosphorus (TP) Protocol, Seagrass Ecosystems Research Laboratory (SERL) Coastlines and Oceans Division, Institute of Environment, Florida International University.\u003c/li\u003e\n\u003cli\u003eSolorzano, L. \u0026amp; Sharp, J. H. (1980) Determination of total dissolved phosphorus and particulate phosphorus in natural waters. Limnology and Oceanography 25: 754\u0026ndash;758.\u003c/li\u003e\n\u003cli\u003eSoranno, P. A., Cheruvelil, K. S., Stevenson, R. J., Rollins, S. L., Holden, S. W., Heaton, S., \u0026amp; Torng, E. (2008). A framework for developing ecosystem-specific nutrient criteria: Integrating biological thresholds with predictive modeling. \u003cem\u003eLimnology and Oceanography, 53\u003c/em\u003e, 773-787. \u003c/li\u003e\n\u003cli\u003eStevenson, R. J. (2006). Refining diatom indicators for valued ecological attributes and development of water quality criteria. In N. Ognjanova-Rumenova \u0026amp; K. Manoylov (Eds.), \u003cem\u003eAdvances in Phycological Studies\u003c/em\u003e (pp. 365-383). Moscow, Russia: Pensoft Publishers.\u003c/li\u003e\n\u003cli\u003eStevenson, R. J. (2011). A revised framework for coupled human and natural systems, propagating thresholds, and managing environmental problems. \u003cem\u003ePhysics and Chemistry of the Earth, 36\u003c/em\u003e, 342-351. \u003c/li\u003e\n\u003cli\u003eStevenson, R. J. (2014). Ecological assessment with algae: a review and synthesis. \u003cem\u003eJournal of Phycology, 50\u003c/em\u003e, 437-461. \u003c/li\u003e\n\u003cli\u003eStevenson, R. J., Bailey, R. C., Harass, M. C., Hawkins, C. P., Alba-Tercedor, J., Couch, C., . . . Johnson, R. K. (2004). Designing data collection for ecological assessments. In M. T. Barbour, S. B. Norton, H. R. Preston, \u0026amp; K. W. Thornton (Eds.), \u003cem\u003eEcological Assessment of Aquatic Resources: Linking Science to Decision-Making\u003c/em\u003e (pp. 55-84). Pensacola, Florida: Society of Environmental Toxicology and Contamination Publication.\u003c/li\u003e\n\u003cli\u003eStevenson, R. J., Hill, B. E., A.T. Herlihy, Yuan, L. L., \u0026amp; Norton, S. B. (2008). Algal-P relationships, thresholds, and frequency distributions guide nutrient criterion development. \u003cem\u003eJournal of the North American Benthological Society, 27\u003c/em\u003e, 783-799. \u003c/li\u003e\n\u003cli\u003eStevenson, R. J., Pan, Y., Manoylov, K., Parker, C., Larsen, D. P., \u0026amp; Herlihy, A. T. (2008). Development of diatom indicators of ecological conditions for streams of the western United States. \u003cem\u003eJournal of the North American Benthological Society, 27\u003c/em\u003e, 1000-1016. \u003c/li\u003e\n\u003cli\u003eStevenson, R. J., Pan, Y., \u0026amp; Vaithiyanathan., P. (2002). Ecological assessment and indicator development in wetlands: the case of algae in the Everglades, USA. \u003cem\u003eVerhandlungen Internationale Vereinigung f\u0026uuml;r Theoretische und Andgewandte Limnologie, 28\u003c/em\u003e, 1248-1252. \u003c/li\u003e\n\u003cli\u003eStevenson, R. J., Pan, Y., \u0026amp; van Dam, H. (2010). Assessing ecological conditions in rivers and streams with diatoms. In J. P. Smol \u0026amp; E. F. Stoermer (Eds.), \u003cem\u003eThe Diatoms: Applications to the Environmental and Earth Sciences\u003c/em\u003e (Second ed., pp. 57-85). Cambridge, UK: Cambridge University Press.\u003c/li\u003e\n\u003cli\u003eStevenson, R. J., Zalack, J., \u0026amp; Wolin, J. (2013). A multimetric index of lake diatom condition using surface sediment assemblages. \u003cem\u003eFreshwater Science, 32\u003c/em\u003e, 1005\u0026ndash;1025. \u003c/li\u003e\n\u003cli\u003eStoddard, J. L., Larsen, D. P., Hawkins, C. P., Johnson, R. K., \u0026amp; Norris, R. H. (2006). Setting expectations for the ecological condition of streams: the concept of reference condition. \u003cem\u003eEcological Applications, 16\u003c/em\u003e, 1267-1276. \u003c/li\u003e\n\u003cli\u003eTang, T., Stevenson, R. J., \u0026amp; Grace, J. B. (2020). The importance of natural versus human factors for ecological conditions of streams and rivers. \u003cem\u003eScience of the Total Environment, 704\u003c/em\u003e, 135268. doi:https://doi.org/10.1016/j.scitotenv.2019.135268\u003c/li\u003e\n\u003cli\u003eTang, T., Stevenson, R. J., \u0026amp; Infante, D. M. (2016). Accounting for regional variation in both natural environment and human disturbance to improve performance of multimetric indices of lotic benthic diatoms. \u003cem\u003eScience of the Total Environment 268\u003c/em\u003e, 1124-1134. \u003c/li\u003e\n\u003cli\u003eTherneau, T, \u0026amp; Atkinson, B. (2025). rpart: Recursive Partioning and Regression Trees (Version R package (4.1.24). DOI 10.32614/CRAN.package.rpart \u003c/li\u003e\n\u003cli\u003eUrgelles, R., Whelan, K. R. T., Muxo, R., Shamblin, R. B., Patterson, J. M., \u0026amp; Atkinson, A. J. (2019) Periphyton monitoring in Big Cypress National Preserve: Protocol narrative. Natural Resource Report. NPS/SFCN/NRR\u0026mdash;2019/1911. National Park Service. Fort Collins, Colorado\u003c/li\u003e\n\u003cli\u003eUSEPA. (1992). \u003cem\u003eFramework for Ecological Risk Assessment\u003c/em\u003e (EPA/630/R-92/001) United States Environmental Protection Agency, Washington, DC: \u003c/li\u003e\n\u003cli\u003eUSEPA. (2000). \u003cem\u003eStressor Identification Guidance Document\u003c/em\u003e (EPA 600/R/96/059) United States Environmental Protection Agency, Washington, DC: \u003c/li\u003e\n\u003cli\u003eVan Dam, H., Mertenes, A., \u0026amp; Sinkeldam, J. (1994). A coded checklist and ecological indicator values of freshwater diatoms from the Netherlands. \u003cem\u003eNetherlands Journal of Aquatic Ecology, 28\u003c/em\u003e, 117-133. \u003c/li\u003e\n\u003cli\u003eWang, Y. K., Stevenson, R., Sweets, P., \u0026amp; DiFranco, J. (2006). Developing and Testing Diatom Indicators for Wetlands in the Casco Bay Watershed, Maine, USA. \u003cem\u003eHydrobiologia, 561\u003c/em\u003e(1), 191-206. Retrieved from http://springerlink.metapress.com/openurl.asp?genre=article\u0026amp;id=doi:10.1007/s10750-005-1614-2\u003c/li\u003e\n\u003cli\u003eWang, Y. K., Stevenson, R. J., \u0026amp; Metzmeier, L. (2005). Development and evaluation of a diatom-based index of biotic integrity for the Interior Plateau Ecoregion. \u003cem\u003eJournal of the North American Benthological Society, 24\u003c/em\u003e, 990-1008. \u003c/li\u003e\n\u003cli\u003eWhelan, K. R. T., Prats, M., \u0026amp; Atkinson, A. J. (2020) Western Big Cypress National Preserve vegetation map. Natural Resource Report NPS/SFCN/NRR\u0026mdash;2020/2170. National Park Service, Fort Collins, Colorado. https://doi.org/10.36967/nrr-2278240. \u003c/li\u003e\n\u003cli\u003eWilson, S. (2020) BCAL phosphorus sample preparation protocol., Blue Carbon Analysis Laboratory and Seagrass Ecosystems Research Laboratory, Coastlines and Oceans\u003c/li\u003e\n\u003cli\u003eWorm, B., Lotze, H. K., Hillebrand, H., \u0026amp; Sommer, U. (2002). Consumer versus resource control of species diversity and ecosystem functioning. \u003cem\u003eNature, 417\u003c/em\u003e(6891), 848-851. \u003c/li\u003e\n\u003cli\u003eWyatt, K. H., Stevenson, R. J., \u0026amp; Turetsky, M. R. (2010). The importance of nutrient co-limitation in regulating algal community composition, productivity, and algal-derived DOC in an oligotrophic marsh in interior Alaska \u003cem\u003eFreshwater Biology 55\u003c/em\u003e, 1845-1860. \u003c/li\u003e\n\u003c/ol\u003e"}],"fulltextSource":"","fullText":"","funders":[],"hasAdminPriorityOnWorkflow":false,"hasManuscriptDocX":true,"hasOptedInToPreprint":true,"hasPassedJournalQc":"","hasAnyPriority":false,"hideJournal":false,"highlight":"","institution":"","isAcceptedByJournal":true,"isAuthorSuppliedPdf":false,"isDeskRejected":"","isHiddenFromSearch":false,"isInQc":false,"isInWorkflow":false,"isPdf":false,"isPdfUpToDate":true,"isWithdrawnOrRetracted":false,"journal":{"display":true,"email":"[email protected]","identity":"environmental-monitoring-and-assessment","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":false,"externalIdentity":"emas","sideBox":"Learn more about [Environmental Monitoring and Assessment](http://link.springer.com/journal/10661)","snPcode":"10661","submissionUrl":"https://submission.nature.com/new-submission/10661/3","title":"Environmental Monitoring and Assessment","twitterHandle":"","acdcEnabled":true,"dfaEnabled":true,"editorialSystem":"em","reportingPortfolio":"Springer Hybrid","inReviewEnabled":true,"inReviewRevisionsEnabled":false},"keywords":"algae, ecology, management, metrics, periphyton, wetlands, phosphorous","lastPublishedDoi":"10.21203/rs.3.rs-6959033/v1","lastPublishedDoiUrl":"https://doi.org/10.21203/rs.3.rs-6959033/v1","license":{"name":"CC BY 4.0","url":"https://creativecommons.org/licenses/by/4.0/"},"manuscriptAbstract":"\u003cp\u003eWe used relationships between diatom species composition and phosphorus concentration in periphyton (mat P) to develop a monitoring and assessment program for wetlands in the Big Cypress National Preserve (BCNP) in Florida, USA. Cluster analysis and regression showed limitation of taxonomic composition of assemblages to a few low P taxa in low mat P conditions, with additional species being able to invade habitats that have higher mat P concentrations. TITAN (Threshold Indicator Taxa ANalysis) and regression respectively identified 11 and 8 low P taxa and both methods identified 47 high P taxa. Congruence of our study results and results of phosphorus experiments confirmed P caused species responses. Metrics using the low and high P taxa traits were highly related to mat P. Metrics calculated with TITAN derived taxa traits and a novel log transformation of relative abundances were most highly related to mat P, differed most among regions with different human disturbance, and were unaffected by natural factors. Benchmarks for management targets were observed for minimally disturbed conditions and for greater than usual changes in assemblages along the P gradient. Diatom metrics were more highly correlated with distance from P sources than mat P, indicating species-based metrics have high value for monitoring and assessment. Our diatom metric development methods, benchmarks in ecological response and minimally disturbed condition, and causal analysis provide multiple new findings integrated for a monitoring and assessment program with effects-based management targets in BCNP. In addition, we advance important examples for applications in other ecological settings.\u003c/p\u003e","manuscriptTitle":"Developing a monitoring and assessment program with diatoms, an improved metric calculation method, and causal analysis for Big Cypress National Preserve, Florida (USA)","msid":"","msnumber":"","nonDraftVersions":[{"code":1,"date":"2025-07-11 15:07:19","doi":"10.21203/rs.3.rs-6959033/v1","editorialEvents":[{"type":"communityComments","content":0},{"type":"decision","content":"Revision requested","date":"2025-08-02T08:20:19+00:00","index":"","fulltext":""},{"type":"editorInvitedReview","content":"","date":"2025-07-13T19:11:31+00:00","index":"hide","fulltext":""},{"type":"reviewerAgreed","content":"203697528487116441559128245335253011624","date":"2025-07-08T15:11:48+00:00","index":"hide","fulltext":""},{"type":"reviewerAgreed","content":"80232767737052398661906096416003307966","date":"2025-07-08T08:28:14+00:00","index":"hide","fulltext":""},{"type":"reviewersInvited","content":"","date":"2025-07-07T16:30:24+00:00","index":"","fulltext":""},{"type":"editorAssigned","content":"","date":"2025-07-02T00:25:38+00:00","index":"","fulltext":""},{"type":"checksComplete","content":"","date":"2025-07-02T00:25:29+00:00","index":"","fulltext":""},{"type":"submitted","content":"Environmental Monitoring and Assessment","date":"2025-06-23T17:30:47+00:00","index":"","fulltext":""}],"status":"published","journal":{"display":true,"email":"[email protected]","identity":"environmental-monitoring-and-assessment","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":false,"externalIdentity":"emas","sideBox":"Learn more about [Environmental Monitoring and Assessment](http://link.springer.com/journal/10661)","snPcode":"10661","submissionUrl":"https://submission.nature.com/new-submission/10661/3","title":"Environmental Monitoring and Assessment","twitterHandle":"","acdcEnabled":true,"dfaEnabled":true,"editorialSystem":"em","reportingPortfolio":"Springer Hybrid","inReviewEnabled":true,"inReviewRevisionsEnabled":false}}],"origin":"","ownerIdentity":"6c0b2a07-e8bf-4970-9fb4-9a7ad2766dbf","owner":[],"postedDate":"July 11th, 2025","published":true,"recentEditorialEvents":[],"rejectedJournal":[],"revision":"","amendment":"","status":"published-in-journal","subjectAreas":[],"tags":[],"updatedAt":"2026-03-16T16:04:42+00:00","versionOfRecord":{"articleIdentity":"rs-6959033","link":"https://doi.org/10.1007/s10661-026-15144-0","journal":{"identity":"environmental-monitoring-and-assessment","isVorOnly":false,"title":"Environmental Monitoring and Assessment"},"publishedOn":"2026-03-15 15:58:41","publishedOnDateReadable":"March 15th, 2026"},"versionCreatedAt":"2025-07-11 15:07:19","video":"","vorDoi":"10.1007/s10661-026-15144-0","vorDoiUrl":"https://doi.org/10.1007/s10661-026-15144-0","workflowStages":[]},"version":"v1","identity":"rs-6959033","journalConfig":"researchsquare"},"__N_SSP":true},"page":"/article/[identity]/[[...version]]","query":{"redirect":"/article/rs-6959033","identity":"rs-6959033","version":["v1"]},"buildId":"8U1c8b4HqxoKbykW_rLl7","isFallback":false,"isExperimentalCompile":false,"dynamicIds":[84888],"gssp":true,"scriptLoader":[]}