Pink shrimp Farfantepenaeus duorarum spatiotemporal abundance trends along an urban, subtropical shoreline slated for restoration

The Biscayne Bay Coastal Wetlands (BBCW) project of the Comprehensive Everglades Restoration Plan (CERP) aims to reduce point-source freshwater discharges and spread freshwater flow along the mainland shoreline of southern Biscayne Bay. These actions will be taken to approximate conditions in the coastal wetlands and bay that existed prior to construction of canals and water control structures. An increase in pink shrimp (Farfantepenaeus duorarum) density to ≥ 2 individuals m-2 during the wet season (i.e., August-October) along the mainland shoreline was previously proposed as an indication of BBCW success. This study examined pre-BBCW baseline densities and compared them with the proposed target. Densities were monitored by seasonal (wet, dry) throw-trapping (1 m2 replicated in triplicate) at 47 sites along ~22 km of the southwestern Biscayne Bay coastline over 10 years (2007–2016). Densities varied across years and were most often higher in dry seasons. Quantile regression revealed density limitation by four habitat attributes: water temperature (°C), depth (m), salinity (ppt), and submerged aquatic vegetation (SAV: % cover). Procrustean analyses that tested for concordance between the spatial and temporal distributions of shrimp densities and habitat metrics found that water temperature, water depth, and salinity explained ~ 28%, 28%, and 22% of density variability, respectively. No significant relationship with SAV was observed. Hierarchical clustering was used to identify spatially and temporally similar groupings of pink shrimp densities by sites or season-years. Significant groupings were then investigated with respect to potentially limiting habitat attributes. Six site and four year-season clusters were identified. Although habitat attributes significantly differed among spatial clusters, within-cluster median pink shrimp densities did not correlate with within-cluster minima, maxima, medians, or standard deviations of habitat attributes. Overall, pink shrimp density (X¯ = 0.86, SD = 1.32 shrimp m-2) was significantly lower (t(α = 0.10,2),939 = -26.53, P <0.0001) than the 2 shrimp m-2 CERP Interim Goal target. Pink shrimp density corresponded significantly with salinity and appeared limited to density < 2 shrimp m-2 by salinity < ~18 ppt. Salinity is an environmental attribute that will be directly influenced by CERP implementation.

). Cluster 6 (red) exhibited similar spatial density patterns that were, on average, significantly higher than all other cluster groups. Conversely, Cluster 1 (light blue) exhibited significant, consistently lower densities than all other groups except for Cluster 4 (dark green). Clusters 2-5 exhibited intermediate average densities (Fig 4). anticipated to enrich estuarine faunal assemblages as well as increase estuarine species distributions and abundances [4,27,28]. Expansion of continuous submerged aquatic vegetation (SAV) habitats dominated by Halodule wrightii, a species commonly associated with low and variable salinity, is foreseen [11,[24][25][26]29,30]. BBCW implementation goals for benthic habitat include increased spatial extent of nearshore seagrass beds, especially seaward expansion of H. wrightii [30]. Increased overlap of optimal salinity conditions with preferred benthic SAV habitats would provide synergistic benefits to estuarine fauna such as pink shrimp [4,5,27,31].
The purpose of this study was to investigate spatiotemporal trends in pink shrimp density along the southwestern Biscayne Bay shoreline. We investigated the plausibility of the post-CERP establishment of � 2 shrimp m -2 IG. Further, we addressed presumptions that (1) pink shrimp peak abundance occurs during the wet season and (2) nearshore mesohaline salinity goals would yield increased pink shrimp abundance in the nearshore zone. Pink shrimp density relationships to species-specific and total SAV % cover and canopy height were also investigated. Our focus was on evaluating temporal (i.e., seasonal and inter-annual) and spatial (i.e., among sites) pink shrimp density trends relative to habitat attributes. This was achieved by (1) using quantile regression to identify habitat attributes that potentially limit pink shrimp density, (2) organizing pink shrimp density and habitat observations via heatmaps to visually assess spatiotemporal variability and trends, (3) using Procrustean analysis to measure concordance between shrimp density and habitat attribute matrices, (4) employing hierarchical clustering analysis to identify spatiotemporal density clusters, and (5) investigating distributional aspects (median, minimum, maximum, and standard deviation) of habitat attribute values (temperature, salinity, water depth, and SAV % cover) within density clusters to link density patterns to the environment. These analyses employ data from wet and dry seasons of 10 years of epifaunal community sampling at 47 sites within 50 m of shore spanning~22 km of shoreline.

Study area
Biscayne Bay is a large (1,110 km 2 ), shallow (depths generally < 3 m), subtropical lagoon system extending approximately 56 km north to south along the southeast coast of Florida, USA (Fig 1). Where coastal urban development is low, the bay's mainland shoreline consists of mangrove-seagrass ecotone punctuated by natural tidal creeks, freshwater canals, and other channels [32]. Overland freshwater discharges and groundwater seepage create a salinity gradient across the bay perpendicular to the shoreline with three salinity zones: (1) western nearshore areas usually affording the lowest salinities; (2) the bay central axis marked by near oceanic salinities; and (3) oceanic salinities near the eastern passes through barrier islands to the open ocean [24,29,33]. Tidal ranges are generally 0.5 to 1 m [34,35].

Field surveys
Epifaunal communities and SAV habitats were surveyed seasonally at fixed sampling sites (n = 47) along the southwestern Biscayne Bay nearshore zone (0-50 m) from Shoal Point to Turkey Point (Fig 1). Surveys were conducted in public waters under authority of Biscayne National Park (Study #: BISC 06016, Permit #: BISC-2017-SCI0022). Surveys were conducted in dry (January-March sampling) and wet (July-September sampling) seasons that characterize south Florida's climate. The primary sampling unit was the 20 m buffer around GPS coordinates that identified permanent sampling sites. Sites were located in the shallow, open water along the western shoreline mangrove-seagrass ecotone, an area likely to be directly affected by CERP implementation. During each survey, the 47 fixed sampling sites were visited within 3 hr of high tide over 5 to 7 days within a 2 week period. Water quality and habitat parameters, including water temperature (˚C), salinity (ppt), pH, dissolved oxygen saturation (%), dissolved oxygen concentration (mg L -1 ), water depth (m), and sediment depth (m), were recorded at each site. This study analyzed farfantepenaid and associated chemical and physical sampling data collected from 2007 through 2016.
Benthic habitats were assessed for species-specific SAV % cover by visual assessment of 10 replicate 0.5 m 2 quadrats per site [24,25]. In addition, canopy height (maximum seagrass blade length) was measured to provide a topography metric. Species-specific and total SAV % cover data following the methods of Lirman et al. [25] were obtained for the period 2008 to 2016.
Epifaunal communities were sub-sampled at each site (n = 3) using an open-ended, rigidsided aluminum box (i.e., throw-trap) measuring 45 cm tall by 1 m 2 [36,37]. Two 3-mm stretch-mesh cover nets affixed to opposite sides of the throw-trap upper surface prevented epifauna escape during deployment in water deeper than 45 cm. Once deployed, the throwtrap was cleared of trapped epifauna by sweeping (n = 4) its interior from alternating directions with a metal-framed seine fitted with 3 mm stretch-mesh, while gently tapping the substrate with the seine frame. Organisms collected from each sub-sample throw-trap deployment were bagged and numbered separately for storing and processing. Samples were frozen during storage until processing. No protected species were sampled.

Epifauna identification and measurement
Taxonomic identifications and size measurements were conducted in the laboratory. Organisms collected from each replicate throw-trap deployment at a site were maintained and processed independently. Where possible, carapace length (CL, mm) and total length (TL, mm) were recorded for each farfantepenaeid shrimp. Shrimps >8.0 mm CL were identified to species primarily using petasma and thelycum (i.e., sexual) morphology, although other characteristics were also used [38][39][40][41]. Shrimps <8.0 mm CL were identified to genus due to low degree of sexual morphological development [39].

Statistical analysis
All statistical analyses were performed using the R statistical package (The R Foundation, https://www.r-project.org/). Statistical analyses were performed with a Type 1 error criterion of α = 0.10 to reduce potential Type 2 errors. Shrimp density (#/ m 2 ) at each site was calculated as the sum of observed shrimps from the triplicate throw-trap sub-samples divided by 3 m 2 . Density data were natural logarithm (x + 1) transformed before analysis to reduce influence of outliers.
Potential habitat limitations on pink shrimp density. As a statistical interpretation of the ecological concept of Leibig's Law of the Minimum [42,43], quantile regression (QR) has been applied to identify species distribution or abundance limitation by specific habitat attributes by focusing specifically on the upper bound of the abundance vs. habitat attribute relationship [44][45][46][47]. Pink shrimp density was first plotted against individual habitat factors to graphically assess potential limiting factors. QRs (function 'rq' of package 'quantreg') fit to the 0.5 and 0.9 density percentiles were used to statistically identify a subset of habitat attributes that suggested limitation at the median and upper edge of the density distribution. Analyses initially considered water temperature (˚C), salinity (ppt), pH, dissolved oxygen saturation (%), dissolved oxygen concentration (mg/L), water depth (m), sediment depth (m), and the following SAV metrics: Thalassia testudinum % cover, H. wrightii % cover, total seagrass % cover, total SAV % cover, and total SAV canopy height. As in previous studies in the same region [24][25][26], Syringodium filiforme was rarely encountered (n = 10, 1.6% of total samples) and thus was not further considered. Ultimately, potential limiting variables were confined to the following four: temperature, salinity, water depth, and SAV % cover. Due to lack of SAV data for 2007, specific analyses that included SAV covered only the period from 2008 forward.
Spatiotemporal relationships. Heatmaps were generated to visualize spatiotemporal trends in pink shrimp density and the habitat attributes found by QR to potentially limit pink shrimp density. Observation data were converted to 47 row by 20 column matrices to display their spatial (47 sampling sites) and temporal (10 year x 2 seasons = 20 year-seasons) patterns, and color gradients were used to represent the magnitude of density.
Procrustean analyses allowed direct testing of statistical concordance between matrices of shrimp density and habitat attributes [51][52][53]. Procrustean analysis minimizes the residual sum of squares between a target matrix (X: here, shrimp density) and a second, fitted matrix (Y: here, habitat attributes), superimposed on it by scaling, rotating, and dilating [51,52]. The Procrustean Sum of Squares (PSS, also known as Gower's Statistic: m 2 X,Y ) represents the minimized residual sum of squares from the fitting procedure and is used to assess Procrustean fit ranging from 0 to 1, with higher values presenting poorer fit [51][52][53]. The PSS metric is equivalent to 1 -r 2 , where r is a Pearson correlation coefficient [52]. Because the method hinges on one-to-one relationships between the matrices being compared, Procrustean analysis [54,55] cannot accommodate missing values. Following Adams et al. [54] and Arbour and Brown [55], missing habitat attribute values were imputed with linear regressions that included site, season, and year as potential factors. PROTEST (function 'protest,' package 'vegan', permutation n = 9999) provided statistical significance of Procrustean fits between density and habitat attribute matrices [51].
Pink shrimp density and habitat attributes among density clusters. Hierarchical clustering procedures were used to identify groups with similar density patterns among sites or year-seasons [56,57]. Bray-Curtis dissimilarity matrices were constructed ('vegdist' function, 'vegan' package) with respect to site (i.e., spatial) and year-season (i.e., temporal) density observations. Hierarchical agglomerative clustering (function 'hclust') using the "Ward.D2" agglomeration method identified spatially and temporally similar density groupings. A prioiri statistical significance of clusters was tested via similarity profiling (function 'simprof' of package 'clustsig') (permutations = 999, number of expected groups = 1000) of identified density cluster memberships [57]. Permutational multivariate ANOVA (PERMANOVA: function 'adonis2' of package 'vegan') testing provided a posteriori cluster significance [58,59]. PERMA-NOVA was also used to investigate inter-annual and seasonal density differences using yearseason cluster membership as a categorical nesting factor. To investigate potential dispersion influences on PERMANOVA significance, multivariate homogeneity of dispersions analysis (function 'betadisper' of package 'vegan') was used to test for inter-cluster differences in dispersion (i.e., distance to centroid) [60]. The density heat map was rearranged to display site and year-season cluster membership.
Pink shrimp density and habitat attributes previously detected via QR as potentially limiting to pink shrimp density were investigated among site and year-season clusters. First, medians and confidence intervals (CI) of density and habitat attributes were computed for each site cluster and year-season cluster. Confidence intervals (CIs) about median values were computed as: where IQR = interquantile ranges and n = sample size, as described in McGill et al. [61] and Chambers et al. [62]. Plots of density and habitat attributes' median, CIs, minimum, and maximum values were used to visualize their distributions within site and year-season clusters. Density and habitat attributes were analyzed with respect to site or year-season clusters. Nonparametric tests were used because parametric normality and equality of variance assumptions were usually violated. Kruskal-Wallis tests were used to investigate differences in distribution shape and range (i.e., location: [63]) of density and habitat attributes among site or year-season clusters. Post-hoc Tukey-type nonparametric Conover multiple comparison tests (function 'posthoc.tukey.conover.test' of package 'PMCMR') were used to test for significant pairwise differences. These tests were implemented as χ 2 distributions to correct for data ties, and p-values were Bonferroni-corrected [63,64]. A series of correlation analyses was used to identify habitat attribute distribution characteristics that associated with site or year-season cluster median densities. Pearson correlation analyses were applied to median, minimum, maximum, and standard deviation of habitat attributes within site or year-season clusters.
Temperatures ranged from 12.49 to 36.06˚C. Average temperatures demonstrated a clear pattern of cooler (22.73 ± 2.65˚C) and warmer (30.65 ± 1.53˚C) values for dry and wet seasons, respectively ( Table 1). The dry season record was punctuated by an extreme cold front event that occurred during the 2010 field sampling. No pattern of variation in average temperatures among sites was readily discernable (Table 2). Salinities ranged from 2.48 to 39.71 ppt; overall average salinity was 23.64 (± 7.05) ppt (Table 2). Spatially averaged wet season salinities were generally lower than those of dry seasons, although 2011, 2014, and 2015 wet seasons were notable exceptions with higher average salinity than both the preceding and following dry seasons (Table 1). Sampling sites' mean salinity and standard deviation of salinity were negatively correlated (Pearson r = -0.63, t = -5.49, d.f. = 45, p < 0.0001; S2 Fig). Prior to epifaunal sampling during the 2011 and 2015 wet seasons, hypersaline (>40 ppt) conditions of considerable duration were recorded by continuous salinity loggers deployed within this sampling domain [65]. Only four temporally averaged site salinities were mesohaline, most (n = 42) were polyhaline, and one was euhaline ( Table 2). Water depths ranged from 0.19 to 1.5 m and averaged 0.74 (± 0.20) overall (Table 2) with no appreciable trends among year-seasons or among sites. Total SAV % cover ranged from 4.57 to 100% and averaged 66.57% (± 21.97) with no clear year-season variation patterns (Tablea 1, 2). A planktonic microalgal bloom event was observed in parts of the Biscayne Bay coastal area during the 2013 wet season [65,66].

Habitat limitations on pink shrimp density
Of the multiple habitat attributes investigated, significant QR analysis results revealed that temperature (˚C), salinity (ppt), water depth (m), and SAV (% cover) potentially limited pink shrimp density (Table 3, Fig 2). QR of density vs. temperature yielded a single-knot natural cubic spline relationship which was roughly dome-shaped and maximized at 26.6˚C, with tails that tapered off at higher and lower temperatures ( Table 3, Fig 2A). Temperatures between 21.08 and 31.33˚C did not appear to limit pink shrimp densities to <2 shrimp m -2 (Table 3, Fig 2A). Although a series of functional shapes was considered for the QR density vs. salinity response curve, only the linear and log-linear responses were both significant and ecologically plausible [67]. The log-linear response, which suggested more severe density limitation below10 ppt and appeared more asymptotic at salinities above 10 ppt (Table 3, Fig 2B), seemed more plausible than the linear response. Salinities <~18 ppt limited shrimp density to <2 shrimp m -2 (Table 3, Fig 2B). QR of pink shrimp density against water depth (m) yielded a 3 knot (0.25, 0.5, 0.75 quantile) splined bimodal relationship with steep increases in limitation below~0.6 m and above~1.0 m ( Table 3, Fig 2C). Apparent constraint of density to <2 shrimp m -2 occurred at water depths less than 0.43 m and greater than 1.05 m ( Table 3, Fig  2C). Shrimp density had a logarithmic linear relationship with SAV % cover (Table 3, Fig 2D). SAV cover less than 45% limited density to <2 shrimp m -2 ( Table 3, Fig 2D). For the four habitat attributes, significant QRs were observed at the 0.9 quantile, but not at the 0.5 quantile (Table 3).

Spatiotemporal relationships
Heatmap visualization of pink shrimp spatiotemporal density trends revealed a general absence of pink shrimp from sites 13 to 28 (approximately Black Point to Fender Point, Fig 1) and sites 45 to 47 (near Turkey Point, Fig 1) across all year-seasons ( Fig 3A). Within these groups of sites, only 16 (4.4%, N = 360) and 4 (6.7%, N = 60) instances of pink shrimp densities >2 shrimp m -2 were observed, respectively. Generally, higher densities were observed at sites   Fig 3A). Other year-seasons (2008 dry, 2009 dry, 2012 dry, 2012 wet, 2014 dry, and 2015 wet: Fig 3A) exhibited relatively high average density (> 1 shrimp m -2 ), because the low and zero-catch observations from Black Point to Fender Point (Fig 1) were offset by higher density observations at most other sites (Fig 3A). Average densities in these six yearseasons yielded the highest year-season average densities, all >1 shrimp m -2 (Table 1). Cluster analyses were performed to organize year-season elements and site-specific elements into groups according to shrimp density. Then, the heatmap of shrimp density (Fig 3A) was rearranged to reflect both the site and year-season cluster groupings (Fig 3B). The cluster groupings will be discussed in the following section. Heatmaps were developed to visualize spatiotemporal trends in temperature, salinity, water depth, and SAV (Fig 3C, 3D, 3E and 3F). Procrustean analyses revealed significant concordance of the shrimp density matrix (Fig 3A) with water depth, temperature, and salinity habitat attribute matrices (Fig 3C, 3D and 3E;  Table 2. https://doi.org/10.1371/journal.pone.0198539.g002 Pink shrimp spatiotemporal abundance along shoreline prior to restoration Table 4) but not with the SAV matrix ( Fig 3F, Table 4). Water depth and temperature exhibited the highest correlations, followed by salinity (Table 4). Each comparison yielded a high residual sum of squares (high m 2 X,Y values), indicating relatively weak explanatory power of individual habitat attributes (Table 4). Procrustean fitting procedures explained 28.3, 27.1, and 22.1% of the variability in density for water depth, temperature, and salinity, respectively.

Discussion
Analysis of 10 years of monitoring data revealed few instances (11.2%) of pink shrimp densities > 2 shrimp m -2 , the IG for focal Biscayne Bay pink shrimp populations [1]. All but one spatially averaged year-season density and all but a few temporally averaged site densities in the 2007-2016 database were significantly below 2 shrimp m -2 . CERP implementation is expected to result in more favorable salinity conditions for pink shrimp, leading to higher shrimp densities [1,5,27]. Reductions in extreme salinity variability in the southern half of the study area (i.e., Black Point to Convoy Point: Fig 1) could lead to average densities > 2 shrimp m -2 across the entire study spatial domain. However, our results suggest that~10ppt salinity (i.e., low mesohaline to oligohaline: Fig 2B) was a threshold below which pink shrimp densities were severely limited. Above this threshold, the positive response to salinity continued more gradually. The shape of the response may have been affected by the lack of pink shrimp density observations at salinities >39.71 ppt in this study. Limitation of pink shrimp densities at salinities <10 ppt does not support CERP post-restoration IGs of >2 shrimp m -2 with reduction of Biscayne Bay nearshore salinity regimes to oligohaline and low mesohaline conditions. Pink shrimp densities reported here represent an underestimate of true density as a substantial number of shrimps (~25%) were removed from analysis due to concern about gear sampling inefficiency of shrimps <5 mm CL. Spatial pink shrimp density patterns from Black Point to Convoy Point (Fig 1) were represented by membership in one low density site cluster (cluster 1) and one intermediate density site cluster (cluster 2). This zone is strongly influenced by canal discharges [11,[23][24][25][26]33,68]. Both rapid (<60 min to 2 d) and extreme (~25 ppt) salinity reductions can occur along this stretch of coastline [22,24,25,29,69,70]. Such salinity fluctuations can alter fish community assemblages [23] and may affect foraging behavior and survival [22,23]. Pink shrimp may avoid these conditions as they have been reported to migrate to avoid large-volume riverine inflows [71]. Rapid salinity reductions of greater than 20 ppt cause near complete pink shrimp mortality in laboratory settings [72][73][74]. Low and intermediate density site clusters 1 and 2 included sites 45, 46, and 47 (Figs 3 and 4), which are located near Turkey Point (Fig 1), well south of the canal zone (Black Point to Convoy Point) and so were not impacted by low and variable salinity conditions. These sites had moderate minimum salinities (�11.06 ppt), the highest temporally-averaged salinity across all sites (�29.57 ppt: Table 2), and generally high SAV cover (Table 2, Fig 3F). The cause of their low shrimp density is unknown, but does not seem to be related to limitation due to salinity or SAV cover. Site cluster 6, comprised of sites located away from canal mouths (Fig 1), had the highest median density and the second highest minimum salinity (9.62 ppt). Previous field [9,75] and modeling [35,76] studies describe the area represented by the more northern sites of this cluster as an area of relatively high shrimp abundance. These northern sites were situated immediately across the bay from a wide ocean inlet known as the Safety Valve. This inlet has been considered a primary postlarval immigration pathway [35,76] and may have contributed to high densities at northern sites [77,78]. Shrimp cumulative size frequency distributions differed between northern (sites 1-17) and southern (sites 18-47) sampling sites (S1 Fig). Comparison of these size distributions suggested that juvenile shrimps were more abundant in the north. This difference does not necessarily mean that greater recruitment was occurring in the north because differential growth and/or mortality also could have caused the size distribution difference. The inclusion of southerly sites (33,34,37,40,41,43, and 44) within site cluster 6 (Figs 1 and 3B) contradicts the notion of southern recruitment limitation. Some southern sites in high density cluster 6 are located near mangrove creeks that drain more natural watersheds. Other cluster 6 sites, both northern and southern, are located near, but not immediately adjacent to, canals with relatively small freshwater discharges (Military Canal: 1994-2003 annual mean canal output = 21.9 cfs; Cutler Drain C-100: 1994-2003 annual mean output = 46.1 cfs; [33]). However, sites immediately adjacent to the low-volume-discharge canals (i.e., 6 and 32) do not have high shrimp density and clustered with intermediate and low-density sites (clusters 1 and 2, respectively: Figs 1 and 3).
Most year-seasons (75%) were aggregated within one large cluster, indicating a general lack of inter-annual and inter-season variability in Biscayne Bay juvenile pink shrimp densities. However, year-season cluster 2, which had the highest median cluster shrimp density, consisted mainly (60%) of dry season sampling events. Smaller shrimp were associated with the higher dry season densities (S1 Fig), presumably because of higher dry season postlarval recruitment. Observation of higher densities in the dry rather than wet season is contrary to the pink shrimp IG, which focused on improving 'peak' fall (wet) season abundance [1]. The pink shrimp IG was based upon previous findings of a summer/fall (i.e., wet season) peak in abundance [8,9]. Others reported peak juvenile abundances in late fall/early winter (i.e., November/December) [12] or estimated maximal Biscayne Bay juvenile pink shrimp populations occurring in November (i.e., late fall) [75,79]. Differences in sampling gear and spatial domain and the short durations (�2 yr) of the four pervious studies [8,9,75,79] complicate comparisons with the present study. Although of greater duration, the present study's biannual sampling effort may have insufficient resolution to precisely identify the period of peak pink shrimp density or its year to year variation.
A lack of understanding of Biscayne Bay pink shrimp recruitment complicates study of their abundance patterns. The only study on recruitment reported a late fall through early winter (i.e., October through March: [80]) peak in recruitment, but was too short in duration (1 yr) to provide any information on inter-annual trends or consistency. This peak agreed with juvenile abundance studies reporting a late fall/early winter peak [75,79]. A model of pink shrimp postlarval recruitment from Tortugas and Marquesas spawning grounds found that oceanographic processes favored potential recruitment pathways up the Atlantic coast of the Florida Keys during late wet season and early dry season months [81]. Modeling of larval permit Trachinotus falcatus originating from spawning grounds near those of pink shrimp also  Pink shrimp spatiotemporal abundance along shoreline prior to restoration found similar recruitment patterns for the Florida Keys and Biscayne Bay [82]. Oceanographic, coastal, and climatic conditions affect pink shrimp adult reproductive activity [83,84], and larval abundances [85] and interact with behavior to influence early life stage recruitment to nearshore areas [81,[86][87][88][89][90][91][92][93]. Use of pink shrimp and other offshore spawning species as indicators of ecological conditions in their nearshore nursery grounds is complicated by life cycles affected by prior external conditions [5,32]. According to QR analysis, pink shrimp densities were significantly limited by four habitat attributes (water temperature, salinity, water depth, and SAV % cover). Two of these (i.e., salinity regime and SAV % cover) can be influenced by freshwater management within this study domain. Procrustean analysis confirmed the influence of each habitat attribute except SAV % cover. These habitat attributes vary at different time scales; for example, water depth can differ by as much as 1.3 m within 12 hr during extreme tidal cycles, whereas SAV % cover may integrate salinity, nutrients, water clarity, and other influential factors from 6 mo. to 1 yr or more. Responses on different time scales may affect the ability to discern relationships.
The correlation between spatiotemporal pink shrimp density and water depth was the strongest in this study. This was unexpected given the narrow spatial sampling domain along the mangrove-seagrass ecotone and the small range of variation in water depth in the data. Associations between nearshore pink shrimp abundance and depth have been previously reported [75,80,94,95]. Other studies that focused on very nearshore areas (<100 m from shore) also found higher abundances there [8][9][10]12,96]. Recruiting postlarval pink shrimp often concentrate in SAV near the low-tide mark along shorelines [9,12,[97][98][99][100][101][102]. Other pink shrimp habitat investigations [80,94] found multiple habitat attributes (e.g., salinity, salinity standard deviation, standard deviation of turbidity, temperature, median sediment size, dissolved oxygen concentration, water depth, and benthic habitat characteristics) can influence pink shrimp abundance. As re-iterated by Zink et al. [6], Costello et al. [12] stated that ". . .factors other than salinity per se control abundance of the euryhaline juveniles. . ." The present results clearly indicate water depth is one of these factors.
Fluctuating tidal depths and/or reduced detection probability at greater depths likely contributed to the domed shape of the QR relationship [103]. Short-term variations in water depth are caused by astronomical tides and meterological/climatological events [104], which may override any effects of variation in freshwater inflow, especially in a relatively open bay area such as south-central Biscayne Bay. Hydrological alterations can both increase inundation [105] or decrease water depth via sediment accretion [106]. The multiple broad openings connecting Biscayne Bay to the open ocean promote dominance of water depth by astronomical tides [107,108]. Given the exposed nature of that part of Biscayne Bay's coastline, changes in freshwater inflow within the plausible range that may occur along this coastline are unlikely to appreciably alter water depth within the domain. On longer time scales vertical movements of land surface and coastal geomorphology also can influence nearshore water depth and hydroperiod [104]. Sea level rise may also interact with hydrological changes to further alter coastal water levels and marsh inundation [109]. Although these processes are slow, sea level rise effects on water depth could influence nearshore pink shrimp density within the CERP time frame.
In this study, temperature emerged as the second most influential habitat attribute. The 90 th QR dome shaped response curve was shifted towards cooler water temperatures, which reflected the previously discussed higher densities observed during the dry (winter) season. As has already been discussed, seasonality, and thus temperature, relationships to density are dependent upon recruitment patterns. However, it is important to note here that pink shrimp are known to 'overwinter' in estuaries as far north as North Carolina [98] and can occur at higher abundance during winter periods at the southern limit of their range [110].
Along with temperature, salinity was indicated by Kinne [111,112] as a major factor influencing organismal abundance in estuaries. Re-analyses by Zink et al. [6] of data presented by Brusher and Ogren [113] and Minello [114] found increasing pink shrimp abundance with increasing salinity and no statistical difference between polyhaline and mesohaline conditions. It was unexpected to find no limitation of pink shrimp density at high salinities (>35 ppt), but a shortcoming of this study is lack of hypersaline observations (>40ppt). Perhaps the upper limit of salinity observations in our data (i.e., none >39.71 ppt) was not sufficient to reflect a reduction of densities that may occur under hypersaline conditions [5].
The Biscayne Bay pink shrimp IG proposed a pink shrimp preference for seagrasses and presumed that increased % cover of seagrasses would increase the seaward spatial extent of H. wrightii [30] and pink shrimp abundance [1,5,30]. Presently, total SAV QRs yielded the most plausible relationship between pink shrimp density and the benthic habitat metrics investigated. Procrustean analysis test results did not support this. The seemingly weak statistical relationships with either total or species-specific SAV metrics was unexpected. Pink shrimp associations with H. wrightii have been previously reported [12,96,97,115], while other studies have reported high pink shrimp densities associated with total SAV biomass or % cover [10,116,117]. Although one study reports negative impacts of drift and attached algal biomass [118], the positive relationships reported by most studies suggest a stronger relationship between pink shrimp density and either species-specific or total SAV than found in this study.
Several environmental perturbations occurred during the period covered by this study, but only one may have affected pink shrimp density. Variability in climatic conditions led to both wetter and drier than normal wet seasons ( Fig 3D). However, departure from typical salinity regimes did not seem to influence temporal density patterns appreciably. For example, the second highest wet-season pink shrimp density (1.29 ± 1.65 shrimp m -2 : Table 1) coincided with 2012 record rainfall that reduced salinities (3.34 to 22.08 ppt) across the spatial domain ( Fig  3D). Conversely, the highest wet season pink shrimp average density (1.45 ± 2.25 shrimp m -2 : Table 1) occurred during the 2015 wet season, previously denoted as a 'hypersaline' period [65]. Record dry season rainfall during 2016 yielded the lowest average dry season salinity ( Table 1, Fig 3D) while pink shrimp densities that dry season were moderate (0.84 ± 1.18 shrimp m -2 : Table 1). Despite their differing salinity conditions, these year-seasons were assigned to the same shrimp density cluster (Fig 3B). The 2013 wet season clustered separately from the others, suggesting a negative impact of microalgal bloom conditions [65,66] on pink shrimp density. No discernable impact on pink shrimp densities was observed after passage of an extreme cold front in the 2010 dry season.
Due to the field sampling design, the present study results may be applicable only to shallow areas < 100 m from the shoreline. Application of study results to areas further offshore should proceed with caution, if at all, due to potential interaction with other habitat attributes that influence trends in pink shrimp density. The present study was also limited by the apparent low catchability of very recently settled pink shrimp by the throw trap gear suggested by reduced numbers of pink shrimp from 3 to 5 mm CL (S1 Fig). Pink shrimp postlarvae are generally considered settled in their nursery habitat by 3 mm CL [9,12,[86][87][88]. Pink shrimp postlarvae settle in shallow (�1 m), calm water areas along shorelines [12,115], which would suggest they should be readily available to the present field sampling program that samples nearshore waters generally < 1 m deep.
The RECOVER Biscayne Bay IG has set >2 shrimp m -2 as a target wet season pink shrimp density to be achieved with CERP BBCW implementation. But achievement of BBCW and CERP salinity IGs, which include low mesohaline (<10 ppt) and even oligohaline conditions (<5 ppt), may not support increased pink shrimp density. The Biscayne Bay pink shrimp IG may need modification to clarify whether the � 2 shrimp m -2 target refers to all monitoring observations or a seasonal or annual average density across the entire shoreline and to reconsider spatial and seasonal abundance patterns. The present 10 yr analysis of pink shrimp density patterns could be used to refine the RECOVER Biscayne Bay IG.